Repository: pymor/pymor

Issue:
Add "operator plus low-rank" operator class
This could be used in Riccati equation solvers when S is nonzero. There, an operator of the form A-SR<sup>-1</sup>C appears.

I would propose to add a class which would implement the operator A+BD<sup>-1</sup>C<sup>T</sup>, where A is an operator, B and C are vector arrays, and D is a NumPy array. This class would use the [Sherman–Morrison–Woodbury formula](https://en.wikipedia.org/wiki/Woodbury_matrix_identity) in `apply_inverse[_adjoint]`:
(A + B D<sup>-1</sup> C<sup>T</sup>)<sup>-1</sup> = A<sup>-1</sup> - A<sup>-1</sup> B (D + C<sup>T</sup> A<sup>-1</sup> B)<sup>-1</sup> C<sup>T</sup> A<sup>-1</sup>

I'm not sure what would be a good name for such a class. Maybe `OperatorPlusLowRank` or `LowRankUpdatedOperator`?

Discussion:
sdrave: I think I would go for adding rules to `algrithms.lincomb` at the moment. I dislike adding such a special case to `LincombOperator.apply_inverse`.

In the long run, this might be an argument for adding an `algorithms.inverse` `RuleTable`, but I am not so sure if there are enough applications. Another might be solving `s*E - A` but I would still prefer `apply_shifted_inverse` as this is an important operation in pyMOR.

Is there a better name than `middle`? Has this some name in some community?
drittelhacker: Just to add some additional complication: As an alternative to the SMW for the low-rank updated operator, there is also the augmented matrix approach (reformulation, such that the lr-updated operator is the schur complement). I recently learnt that at Virginia Tech they have successfully used that to avoid the numerical stability issues of SMW in some cases. 

Also, the `s*E-A` should maybe be a `p*E+q*A` to cover both the case for the TFM evaluation and the ADI/RKSM style linear system solves.

In principle I am in support of splitting the operators and algorithm driven inverses, in case we can not map the later to actual operators in the discretization backend.
pmli: > I think I would go for adding rules to `algrithms.lincomb` at the moment. I dislike adding such a special case to `LincombOperator.apply_inverse`.

Is the input to `assemble_lincomb` guaranteed to have nonzero coefficients? This would simplify the implementation.

> Is there a better name than `middle`? Has this some name in some community?

Not sure. In HOSVD, there is a "core" tensor. There is also CUR decomposition, but I don't know what "U" stands for. I'm open to other suggestions.
sdrave: > Is the input to assemble_lincomb guaranteed to have nonzero coefficients? This would simplify the implementation.

Not ATM, and I am not completely convinced how useful it is in general. It would be quite easy to add this as the first rule to `AssembleLincombRules`. Where would you need it?

> Not sure. In HOSVD, there is a "core" tensor. There is also CUR decomposition, but I don't know what "U" stands for. I'm open to other suggestions.

The inverse at the `M` seems a little strange. Where does this come from? Would it be an option to pass an `InverseOperator` where needed? - Without the inverse, I would probably vote for 'core'. Maybe even with the inverse.
pmli: > Not ATM, and I am not completely convinced how useful it is in general. It would be quite easy to add this as the first rule to `AssembleLincombRules`. Where would you need it?

This would avoid checking in individual rules for zero coefficients, e.g., if `LowRankOperators` appear only with zero coefficients, a new `LowRankUpdatedOperator` does not have to be created.

> The inverse at the `M` seems a little strange. Where does this come from?

It appears in some variants of balanced truncation and the inverse does not have to be computed in the SWM formula. Also, there is no pyMOR interface for `L.lincomb(np.linalg.inv(M))` which would not invert `M` directly, so it cannot be "hidden" in one of the low-rank factors.

> Would it be an option to pass an `InverseOperator` where needed?

The implementation of `LowRankOperator.apply` is simpler if `M` is a NumPy array...

sdrave: >>    Not ATM, and I am not completely convinced how useful it is in general. It would be quite easy to add this as the first rule to AssembleLincombRules. Where would you need it?

> This would avoid checking in individual rules for zero coefficients, e.g., if LowRankOperators appear only with zero coefficients, a new LowRankUpdatedOperator does not have to be created.

Ok, feel free to add a rule if you feel the need.

> >    Would it be an option to pass an InverseOperator where needed?

Oh, sorry, I wasn't aware this is a NumPy array (which of course makes sense). Maybe it would be an option to have the inverse as an option. I assume there might be other applications for a `LowRankOperator` (think of SVD).
pmli: > Maybe it would be an option to have the inverse as an option. I assume there might be other applications for a `LowRankOperator` (think of SVD).

I was thinking about this. In that case, a linear combination of `LowRankOperators` becomes complicated and it's not clear what to do if one operator has an inverse and one doesn't.

Patch:
diff --git a/src/pymor/algorithms/lincomb.py b/src/pymor/algorithms/lincomb.py
index 6cdb9d2b4..cb086cdeb 100644
--- a/src/pymor/algorithms/lincomb.py
+++ b/src/pymor/algorithms/lincomb.py
@@ -2,13 +2,18 @@
 # Copyright 2013-2019 pyMOR developers and contributors. All rights reserved.
 # License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
 
+from itertools import chain
+
 import numpy as np
+import scipy.linalg as spla
 
 from pymor.algorithms.rules import RuleTable, match_generic, match_class_all, match_class_any, match_always
 from pymor.core.exceptions import RuleNotMatchingError
 from pymor.operators.block import (BlockOperator, BlockRowOperator, BlockColumnOperator, BlockOperatorBase,
                                    BlockDiagonalOperator, SecondOrderModelOperator, ShiftedSecondOrderModelOperator)
-from pymor.operators.constructions import ZeroOperator, IdentityOperator, VectorArrayOperator, LincombOperator
+from pymor.operators.constructions import (ZeroOperator, IdentityOperator, VectorArrayOperator, LincombOperator,
+                                           LowRankOperator, LowRankUpdatedOperator)
+from pymor.vectorarrays.constructions import cat_arrays
 
 
 def assemble_lincomb(operators, coefficients, solver_options=None, name=None):
@@ -55,6 +60,20 @@ class AssembleLincombRules(RuleTable):
         super().__init__(use_caching=False)
         self.__auto_init(locals())
 
+    @match_always
+    def action_zero_coeff(self, ops):
+        if all(coeff != 0 for coeff in self.coefficients):
+            raise RuleNotMatchingError
+        without_zero = [(op, coeff)
+                        for op, coeff in zip(ops, self.coefficients)
+                        if coeff != 0]
+        if len(without_zero) == 0:
+            return ZeroOperator(ops[0].range, ops[0].source, name=self.name)
+        else:
+            new_ops, new_coeffs = zip(*without_zero)
+            return assemble_lincomb(new_ops, new_coeffs,
+                                    solver_options=self.solver_options, name=self.name)
+
     @match_class_any(ZeroOperator)
     def action_ZeroOperator(self, ops):
         without_zero = [(op, coeff)
@@ -167,6 +186,64 @@ class AssembleLincombRules(RuleTable):
                 blocks[i, j] = ops[0].blocks[i, j] * c
             return operator_type(blocks)
 
+    @match_generic(lambda ops: sum(1 for op in ops if isinstance(op, LowRankOperator)) >= 2)
+    def action_merge_low_rank_operators(self, ops):
+        low_rank = []
+        not_low_rank = []
+        for op, coeff in zip(ops, self.coefficients):
+            if isinstance(op, LowRankOperator):
+                low_rank.append((op, coeff))
+            else:
+                not_low_rank.append((op, coeff))
+        inverted = [op.inverted for op, _ in low_rank]
+        if len(inverted) >= 2 and any(inverted) and any(not _ for _ in inverted):
+            return None
+        inverted = inverted[0]
+        left = cat_arrays([op.left for op, _ in low_rank])
+        right = cat_arrays([op.right for op, _ in low_rank])
+        core = []
+        for op, coeff in low_rank:
+            core.append(op.core)
+            if inverted:
+                core[-1] /= coeff
+            else:
+                core[-1] *= coeff
+        core = spla.block_diag(*core)
+        new_low_rank_op = LowRankOperator(left, core, right, inverted=inverted)
+        if len(not_low_rank) == 0:
+            return new_low_rank_op
+        else:
+            new_ops, new_coeffs = zip(*not_low_rank)
+            return assemble_lincomb(chain(new_ops, [new_low_rank_op]), chain(new_coeffs, [1]),
+                                    solver_options=self.solver_options, name=self.name)
+
+    @match_generic(lambda ops: len(ops) >= 2)
+    @match_class_any(LowRankOperator, LowRankUpdatedOperator)
+    def action_merge_into_low_rank_updated_operator(self, ops):
+        new_ops = []
+        new_lr_ops = []
+        new_coeffs = []
+        new_lr_coeffs = []
+        for op, coeff in zip(ops, self.coefficients):
+            if isinstance(op, LowRankOperator):
+                new_lr_ops.append(op)
+                new_lr_coeffs.append(coeff)
+            elif isinstance(op, LowRankUpdatedOperator):
+                new_ops.append(op.operators[0])
+                new_coeffs.append(coeff * op.coefficients[0])
+                new_lr_ops.append(op.operators[1])
+                new_lr_coeffs.append(coeff * op.coefficients[1])
+            else:
+                new_ops.append(op)
+                new_coeffs.append(coeff)
+        lru_op = assemble_lincomb(new_ops, new_coeffs)
+        lru_lr_op = assemble_lincomb(new_lr_ops, new_lr_coeffs)
+        lru_lr_coeff = 1
+        if isinstance(lru_lr_op, LincombOperator):
+            lru_lr_op, lru_lr_coeff = lru_lr_op.operators[0], lru_lr_op.coefficients[0]
+        return LowRankUpdatedOperator(lru_op, lru_lr_op, 1, lru_lr_coeff,
+                                      solver_options=self.solver_options, name=self.name)
+
     @match_always
     def action_call_assemble_lincomb_method(self, ops):
         id_coeffs, ops_without_id, coeffs_without_id = [], [], []
diff --git a/src/pymor/algorithms/to_matrix.py b/src/pymor/algorithms/to_matrix.py
index b5daf44a7..0b6f39ce0 100644
--- a/src/pymor/algorithms/to_matrix.py
+++ b/src/pymor/algorithms/to_matrix.py
@@ -10,7 +10,8 @@ import scipy.sparse.linalg as spsla
 from pymor.algorithms.rules import RuleTable, match_class
 from pymor.operators.block import BlockOperatorBase
 from pymor.operators.constructions import (AdjointOperator, ComponentProjection, Concatenation, IdentityOperator,
-                                           LincombOperator, VectorArrayOperator, ZeroOperator)
+                                           LincombOperator, LowRankOperator, LowRankUpdatedOperator,
+                                           VectorArrayOperator, ZeroOperator)
 from pymor.operators.numpy import NumpyMatrixOperator
 
 
@@ -139,6 +140,21 @@ class ToMatrixRules(RuleTable):
             res = res + op_coefficients[i] * self.apply(op.operators[i])
         return res
 
+    @match_class(LowRankOperator)
+    def action_LowRankOperator(self, op):
+        format = self.format
+        if not op.inverted:
+            res = op.left.to_numpy().T @ op.core @ op.right.to_numpy()
+        else:
+            res = op.left.to_numpy().T @ spla.solve(op.core, op.right.to_numpy())
+        if format is not None and format != 'dense':
+            res = getattr(sps, format + '_matrix')(res)
+        return res
+
+    @match_class(LowRankUpdatedOperator)
+    def action_LowRankUpdatedOperator(self, op):
+        return op.coeff * self.apply(op.operator) + op.lr_coeff * self.apply(op.lr_operator)
+
     @match_class(VectorArrayOperator)
     def action_VectorArrayOperator(self, op):
         format = self.format
diff --git a/src/pymor/operators/constructions.py b/src/pymor/operators/constructions.py
index ec82469c6..fa486a4ab 100644
--- a/src/pymor/operators/constructions.py
+++ b/src/pymor/operators/constructions.py
@@ -9,6 +9,7 @@ from itertools import chain
 from numbers import Number
 
 import numpy as np
+import scipy.linalg as spla
 
 from pymor.core.defaults import defaults
 from pymor.core.exceptions import InversionError
@@ -161,7 +162,7 @@ class LincombOperator(OperatorBase):
                 raise NotImplementedError
         derivative_coefficients = []
         for coef in self.coefficients:
-            if isinstance(coef,Parametric):
+            if isinstance(coef, Parametric):
                 derivative_coefficients.append(coef.d_mu(component, index))
             else:
                 derivative_coefficients.append(0.)
@@ -352,6 +353,156 @@ class Concatenation(OperatorBase):
         return Concatenation(operators, solver_options=other.solver_options)
 
 
+class LowRankOperator(OperatorBase):
+    """Non-parametric low-rank operator.
+
+    Represents an operator of the form :math:`L C R^H` or
+    :math:`L C^{-1} R^H` where :math:`L` and :math:`R` are
+    |VectorArrays| of column vectors and :math:`C` a 2D |NumPy array|.
+
+    Parameters
+    ----------
+    left
+        |VectorArray| representing :math:`L`.
+    core
+        |NumPy array| representing :math:`C`.
+    right
+        |VectorArray| representing :math:`R`.
+    inverted
+        Whether :math:`C` is inverted.
+    solver_options
+        The |solver_options| for the operator.
+    name
+        Name of the operator.
+    """
+
+    linear = True
+
+    def __init__(self, left, core, right, inverted=False, solver_options=None, name=None):
+        assert isinstance(left, VectorArrayInterface)
+        assert isinstance(right, VectorArrayInterface)
+        assert len(left) == len(right)
+        assert (isinstance(core, np.ndarray)
+                and core.ndim == 2
+                and core.shape[0] == core.shape[1] == len(left))
+
+        self.__auto_init(locals())
+        self.source = right.space
+        self.range = left.space
+
+    @property
+    def H(self):
+        options = {
+            'inverse': self.solver_options.get('inverse_adjoint'),
+            'inverse_adjoint': self.solver_options.get('inverse'),
+        } if self.solver_options else None
+        return type(self)(self.right,
+                          self.core.T.conj(),
+                          self.left,
+                          inverted=self.inverted,
+                          solver_options=options,
+                          name=self.name + '_adjoint')
+
+    def apply(self, U, mu=None):
+        assert U in self.source
+        V = self.right.dot(U)
+        if self.inverted:
+            V = spla.solve(self.core, V)
+        else:
+            V = self.core @ V
+        return self.left.lincomb(V.T)
+
+    def apply_adjoint(self, V, mu=None):
+        assert V in self.range
+        U = self.left.dot(V)
+        if self.inverted:
+            U = spla.solve(self.core.T.conj(), U)
+        else:
+            U = self.core.T.conj() @ U
+        return self.right.lincomb(U.T)
+
+
+class LowRankUpdatedOperator(LincombOperator):
+    r"""|Operator| plus :class:`LowRankOperator`.
+
+    Represents a linear combination of an |Operator| and
+    :class:`LowRankOperator`. Uses the Sherman-Morrison-Woodbury formula
+    in `apply_inverse` and `apply_inverse_adjoint`:
+
+    .. math::
+        \left(\alpha A + \beta L C R^H \right)^{-1}
+        & = \alpha^{-1} A^{-1}
+            - \alpha^{-1} \beta A^{-1} L C
+            \left(\alpha C + \beta C R^H A^{-1} L C \right)^{-1}
+            C R^H A^{-1}, \\
+        \left(\alpha A + \beta L C^{-1} R^H \right)^{-1}
+        & = \alpha^{-1} A^{-1}
+            - \alpha^{-1} \beta A^{-1} L
+            \left(\alpha C + \beta R^H A^{-1} L \right)^{-1}
+            R^H A^{-1}.
+
+    Parameters
+    ----------
+    operator
+        |Operator|.
+    lr_operator
+        :class:`LowRankOperator`.
+    coeff
+        A linear coefficient for `operator`. Can either be a fixed
+        number or a |ParameterFunctional|.
+    lr_coeff
+        A linear coefficient for `lr_operator`. Can either be a fixed
+        number or a |ParameterFunctional|.
+    solver_options
+        The |solver_options| for the operator.
+    name
+        Name of the operator.
+    """
+
+    def __init__(self, operator, lr_operator, coeff, lr_coeff,
+                 solver_options=None, name=None):
+        assert isinstance(lr_operator, LowRankOperator)
+        super().__init__([operator, lr_operator], [coeff, lr_coeff],
+                         solver_options=solver_options, name=name)
+        self.__auto_init(locals())
+
+    def apply_inverse(self, V, mu=None, least_squares=False):
+        if least_squares:
+            return super().apply_inverse(V, mu=mu, least_squares=True)
+        A, LR = self.operators
+        L, C, R = LR.left, LR.core, LR.right
+        if not LR.inverted:
+            L = L.lincomb(C.T)
+            R = R.lincomb(C.conj())
+        alpha, beta = self.evaluate_coefficients(mu)
+        AinvV = A.apply_inverse(V)
+        AinvL = A.apply_inverse(L)
+        mat = alpha * C + beta * R.dot(AinvL)
+        RhAinvV = R.dot(AinvV)
+        U = AinvV
+        U.axpy(-beta, AinvL.lincomb(spla.solve(mat, RhAinvV).T))
+        U.scal(1 / alpha)
+        return U
+
+    def apply_inverse_adjoint(self, U, mu=None, least_squares=False):
+        if least_squares:
+            return super().apply_inverse_adjoint(U, mu=mu, least_squares=True)
+        A, LR = self.operators
+        L, C, R = LR.left, LR.core, LR.right
+        if not LR.inverted:
+            L = L.lincomb(C.T)
+            R = R.lincomb(C.conj())
+        alpha, beta = (c.conjugate() for c in self.evaluate_coefficients(mu))
+        AinvhU = A.apply_inverse_adjoint(U)
+        AinvhR = A.apply_inverse_adjoint(R)
+        mat = alpha * C.T.conj() + beta * L.dot(AinvhR)
+        LhAinvhU = L.dot(AinvhU)
+        V = AinvhU
+        V.axpy(-beta, AinvhR.lincomb(spla.solve(mat, LhAinvhU).T))
+        V.scal(1 / alpha)
+        return V
+
+
 class ComponentProjection(OperatorBase):
     """|Operator| representing the projection of a |VectorArray| onto some of its components.
 


Test patch:
diff --git a/src/pymortests/low_rank_op.py b/src/pymortests/low_rank_op.py
new file mode 100644
index 000000000..27f074499
--- /dev/null
+++ b/src/pymortests/low_rank_op.py
@@ -0,0 +1,135 @@
+# 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
+import scipy.linalg as spla
+
+from pymor.algorithms.lincomb import assemble_lincomb
+from pymor.operators.constructions import LowRankOperator, LowRankUpdatedOperator
+from pymor.operators.numpy import NumpyMatrixOperator
+from pymor.vectorarrays.numpy import NumpyVectorSpace
+
+
+def construct_operators_and_vectorarrays(m, n, r, k, seed=0):
+    space_m = NumpyVectorSpace(m)
+    space_n = NumpyVectorSpace(n)
+    rng = np.random.RandomState(seed)
+    A = NumpyMatrixOperator(rng.randn(m, n))
+    L = space_m.random(r, distribution='normal', random_state=rng)
+    C = rng.randn(r, r)
+    R = space_n.random(r, distribution='normal', random_state=rng)
+    U = space_n.random(k, distribution='normal', random_state=rng)
+    V = space_m.random(k, distribution='normal', random_state=rng)
+    return A, L, C, R, U, V
+
+
+def test_low_rank_apply():
+    _, L, C, R, U, _ = construct_operators_and_vectorarrays(6, 5, 2, 3)
+
+    LR = LowRankOperator(L, C, R)
+    V = LR.apply(U)
+    assert np.allclose(V.to_numpy().T, L.to_numpy().T @ C @ (R.to_numpy() @ U.to_numpy().T))
+
+    LR = LowRankOperator(L, C, R, inverted=True)
+    V = LR.apply(U)
+    assert np.allclose(V.to_numpy().T,
+                       L.to_numpy().T @ spla.solve(C, R.to_numpy() @ U.to_numpy().T))
+
+
+def test_low_rank_apply_adjoint():
+    _, L, C, R, _, V = construct_operators_and_vectorarrays(6, 5, 2, 3)
+
+    LR = LowRankOperator(L, C, R)
+    U = LR.apply_adjoint(V)
+    assert np.allclose(U.to_numpy().T, R.to_numpy().T @ C.T @ (L.to_numpy() @ V.to_numpy().T))
+
+    LR = LowRankOperator(L, C, R, inverted=True)
+    U = LR.apply_adjoint(V)
+    assert np.allclose(U.to_numpy().T,
+                       R.to_numpy().T @ spla.solve(C.T, L.to_numpy() @ V.to_numpy().T))
+
+
+def test_low_rank_updated_apply_inverse():
+    A, L, C, R, _, V = construct_operators_and_vectorarrays(5, 5, 2, 3)
+
+    LR = LowRankOperator(L, C, R)
+    op = LowRankUpdatedOperator(A, LR, 1, 1)
+    U = op.apply_inverse(V)
+    mat = A.matrix + L.to_numpy().T @ C @ R.to_numpy()
+    assert np.allclose(U.to_numpy().T, spla.solve(mat, V.to_numpy().T))
+
+    LR = LowRankOperator(L, C, R, inverted=True)
+    op = LowRankUpdatedOperator(A, LR, 1, 1)
+    U = op.apply_inverse(V)
+    mat = A.matrix + L.to_numpy().T @ spla.solve(C, R.to_numpy())
+    assert np.allclose(U.to_numpy().T, spla.solve(mat, V.to_numpy().T))
+
+
+def test_low_rank_updated_apply_inverse_adjoint():
+    A, L, C, R, U, _ = construct_operators_and_vectorarrays(5, 5, 2, 3)
+
+    LR = LowRankOperator(L, C, R)
+    op = LowRankUpdatedOperator(A, LR, 1, 1)
+    V = op.apply_inverse_adjoint(U)
+    mat = A.matrix + L.to_numpy().T @ C @ R.to_numpy()
+    assert np.allclose(V.to_numpy().T, spla.solve(mat.T, U.to_numpy().T))
+
+    LR = LowRankOperator(L, C, R, inverted=True)
+    op = LowRankUpdatedOperator(A, LR, 1, 1)
+    V = op.apply_inverse_adjoint(U)
+    mat = A.matrix + L.to_numpy().T @ spla.solve(C, R.to_numpy())
+    assert np.allclose(V.to_numpy().T, spla.solve(mat.T, U.to_numpy().T))
+
+
+def test_low_rank_assemble():
+    r1, r2 = 2, 3
+    _, L1, C1, R1, _, _ = construct_operators_and_vectorarrays(5, 5, r1, 0)
+    _, L2, C2, R2, _, _ = construct_operators_and_vectorarrays(5, 5, r2, 0, seed=1)
+
+    LR1 = LowRankOperator(L1, C1, R1)
+    LR2 = LowRankOperator(L2, C2, R2)
+    op = assemble_lincomb([LR1, LR2], [1, 1])
+    assert isinstance(op, LowRankOperator)
+    assert len(op.left) == r1 + r2
+    assert not op.inverted
+
+    op = (LR1 + (LR1 + LR2) + LR2).assemble()
+    assert isinstance(op, LowRankOperator)
+
+    LR1 = LowRankOperator(L1, C1, R1, inverted=True)
+    LR2 = LowRankOperator(L2, C2, R2, inverted=True)
+    op = assemble_lincomb([LR1, LR2], [1, 1])
+    assert isinstance(op, LowRankOperator)
+    assert len(op.left) == r1 + r2
+    assert op.inverted
+
+    LR1 = LowRankOperator(L1, C1, R1, inverted=True)
+    LR2 = LowRankOperator(L2, C2, R2)
+    op = assemble_lincomb([LR1, LR2], [1, 1])
+    assert op is None
+
+
+def test_low_rank_updated_assemble():
+    A, L, C, R, _, _ = construct_operators_and_vectorarrays(5, 5, 2, 0)
+    LR = LowRankOperator(L, C, R)
+
+    op = (A + LR).assemble()
+    assert isinstance(op, LowRankUpdatedOperator)
+
+    op = (A + LR + LR).assemble()
+    assert isinstance(op, LowRankUpdatedOperator)
+
+    op = (A + (A + LR).assemble() + LR).assemble()
+    assert isinstance(op, LowRankUpdatedOperator)
+
+
+def test_low_rank_updated_assemble_apply():
+    A, L, C, R, U, _ = construct_operators_and_vectorarrays(5, 5, 2, 3)
+
+    LR = LowRankOperator(L, C, R)
+    op = (A + (A + LR).assemble() + LR).assemble()
+    V = op.apply(U)
+    assert np.allclose(V.to_numpy().T,
+                       2 * A.matrix @ U.to_numpy().T
+                       + 2 * L.to_numpy().T @ C @ (R.to_numpy() @ U.to_numpy().T))
diff --git a/src/pymortests/to_matrix.py b/src/pymortests/to_matrix.py
index a85088c5c..f12d400bd 100644
--- a/src/pymortests/to_matrix.py
+++ b/src/pymortests/to_matrix.py
@@ -8,8 +8,8 @@ import scipy.sparse as sps
 
 from pymor.algorithms.to_matrix import to_matrix
 from pymor.operators.block import BlockOperator, BlockDiagonalOperator
-from pymor.operators.constructions import (AdjointOperator, ComponentProjection, IdentityOperator, VectorArrayOperator,
-                                           ZeroOperator)
+from pymor.operators.constructions import (AdjointOperator, ComponentProjection, IdentityOperator, LowRankOperator,
+                                           LowRankUpdatedOperator, VectorArrayOperator, ZeroOperator)
 from pymor.operators.numpy import NumpyMatrixOperator
 from pymor.vectorarrays.numpy import NumpyVectorSpace
 
@@ -200,6 +200,42 @@ def test_to_matrix_LincombOperator():
     assert_type_and_allclose(C, Cop, 'sparse')
 
 
+def test_to_matrix_LowRankOperator():
+    np.random.seed(0)
+    m = 6
+    n = 5
+    r = 2
+    L = np.random.randn(m, r)
+    Lva = NumpyVectorSpace.make_array(L.T)
+    C = np.random.randn(r, r)
+    R = np.random.randn(n, r)
+    Rva = NumpyVectorSpace.make_array(R.T)
+
+    LR = LowRankOperator(Lva, C, Rva)
+    assert_type_and_allclose(L @ C @ R.T, LR, 'dense')
+
+    LR = LowRankOperator(Lva, C, Rva, inverted=True)
+    assert_type_and_allclose(L @ spla.solve(C, R.T), LR, 'dense')
+
+
+def test_to_matrix_LowRankUpdatedOperator():
+    np.random.seed(0)
+    m = 6
+    n = 5
+    r = 2
+    A = np.random.randn(m, n)
+    Aop = NumpyMatrixOperator(A)
+    L = np.random.randn(m, r)
+    Lva = NumpyVectorSpace.make_array(L.T)
+    C = np.random.randn(r, r)
+    R = np.random.randn(n, r)
+    Rva = NumpyVectorSpace.make_array(R.T)
+    LR = LowRankOperator(Lva, C, Rva)
+
+    op = LowRankUpdatedOperator(Aop, LR, 1, 1)
+    assert_type_and_allclose(A + L @ C @ R.T, op, 'dense')
+
+
 def test_to_matrix_VectorArrayOperator():
     np.random.seed(0)
     V = np.random.randn(10, 2)