Unverified Commit a478fa00 authored by Stephan Rave's avatar Stephan Rave Committed by GitHub
Browse files

Merge pull request #1552 from ullmannsven/main

Add randomized linear algebra methods for svd and ghep
parents f5137c40 0bccad1e
Pipeline #141376 failed with stages
in 83 minutes and 54 seconds
......@@ -22,6 +22,9 @@
* support for discrete-time LTI systems, Lyapunov equations and balanced truncation
* dedicated Hankel operator class
* Sven Ullmann, ullmannsven@gmx.de
* randomized algorithms for generalized SVD and generalized Hermitian eigenvalue problem
### pyMOR 2021.2
* Tim Keil, tim.keil@uni-muenster.de
......
......@@ -387,6 +387,23 @@
publisher={IEEE}
}
@article{SHB21,
title = {Randomized algorithms for generalized singular value decomposition with application to sensitivity analysis},
author = {Saibaba AK, Hart J, van BloemenWaanders B.},
journal = {Numerical Linear Algebra with Applications},
volume = {8},
issue = {4},
year = {2021}
}
@article{SJK16,
title = {Randomized algorithms for generalized {H}ermitian eigenvalue problems with application to computing Karhunen Loeve expansion},
author = {Saibaba AK, Lee J. and Kitanidis PK},
journal = {Numerical Linear Algebra with Applications},
year = {2016},
DOI = { 10.1002/nla.2026}
}
@article{S11,
title={Symplectic Gram-Schmidt Algorithm with Re-Orthogonalization},
author={Al-Aidarous, Eman Salem},
......@@ -395,7 +412,7 @@
number={1},
pages={11--20},
year = {2011},
doi = {10.4197/Sci.23-1.2},
doi = {10.4197/Sci.23-1.2}
}
@article{TRLBK14,
......
......@@ -3,12 +3,16 @@
# License: BSD 2-Clause License (https://opensource.org/licenses/BSD-2-Clause)
import numpy as np
import scipy as sp
from scipy.sparse.linalg import eigsh, LinearOperator
from scipy.linalg import lu_factor, lu_solve
from scipy.special import erfinv
from pymor.algorithms.gram_schmidt import gram_schmidt
from pymor.core.defaults import defaults
from pymor.operators.interface import Operator
from pymor.core.logger import getLogger
from pymor.operators.constructions import InverseOperator, IdentityOperator
@defaults('tol', 'failure_tolerance', 'num_testvecs')
......@@ -94,13 +98,13 @@ def adaptive_rrf(A, source_product=None, range_product=None, tol=1e-4,
@defaults('q', 'l')
def rrf(A, source_product=None, range_product=None, q=2, l=8, iscomplex=False):
"""Randomized range approximation of `A`.
This is an implementation of Algorithm 4.4 in :cite:`HMT11`.
def rrf(A, source_product=None, range_product=None, q=2, l=8, return_rand=False, iscomplex=False):
r"""Randomized range approximation of `A`.
Given the |Operator| `A`, the return value of this method is the |VectorArray|
`Q` whose vectors form an approximate orthonomal basis for the range of `A`.
`Q` whose vectors form an approximate orthonormal basis for the range of `A`.
This method is based on algorithm 2 in :cite:`SHB21`.
Parameters
----------
......@@ -114,6 +118,8 @@ def rrf(A, source_product=None, range_product=None, q=2, l=8, iscomplex=False):
The number of power iterations.
l
The block size of the normalized power iterations.
return_rand
If `True`, the randomly sampled |VectorArray| R is returned.
iscomplex
If `True`, the random vectors are chosen complex.
......@@ -121,21 +127,206 @@ def rrf(A, source_product=None, range_product=None, q=2, l=8, iscomplex=False):
-------
Q
|VectorArray| which contains the basis, whose span approximates the range of A.
R
The randomly sampled |VectorArray| (if `return_rand` is `True`).
"""
assert source_product is None or isinstance(source_product, Operator)
assert range_product is None or isinstance(range_product, Operator)
assert isinstance(A, Operator)
if range_product is None:
range_product = IdentityOperator(A.range)
else:
assert isinstance(range_product, Operator)
if source_product is None:
source_product = IdentityOperator(A.source)
else:
assert isinstance(source_product, Operator)
assert 0 <= l <= min(A.source.dim, A.range.dim) and isinstance(l, int)
assert q >= 0 and isinstance(q, int)
R = A.source.random(l, distribution='normal')
if iscomplex:
R += 1j*A.source.random(l, distribution='normal')
Q = A.apply(R)
gram_schmidt(Q, range_product, atol=0, rtol=0, copy=False)
for i in range(q):
Q = A.apply_adjoint(Q)
Q = A.apply_adjoint(range_product.apply(Q))
Q = source_product.apply_inverse(Q)
gram_schmidt(Q, source_product, atol=0, rtol=0, copy=False)
Q = A.apply(Q)
gram_schmidt(Q, range_product, atol=0, rtol=0, copy=False)
return Q
if return_rand:
return Q, R
else:
return Q
@defaults('p', 'q', 'modes')
def random_generalized_svd(A, range_product=None, source_product=None, modes=6, p=20, q=2):
r"""Randomized SVD of an |Operator|.
Viewing the `A` as a `A.dim x len(A)` matrix, the return value
of this method is the randomized singular value decomposition of `A`, where the
inner product on :math:`\mathbb{R}^{\mathtt{A.dim}}` is given by 'range_product' and
the inner product on :math:`\mathbb{R}^{\mathtt{len(A)}}` is given by `source_product`.
.. math::
A = U \Sigma V^H \mathtt{source_product}
This method is based on :cite:`SHB21`.
Parameters
----------
A
The |Operator| for which the randomized SVD is to be computed.
range_product
Range product |Operator| w.r.t which the randomized SVD is computed.
source_product
Source product |Operator| w.r.t which the randomized SVD is computed.
modes
The first `modes` approximated singular values and vectors are returned.
p
If not `0`, adds `p` columns to the randomly sampled matrix (oversampling parameter).
q
If not `0`, performs `q` so-called power iterations to increase the relative weight
of the first singular values.
Returns
-------
U
|VectorArray| of approximated left singular vectors.
s
One-dimensional |NumPy array| of the approximated singular values.
Vh
|VectorArray| of the approximated right singular vectors.
"""
logger = getLogger('pymor.algorithms.rand_la')
assert isinstance(A, Operator)
assert 0 <= modes <= max(A.source.dim, A.range.dim) and isinstance(modes, int)
assert 0 <= p <= max(A.source.dim, A.range.dim) - modes and isinstance(p, int)
assert q >= 0 and isinstance(q, int)
if range_product is None:
range_product = IdentityOperator(A.range)
else:
assert isinstance(range_product, Operator)
assert range_product.source == range_product.range == A.range
if source_product is None:
source_product = IdentityOperator(A.source)
else:
assert isinstance(source_product, Operator)
assert source_product.source == source_product.range == A.source
if A.source.dim == 0 or A.range.dim == 0:
return A.source.empty(), np.array([]), A.range.empty()
Q = rrf(A, source_product=source_product, range_product=range_product, q=q, l=modes+p)
B = A.apply_adjoint(range_product.apply(Q))
Q_B, R_B = gram_schmidt(source_product.apply_inverse(B), product=source_product, return_R=True)
U_b, s, Vh_b = sp.linalg.svd(R_B.T, full_matrices=False)
with logger.block(f'Computing generalized left-singular vectors ({modes} vectors) ...'):
U = Q.lincomb(U_b.T)
with logger.block(f'Computung generalized right-singular vector ({modes} vectors) ...'):
Vh = Q_B.lincomb(Vh_b)
return U[:modes], s[:modes], Vh[:modes]
@defaults('modes', 'p', 'q')
def random_ghep(A, E=None, modes=6, p=20, q=2, single_pass=False):
r"""Approximates a few eigenvalues of a symmetric linear |Operator| with randomized methods.
Approximates `modes` eigenvalues `w` with corresponding eigenvectors `v` which solve
the eigenvalue problem
.. math::
A v_i = w_i v_i
or the generalized eigenvalue problem
.. math::
A v_i = w_i E v_i
if `E` is not `None`.
This method is an implementation of algorithm 6 and 7 in :cite:`SJK16`.
Parameters
----------
A
The Hermitian linear |Operator| for which the eigenvalues are to be computed.
E
The Hermitian |Operator| which defines the generalized eigenvalue problem.
modes
The number of eigenvalues and eigenvectors which are to be computed.
p
If not `0`, adds `p` columns to the randomly sampled matrix in the :func:`rrf` method
(oversampling parameter).
q
If not `0`, performs `q` power iterations to increase the relative weight
of the larger singular values. Ignored when `single_pass` is `True`.
single_pass
If `True`, computes the GHEP where only one set of matvecs Ax is required, but at the
expense of lower numerical accuracy.
If `False`, the methods require two sets of matvecs Ax.
Returns
-------
w
A 1D |NumPy array| which contains the computed eigenvalues.
V
A |VectorArray| which contains the computed eigenvectors.
"""
logger = getLogger('pymor.algorithms.rand_la')
assert isinstance(A, Operator) and A.linear
assert not A.parametric
assert A.source == A.range
assert 0 <= modes <= max(A.source.dim, A.range.dim) and isinstance(modes, int)
assert 0 <= p <= max(A.source.dim, A.range.dim) - modes and isinstance(p, int)
assert q >= 0 and isinstance(q, int)
if E is None:
E = IdentityOperator(A.source)
else:
assert isinstance(E, Operator) and E.linear
assert not E.parametric
assert E.source == E.range
assert E.source == A.source
if A.source.dim == 0 or A.range.dim == 0:
return A.source.empty(), np.array([]), A.range.empty()
if single_pass:
Omega = A.source.random(modes+p, distribution='normal')
Y_bar = A.apply(Omega)
Y = E.apply_inverse(Y_bar)
Q, R = gram_schmidt(Y, product=E, return_R=True)
X = E.apply2(Omega, Q)
X_lu = lu_factor(X)
T = lu_solve(X_lu, lu_solve(X_lu, Omega.inner(Y_bar)).T).T
else:
C = InverseOperator(E) @ A
Y, Omega = rrf(C, q=q, l=modes+p, return_rand=True)
Q = gram_schmidt(Y, product=E)
T = A.apply2(Q, Q)
w, S = sp.linalg.eigh(T)
w = w[::-1]
S = S[:, ::-1]
with logger.block(f'Computing eigenvectors ({modes} vectors) ...'):
V = Q.lincomb(S)
return w[:modes], V[:modes]
# This file is part of the pyMOR project (https://www.pymor.org).
# Copyright pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (https://opensource.org/licenses/BSD-2-Clause)
import numpy as np
import scipy as sp
from numpy.random import uniform
from pymor.algorithms.rand_la import rrf, adaptive_rrf, random_ghep, random_generalized_svd
from pymor.operators.numpy import NumpyMatrixOperator
from pymor.operators.constructions import VectorArrayOperator
def test_adaptive_rrf():
np.random.seed(0)
A = uniform(low=-1.0, high=1.0, size=(100, 100))
A = A @ A.T
range_product = NumpyMatrixOperator(A)
np.random.seed(1)
B = uniform(low=-1.0, high=1.0, size=(10, 10))
B = B.dot(B.T)
source_product = NumpyMatrixOperator(B)
C = range_product.range.random(10, seed=10)
op = VectorArrayOperator(C)
D = range_product.range.random(10, seed=11)+1j*range_product.range.random(10, seed=12)
op_complex = VectorArrayOperator(D)
Q1 = adaptive_rrf(op, source_product, range_product)
assert Q1 in op.range
Q2 = adaptive_rrf(op_complex, iscomplex=True)
assert np.iscomplexobj(Q2.to_numpy())
assert Q2 in op.range
def test_rrf():
np.random.seed(2)
A = uniform(low=-1.0, high=1.0, size=(100, 100))
A = A @ A.T
range_product = NumpyMatrixOperator(A)
np.random.seed(3)
B = uniform(low=-1.0, high=1.0, size=(10, 10))
B = B @ B.T
source_product = NumpyMatrixOperator(B)
C = range_product.range.random(10, seed=10)
op = VectorArrayOperator(C)
D = range_product.range.random(10, seed=11)+1j*range_product.range.random(10, seed=12)
op_complex = VectorArrayOperator(D)
Q1 = rrf(op, source_product, range_product)
assert Q1 in op.range
assert len(Q1) == 8
Q2 = rrf(op_complex, iscomplex=True)
assert np.iscomplexobj(Q2.to_numpy())
assert Q2 in op.range
assert len(Q2) == 8
def test_random_generalized_svd():
np.random.seed(4)
E = uniform(low=-1.0, high=1.0, size=(5, 5))
E_op = NumpyMatrixOperator(E)
modes = 3
U, s, Vh = random_generalized_svd(E_op, modes=modes, p=1)
U_real, s_real, Vh_real = sp.linalg.svd(E)
assert abs(np.linalg.norm(s-s_real[:modes])) <= 1e-2
assert len(U) == modes
assert len(Vh) == modes
assert len(s) == modes
assert U in E_op.range
assert Vh in E_op.source
def test_random_ghep():
np.random.seed(5)
D = uniform(low=-1.0, high=1.0, size=(5, 5))
D = D @ D.T
D_op = NumpyMatrixOperator(D)
modes = 3
w1, V1 = random_ghep(D_op, modes=modes, p=1, single_pass=False)
w2, V2 = random_ghep(D_op, modes=modes, p=1, single_pass=True)
w_real, V_real = sp.linalg.eigh(D)
w_real = w_real[::-1]
V_real = V_real[:, ::-1]
assert abs(np.linalg.norm(w1-w_real[:modes])) <= 1e-2
assert abs(np.linalg.norm(w2-w_real[:modes])) <= 1
for i in range(0, modes):
assert np.linalg.norm(abs(V1.to_numpy()[i, :]) - abs(V_real[:, i])) <= 1
for i in range(0, modes):
assert np.linalg.norm(abs(V2.to_numpy()[i, :]) - abs(V_real[:, i])) <= 1
assert len(w1) == modes
assert len(V1) == modes
assert V1.dim == D_op.source.dim
# This file is part of the pyMOR project (https://www.pymor.org).
# Copyright pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (https://opensource.org/licenses/BSD-2-Clause)
import numpy as np
from numpy.random import uniform
from pymor.algorithms.randrangefinder import rrf, adaptive_rrf
from pymor.operators.numpy import NumpyMatrixOperator
from pymor.operators.constructions import VectorArrayOperator
np.random.seed(0)
A = uniform(low=-1.0, high=1.0, size=(100, 100))
A = A.dot(A.T)
range_product = NumpyMatrixOperator(A)
np.random.seed(1)
A = uniform(low=-1.0, high=1.0, size=(10, 10))
A = A.dot(A.T)
source_product = NumpyMatrixOperator(A)
B = range_product.range.random(10, seed=10)
op = VectorArrayOperator(B)
C = range_product.range.random(10, seed=11)+1j*range_product.range.random(10, seed=12)
op_complex = VectorArrayOperator(C)
def test_rrf():
Q = rrf(op, source_product, range_product)
assert Q in op.range
assert len(Q) == 8
Q = rrf(op_complex, iscomplex=True)
assert np.iscomplexobj(Q.to_numpy())
assert Q in op.range
assert len(Q) == 8
def test_adaptive_rrf():
B = adaptive_rrf(op, source_product, range_product)
assert B in op.range
B = adaptive_rrf(op_complex, iscomplex=True)
assert np.iscomplexobj(B.to_numpy())
assert B in op.range
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment