Page MenuHomec4science
Files F61756924

pod_engine.py
No OneTemporary
Actions

File Metadata

Created
Wed, May 8, 19:01

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 deepcopy as copy
from rrompy.utilities.base.types import Np1D, Tuple, HFEng, sampList
from rrompy.sampling import sampleList
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['PODEngine']
class PODEngine:
"""
POD engine for general matrix orthogonalization.
"""
def __init__(self, HFEngine:HFEng):
self.HFEngine = HFEngine
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
def GS(self, a:Np1D, Q:sampList, n : int = None,
aA:Np1D = None, QA:sampList = 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 RROMPyException(("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 = self.HFEngine.innerProduct(a, Q)
a = a - Q.dot(nu)
if aA is not None:
aA = aA - QA.dot(nu)
r[:-1] = r[:-1] + nu
r[-1] = self.HFEngine.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:sampList,
only_R : bool = False) -> Tuple[sampList, 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 = copy(A)
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:sampList, Q0 : sampList = None,
only_R : bool = False) -> Tuple[sampList, 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).
"""
N = A.shape[1]
B = copy(A)
V = copy(A)
R = np.zeros((N, N), dtype = A.dtype)
if Q0 is None:
Q = sampleList(np.zeros(A.shape, dtype = A.dtype)
+ np.random.randn(*(A.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.HFEngine.norm(a)
alpha = self.HFEngine.innerProduct(a, Q[k])
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 = self.HFEngine.innerProduct(B[J], V[k])
B[J] = B[J] - 2 * np.outer(V[k], vtB)
R[k, J] = self.HFEngine.innerProduct(B[J], Q[k])
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 = self.HFEngine.innerProduct(Q[J], V[k])
Q[J] = Q[J] - 2 * np.outer(V[k], vtQ)
return Q, R

Event Timeline

c4science · Help