# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
import numpy as np
from scipy.sparse.linalg import eigsh, LinearOperator
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
[docs]@defaults('tol', 'failure_tolerance', 'num_testvecs')
def adaptive_rrf(A, source_product=None, range_product=None, tol=1e-4,
failure_tolerance=1e-15, num_testvecs=20, lambda_min=None, iscomplex=False):
r"""Adaptive randomized range approximation of `A`.
This is an implementation of Algorithm 1 in [BS18]_.
Given the |Operator| `A`, the return value of this method is the |VectorArray|
`B` with the property
.. math::
\Vert A - P_{span(B)} A \Vert \leq tol
with a failure probability smaller than `failure_tolerance`, where the norm denotes the
operator norm. The inner product of the range of `A` is given by `range_product` and
the inner product of the source of `A` is given by `source_product`.
Parameters
----------
A
The |Operator| A.
source_product
Inner product |Operator| of the source of A.
range_product
Inner product |Operator| of the range of A.
tol
Error tolerance for the algorithm.
failure_tolerance
Maximum failure probability.
num_testvecs
Number of test vectors.
lambda_min
The smallest eigenvalue of source_product.
If `None`, the smallest eigenvalue is computed using scipy.
iscomplex
If `True`, the random vectors are chosen complex.
Returns
-------
B
|VectorArray| which contains the basis, whose span approximates the range of A.
"""
assert source_product is None or isinstance(source_product, Operator)
assert range_product is None or isinstance(range_product, Operator)
assert isinstance(A, Operator)
B = A.range.empty()
R = A.source.random(num_testvecs, distribution='normal')
if iscomplex:
R += 1j*A.source.random(num_testvecs, distribution='normal')
if source_product is None:
lambda_min = 1
elif lambda_min is None:
def mv(v):
return source_product.apply(source_product.source.from_numpy(v)).to_numpy()
def mvinv(v):
return source_product.apply_inverse(source_product.range.from_numpy(v)).to_numpy()
L = LinearOperator((source_product.source.dim, source_product.range.dim), matvec=mv)
Linv = LinearOperator((source_product.range.dim, source_product.source.dim), matvec=mvinv)
lambda_min = eigsh(L, sigma=0, which="LM", return_eigenvectors=False, k=1, OPinv=Linv)[0]
testfail = failure_tolerance / min(A.source.dim, A.range.dim)
testlimit = np.sqrt(2. * lambda_min) * erfinv(testfail**(1. / num_testvecs)) * tol
maxnorm = np.inf
M = A.apply(R)
while(maxnorm > testlimit):
basis_length = len(B)
v = A.source.random(distribution='normal')
if iscomplex:
v += 1j*A.source.random(distribution='normal')
B.append(A.apply(v))
gram_schmidt(B, range_product, atol=0, rtol=0, offset=basis_length, copy=False)
M -= B.lincomb(B.inner(M, range_product).T)
maxnorm = np.max(M.norm(range_product))
return B
[docs]@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 [HMT11]_.
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`.
Parameters
----------
A
The |Operator| A.
source_product
Inner product |Operator| of the source of A.
range_product
Inner product |Operator| of the range of A.
q
The number of power iterations.
l
The block size of the normalized power iterations.
iscomplex
If `True`, the random vectors are chosen complex.
Returns
-------
Q
|VectorArray| which contains the basis, whose span approximates the range of A.
"""
assert source_product is None or isinstance(source_product, Operator)
assert range_product is None or isinstance(range_product, Operator)
assert isinstance(A, Operator)
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)
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