adding randomized algorithms

AUTHORS.md

@@ -14,6 +14,7 @@
     * improve projection shifts for low-rank ADI method for Lyapunov equations
 
 * Dennis Eickhorn, d.eickhorn@uni-muenster.de
+    * randomized range approximation algorithms in algorithms.randrangefinder
     * fixed complex norms in vectorarrays.interfaces
 
 

docs/source/bibliography.rst

@@ -33,6 +33,10 @@ Bibliography
             Equations, Proceedings of the 11th World Congress on Computational
             Mechanics, 2014.
 
+.. [BS18] A. Buhr, K. Smetana, 
+          Randomized Local Model Order Reduction. 
+          SIAM Journal on Scientific Computing, 40(4), A2120–A2151, 2018.
+
 .. [CLVV06] Y. Chahlaoui, D. Lemonnier, A. Vandendorpe, P. Van
             Dooren,
             Second-order balanced truncation,
@@ -59,6 +63,12 @@ Bibliography
            approximations of parametrized evolution equations,
            M2AN 42(2), 277-302, 2008.
 
+.. [HMT11] N. Halko, P. G. Martinsson and J. A. Tropp,
+           Finding structure with randomness: probabilistic
+           algorithms for constructing approximate matrix
+           decompositions, 
+           SIAM Review, 53(2), 217–288, 2011.
+
 .. [HLR18] C. Himpe, T. Leibner, S. Rave,
            Hierarchical Approximate Proper Orthogonal Decomposition,
            SIAM J. Sci. Comput. 40, A3267-A3292, 2018.

src/pymor/algorithms/randrangefinder.py

@@ -0,0 +1,143 @@
+# This file is part of the pyMOR project (http://www.pymor.org).
+# Copyright 2013-2019 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.interfaces import OperatorInterface
+
+
+@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 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, matmat=mv)
+        Linv = LinearOperator((source_product.range.dim, source_product.source.dim), matvec=mvinv, matmat=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
+
+
+@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 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):
+        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

src/pymortests/randrangefinder.py

@@ -0,0 +1,47 @@
+# This file is part of the pyMOR project (http://www.pymor.org).
+# Copyright 2013-2019 pyMOR developers and contributors. All rights reserved.
+# License: BSD 2-Clause License (http://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.data)
+    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.data)
+    assert B in op.range