add docstring and copy argument

docs/source/bibliography.bib

@@ -204,6 +204,19 @@ @Article{DPW13
  year = {2013}
 }
 
+@article{FKNYY20,
+  title = {Shifted {{Cholesky QR}} for {{Computing}} the {{QR Factorization}} of {{Ill-Conditioned Matrices}}},
+  author = {Fukaya, Takeshi and Kannan, Ramaseshan and Nakatsukasa, Yuji and Yamamoto, Yusaku and Yanagisawa, Yuka},
+  year = {2020},
+  journal = {SIAM Journal on Scientific Computing},
+  volume = {42},
+  number = {1},
+  pages = {A477-A503},
+  issn = {1064-8275, 1095-7197},
+  doi = {10.1137/18M1218212},
+  url = {https://epubs.siam.org/doi/10.1137/18M1218212},
+}
+
 @article{GP05,
   title={A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations},
   author={Grepl, Martin A and Patera, Anthony T},

docs/source/substitutions.py

@@ -145,6 +145,7 @@
 .. |Array| replace:: :class:`NumPy array <numpy.ndarray>`
 
 .. |SciPy| replace:: :mod:`SciPy <scipy>`
+.. |SciPy linalg| replace:: :mod:`SciPy linalg <scipy.linalg>`
 .. |SciPy spmatrix| replace:: :class:`SciPy spmatrix <scipy.sparse.spmatrix>`
 .. |SciPy spmatrices| replace:: :class:`SciPy spmatrices <scipy.sparse.spmatrix>`
 .. |Scipy spmatrix| replace:: :class:`SciPy spmatrix <scipy.sparse.spmatrix>`

src/pymor/algorithms/cholesky_qr.py

@@ -7,7 +7,7 @@
 from pymor.core.cache import CacheableObject, cached
 from pymor.core.logger import getLogger
 
-logger = getLogger('pymor.algorithms.cholesky_qr.cholesky_qr')
+logger = getLogger('pymor.algorithms.cholesky_qr.shifted_cholqr')
 
 
 class _ProductNorm(CacheableObject):
@@ -56,30 +56,89 @@ def _shifted_cholesky(X, shifter, check_finite=True):
     return R
 
 
-def cholesky_qr(A, product=None, maxiter=5, tol=None, check_finite=True, return_R=True, offset=0):
-    """No copy kwarg, because we cannot work in place anyway."""
+def shifted_cholqr(A, product=None, return_R=True, maxiter=3, offset=0, check_tol=None, check_finite=True, copy=True):
+    r"""Orthonormalize a |VectorArray| using the shifted CholeskyQR algorithm.
+
+    This method computes a QR decomposition of a |VectorArray| via Cholesky factorizations
+    of its Gramian matrix according to :cite:`FKNYY20`. For ill-conditioned matrices, the Cholesky
+    factorization will break down. In this case a diagonal shift will be applied to the Gramian.
+
+    - `shifted_cholqr(A, maxiter=3, check_tol=None)` is equivalent to the shifted CholeskyQR3
+    algorithm (Algorithm 4.2 in :cite:`FKNYY20`).
+    - `shifted_cholqr(A, maxiter=np.inf, check_tol=<some_number>)` is equivalent to the shifted
+    CholeskyQR algorithm (Algorithm 4.1 in :cite:`FKNYY20`).
+    - `shifted_cholqr(A, product=<some_product>, maxiter=3, check_tol=None)` is equivalent to the
+    shifted CholeskyQR3 algorithm in an oblique inner product (Algorithm 5.1 in :cite:`FKNYY20`).
+
+    .. note::
+        Note that the regular single iteration CholeskyQR algorithm is unstable in generel.
+        Setting `maxiter=1` will perform the simple shifted CholeskyQR algorithm (Algorithm 2.1
+        in :cite:`FKNYY20`) which is stable but leads to poor orthogonality in general, i.e.
+
+        .. math::
+            \lVert Q^TQ - I \rVert_2 < 2,
+
+        (see Lemma 3.1 in :cite:`FKNYY20`).
+
+    Parameters
+    ----------
+    A
+        The |VectorArray| which is to be orthonormalized.
+    product
+        The inner product |Operator| w.r.t. which to orthonormalize.
+        If `None`, the Euclidean product is used.
+    return_R
+        If `True`, the R matrix from QR decomposition is returned.
+    maxiter
+        Maximum number of iterations. Defaults to 3.
+    offset
+        Assume that the first `offset` vectors are already orthonormal and apply the
+        algorithm for the vectors starting at `offset + 1`.
+    check_tol
+        If not `None`, check if the resulting |VectorArray| is really orthornormal and
+        repeat the algorithm until the check passes
+        or `maxiter` is reached.
+    check_finite
+        This argument is passed down to |SciPy linalg| functions. Disabling may give a
+        performance gain, but may result in problems (crashes, non-termination) if the
+        inputs do contain infinities or NaNs.
+    copy
+        If `True`, create a copy of `A` instead of modifying `A` in-place.
+
+    Returns
+    -------
+    Q
+        The orthonormalized |VectorArray|.
+    R
+        The upper-triangular/trapezoidal matrix (if `compute_R` is `True`).
+    """
     assert 0 <= offset < len(A)
+    assert 0 < maxiter
+    assert check_tol is None or 0 < check_tol
+
+    if maxiter == 1:
+        logger.warning('Single iteration shifted CholeskyQR can lead to poor orthogonality!')
+
+    if copy:
+        A = A.copy()
 
     iter = 1
     product_norm = _ProductNorm(product)
     while iter <= maxiter:
         with logger.block(f'Iteration {iter}'):
-            if tol is None or iter == 1:
+            if check_tol is None or iter == 1:
                 # compute only the necessary parts of the Gramian
-                A_orth, A_todo = A[:offset].copy(), A[offset:]
-                B, X = np.split(A_todo.inner(A, product=product), [offset], axis=1)
+                B, X = np.split(A[offset:].inner(A, product=product), [offset], axis=1)
 
             # compute Cholesky factor of lower right block
-            shifter = _CholeskyShifter(A_todo, product, product_norm)
+            shifter = _CholeskyShifter(A[offset:], product, product_norm)
             Rx = _shifted_cholesky(X, shifter, check_finite=check_finite)
 
             # orthogonalize
             Rinv = spla.lapack.dtrtri(Rx)[0].T
-            if offset == 0:
-                A = A.lincomb(Rinv)
-            else:
-                A_orth.append(A_orth.lincomb(-Rinv@B) + A_todo.lincomb(Rinv))
-                A = A_orth
+            A_todo = A[:offset].lincomb(-Rinv@B) + A[offset:].lincomb(Rinv)
+            del A[offset:]
+            A.append(A_todo)
 
             # update blocks of R
             if return_R:
@@ -91,12 +150,11 @@ def cholesky_qr(A, product=None, maxiter=5, tol=None, check_finite=True, return_
                     spla.blas.dtrmm(1, Rx, Ri, overwrite_b=True)
 
             # check orthogonality
-            if tol is not None:
-                A_orth, A_todo = A[:offset].copy(), A[offset:]
-                B, X = np.split(A_todo.inner(A, product=product), [offset], axis=1)
-                res = spla.norm(X - np.eye(len(A_todo)), ord='fro', check_finite=check_finite)
+            if check_tol is not None:
+                B, X = np.split(A[offset:].inner(A, product=product), [offset], axis=1)
+                res = spla.norm(X - np.eye(len(A) - offset), ord='fro', check_finite=check_finite)
                 logger.info(f'Residual = {res}')
-                if res <= tol*np.sqrt(len(A)):
+                if res <= check_tol*np.sqrt(len(A)):
                     break
 
             iter += 1