# Source code for pymor.algorithms.randrangefinder

# This file is part of the pyMOR project (http://www.pymor.org).

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.interfaces import OperatorInterface

[docs]@defaults('tol', 'failure_tolerance', 'num_testvecs')
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, OperatorInterface)
assert range_product is None or isinstance(range_product, OperatorInterface)
assert isinstance(A, OperatorInterface)

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, OperatorInterface)
assert range_product is None or isinstance(range_product, OperatorInterface)
assert isinstance(A, OperatorInterface)

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):