Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/rrompy/linalg/pod_engine.py
c11471199317master

Last Change
Raw File

pod_engine.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
from copy import copy
from rrompy.utilities.base.types import Np1D, Np2D, Tuple
__all__ = ['PODEngine']
class PODEngine:
"""
POD engine for general matrix orthogonalization.
Args/Attributes:
M: Square matrix associated to scalar product.
"""
def __init__(self, M:Np2D):
self.M = M
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
def norm(self, a:Np1D) -> float:
"""Compute norm of a vector."""
return np.sqrt(np.abs(a.conj().T.dot(self.M.dot(a))))
def GS(self, a:Np1D, Q:Np2D, n : int = None,
aA:Np1D = None, QA:Np2D = None) -> Tuple[Np1D, Np1D, Np1D]:
"""
Compute 1 Gram-Schmidt step with given projector.
Args:
a: vector to be projected;
Q: orthogonal projection matrix;
n: number of columns of Q to be considered;
aA: augmented components of vector to be projected;
QA: augmented components of projection matrix.
Returns:
Resulting normalized vector, coefficients of a wrt the updated
basis.
"""
if n is None:
n = Q.shape[1]
if aA is None != QA is None:
raise Exception(("Either both or none of augmented components "
"must be provided."))
r = np.zeros((n + 1,), dtype = a.dtype)
if n > 0:
Q = Q[:, : n]
for j in range(2): # twice is enough!
nu = Q.conj().T.dot(self.M.dot(a))
a = a - Q.dot(nu)
if aA is not None:
aA = aA - QA.dot(nu)
r[:-1] = r[:-1] + nu
r[-1] = self.norm(a)
if np.isclose(np.abs(r[-1]), 0.):
r[-1] = 1.
a = a / r[-1]
if aA is not None:
aA = aA / r[-1]
return a, r, aA
def QRGramSchmidt(self, A:Np2D,
only_R : bool = False) -> Tuple[Np1D, Np1D]:
"""
Compute QR decomposition of a matrix through Gram-Schmidt method.
Args:
A: matrix to be decomposed;
only_R(optional): whether to skip reconstruction of Q; defaults to
False.
Returns:
Resulting orthogonal and upper-triangular factors.
"""
N = A.shape[1]
Q = np.zeros_like(A, dtype = A.dtype)
R = np.zeros((N, N), dtype = A.dtype)
for k in range(N):
Q[:, k], R[: k + 1, k], _ = self.GS(A[:, k], Q, k)
if only_R:
return R
return Q, R
def QRHouseholder(self, A:Np2D, Q0 : Np2D = None,
only_R : bool = False) -> Tuple[Np1D, Np1D]:
"""
Compute QR decomposition of a matrix through Householder method.
Args:
A: matrix to be decomposed;
Q0(optional): initial orthogonal guess for Q; defaults to random;
only_R(optional): whether to skip reconstruction of Q; defaults to
False.
Returns:
Resulting (orthogonal and )upper-triangular factor(s).
"""
B = copy(A)
N = B.shape[1]
V = np.zeros_like(B, dtype = B.dtype)
R = np.zeros((N, N), dtype = B.dtype)
if Q0 is None:
Q = np.zeros_like(B, dtype = B.dtype) + np.random.randn(*(B.shape))
else:
Q = copy(Q0)
for k in range(N):
if Q0 is None:
Q[:, k], _, _ = self.GS(Q[:, k], Q, k)
a = B[:, k]
R[k, k] = self.norm(a)
alpha = Q[:, k].conj().T.dot(self.M.dot(a))
if np.isclose(np.abs(alpha), 0.): s = 1.
else: s = - alpha / np.abs(alpha)
Q[:, k] = s * Q[:, k]
V[:, k], _, _ = self.GS(R[k, k] * Q[:, k] - a, Q, k)
J = np.arange(k + 1, N)
vtB = V[:, k].conj().dot(self.M.dot(B[:, J]))
B[:, J] = B[:, J] - 2 * np.outer(V[:, k], vtB)
R[k, J] = Q[:, k].conj().T.dot(self.M.dot(B[:, J]))
B[:, J] = B[:, J] - np.outer(Q[:, k], R[k, J])
if only_R:
return R
for k in range(N - 1, -1, -1):
J = np.arange(k, N)
vtQ = V[:, k].conj().T.dot(self.M.dot(Q[:, J]))
Q[:, J] = Q[:, J] - 2 * np.outer(V[:, k], vtQ)
return Q, R
c4science · Help