Source code for pymor.algorithms.svd_va

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

"""Module for SVD method of operators represented by |VectorArrays|."""

import numpy as np
import scipy.linalg as spla

from pymor.algorithms.gram_schmidt import gram_schmidt
from pymor.core.defaults import defaults
from pymor.core.logger import getLogger
from pymor.operators.interface import Operator
from pymor.vectorarrays.interface import VectorArray


[docs]@defaults('rtol', 'atol', 'l2_err') def method_of_snapshots(A, product=None, modes=None, rtol=1e-7, atol=0., l2_err=0.): """SVD of a |VectorArray| using the method of snapshots. Viewing the |VectorArray| `A` as a `A.dim` x `len(A)` matrix, the return value of this method is the singular value decomposition of `A`, where the inner product on R^(`dim(A)`) is given by `product` and the inner product on R^(`len(A)`) is the Euclidean inner product. .. warning:: The left singular vectors may not be numerically orthonormal for ill-conditioned `A`. Parameters ---------- A The |VectorArray| for which the SVD is to be computed. product Inner product |Operator| w.r.t. which the SVD is computed. modes If not `None`, at most the first `modes` singular values and vectors are returned. rtol Singular values smaller than this value multiplied by the largest singular value are ignored. atol Singular values smaller than this value are ignored. l2_err Do not return more modes than needed to bound the l2-approximation error by this value. I.e. the number of returned modes is at most :: argmin_N { sum_{n=N+1}^{infty} s_n^2 <= l2_err^2 } where `s_n` denotes the n-th singular value. Returns ------- U |VectorArray| of left singular vectors. s One-dimensional |NumPy array| of singular values. Vh |NumPy array| of right singular vectors. """ assert isinstance(A, VectorArray) assert product is None or isinstance(product, Operator) assert modes is None or modes <= len(A) if A.dim == 0 or len(A) == 0: return A.space.empty(), np.array([]), np.zeros((0, len(A))) logger = getLogger('pymor.algorithms.svd_va.method_of_snapshots') with logger.block(f'Computing Gramian ({len(A)} vectors) ...'): B = A.gramian(product) with logger.block('Computing eigenvalue decomposition ...'): eigvals = (None if modes is None or l2_err > 0. else (len(B) - modes, len(B) - 1)) evals, V = spla.eigh(B, overwrite_a=True, turbo=True, eigvals=eigvals) evals = evals[::-1] V = V.T[::-1, :] tol = max(rtol ** 2 * evals[0], atol ** 2) above_tol = np.where(evals >= tol)[0] if len(above_tol) == 0: return A.space.empty(), np.array([]), np.zeros((0, len(A))) last_above_tol = above_tol[-1] errs = np.concatenate((np.cumsum(evals[::-1])[::-1], [0.])) below_err = np.where(errs <= l2_err**2)[0] first_below_err = below_err[0] selected_modes = min(first_below_err, last_above_tol + 1) if modes is not None: selected_modes = min(selected_modes, modes) s = np.sqrt(evals[:selected_modes]) V = V[:selected_modes] Vh = V.conj() with logger.block(f'Computing left-singular vectors ({len(V)} vectors) ...'): U = A.lincomb(V / s[:, np.newaxis]) return U, s, Vh
[docs]@defaults('rtol', 'atol', 'l2_err') def qr_svd(A, product=None, modes=None, rtol=4e-8, atol=0., l2_err=0.): """SVD of a |VectorArray| using Gram-Schmidt orthogonalization. Viewing the |VectorArray| `A` as a `A.dim` x `len(A)` matrix, the return value of this method is the singular value decomposition of `A`, where the inner product on R^(`dim(A)`) is given by `product` and the inner product on R^(`len(A)`) is the Euclidean inner product. Parameters ---------- A The |VectorArray| for which the SVD is to be computed. product Inner product |Operator| w.r.t. which the left singular vectors are computed. modes If not `None`, at most the first `modes` singular values and vectors are returned. rtol Singular values smaller than this value multiplied by the largest singular value are ignored. atol Singular values smaller than this value are ignored. l2_err Do not return more modes than needed to bound the l2-approximation error by this value. I.e. the number of returned modes is at most :: argmin_N { sum_{n=N+1}^{infty} s_n^2 <= l2_err^2 } where `s_n` denotes the n-th singular value. Returns ------- U |VectorArray| of left singular vectors. s One-dimensional |NumPy array| of singular values. Vh |NumPy array| of right singular vectors. """ assert isinstance(A, VectorArray) assert product is None or isinstance(product, Operator) assert modes is None or modes <= len(A) if A.dim == 0 or len(A) == 0: return A.space.empty(), np.array([]), np.zeros((0, len(A))) logger = getLogger('pymor.algorithms.svd_va.qr_svd') with logger.block('Computing QR decomposition ...'): Q, R = gram_schmidt(A, product=product, return_R=True, check=False) with logger.block('Computing SVD of R ...'): U2, s, Vh = spla.svd(R, lapack_driver='gesvd') with logger.block('Choosing the number of modes ...'): tol = max(rtol * s[0], atol) above_tol = np.where(s >= tol)[0] if len(above_tol) == 0: return A.space.empty(), np.array([]), np.zeros((0, len(A))) last_above_tol = above_tol[-1] errs = np.concatenate((np.cumsum(s[::-1] ** 2)[::-1], [0.])) below_err = np.where(errs <= l2_err**2)[0] first_below_err = below_err[0] selected_modes = min(first_below_err, last_above_tol + 1) if modes is not None: selected_modes = min(selected_modes, modes) U2 = U2[:, :selected_modes] s = s[:selected_modes] Vh = Vh[:selected_modes] with logger.block(f'Computing left singular vectors ({selected_modes} modes) ...'): U = Q.lincomb(U2.T) return U, s, Vh