pymor.algorithms.rand_la

Module Contents

class pymor.algorithms.rand_la.RandomizedRangeFinder(A, range_product=None, source_product=None, A_adj=None, power_iterations=0, block_size=None, failure_tolerance=1e-15, num_testvecs=20, lambda_min=None, iscomplex=False, qr_method='gram_schmidt', error_estimator='bs18')[source]

Bases: pymor.core.base.BasicObject

Adaptive randomized range approximation of A.

Algorithm 1 in [BS18] is used as the default error estimator.

Given the Operator A, the return value of this method is the VectorArray B with the property

\[\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.

A faster alternative is the leave-one-out error estimator from [ET24] which can be activated by setting error_estimator='loo'. The leave-one-out error estimator estimates the Frobenius norm of the approximation error. Note that while it can be faster and more accurate in practice, this option does only support Euclidian inner products and there are no results on its failure probability so far.

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.

  • power_iterations – Number of power iterations.

  • A_adj – Adjoint Operator to use for power iterations. If None the adjoint is computed using A, source_product and range_product. Set to A for a self for a known self-adjoint operator.

  • failure_tolerance – Maximum failure probability. Only needed for error_estimator=’bs18’.

  • num_testvecs – Number of test vectors. Only needed for error_estimator=’bs18’.

  • lambda_min – The smallest eigenvalue of source_product. Only needed for error_estimator=’bs18’. If None, the smallest eigenvalue is computed using scipy.

  • block_size – Number of basis vectors to add per iteration.

  • iscomplex – If True, the random vectors are chosen complex.

  • qr_method – QR method used for orthogonalization. Either ‘gram_schmidt’ (default) or ‘shifted_chol_qr’.

  • error_estimator – The error estimator used. Either ‘bs18’ (default) or ‘loo’.

Methods

find_range

Find the range of A.
find_range(basis_size=None, tol=None)[source]

Find the range of A.

Parameters:

  • basis_size – Maximum dimension of range approximation.

  • tol – Error tolerance for the algorithm.

Returns:

|VectorArray| which contains the basis, whose span approximates the range of A.

pymor.algorithms.rand_la.randomized_ghep(A, E=None, n=6, power_iterations=0, oversampling=20, single_pass=False, return_evecs=False)[source]

Approximates a few eigenvalues of a Hermitian linear Operator with randomized methods.

Approximates modes eigenvalues w with corresponding eigenvectors v which solve the eigenvalue problem

\[A v_i = w_i v_i\]

or the generalized eigenvalue problem

\[A v_i = w_i E v_i\]

if E is not None.

This method is an implementation of algorithm 6 and 7 in [SLK16].

Parameters:

  • A – The Hermitian linear Operator for which the eigenvalues are to be computed.

  • E – The Hermitian Operator which defines the generalized eigenvalue problem.

  • n – The number of eigenvalues and eigenvectors which are to be computed. Defaults to 6.

  • power_iterations – The number of power iterations to increase the relative weight of the larger singular values. Ignored when single_pass is True.

  • oversampling – The number of samples that are drawn in addition to the desired basis size in the randomized range approximation process.

  • 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 performs two sets of matvecs Ax (default).

  • return_evecs – If True, the eigenvectors are computed and returned. Defaults to False.

Returns:

  • w – A 1D NumPy array which contains the computed eigenvalues.

  • V – A VectorArray which contains the computed eigenvectors.

pymor.algorithms.rand_la.randomized_svd(A, n, source_product=None, range_product=None, power_iterations=0, oversampling=20)[source]

Randomized SVD of an Operator based on [SHvBW21].

Viewing the Operator \(A\) as an \(m\) by \(n\) matrix, this methods computes and returns the randomized generalized singular value decomposition of A according [SHvBW21]:

\[A = U \Sigma V^{-1},\]

where the inner product on the range \(\mathbb{R}^m\) is given by

\[(x, y)_R = x^TRy\]

and the inner product on the source \(\mathbb{R}^n\) is given by

\[(x, y)_S = x^TSy.\]

Note that \(U\) is \(R\)-orthogonal, i.e. \(U^TRU=I\) and \(V\) is \(S\)-orthogonal, i.e. \(V^TSV=I\). In particular, V^{-1}=V^TS.

Parameters:

  • A – The Operator for which the randomized SVD is to be computed.

  • n – The number of eigenvalues and eigenvectors which are to be computed.

  • source_product – Source product Operator \(S\) w.r.t. which the randomized SVD is computed.

  • range_product – Range product Operator \(R\) w.r.t. which the randomized SVD is computed.

  • power_iterations – The number of power iterations to increase the relative weight of the larger singular values.

  • oversampling – The number of samples that are drawn in addition to the desired basis size in the randomized range approximation process.

Returns:

  • U – VectorArray containing the approximated left singular vectors.

  • s – One-dimensional NumPy array of the approximated singular values.

  • V – VectorArray containing the approximated right singular vectors.