Low-rank updated operators and Sherman-Morrison-Woodbury formula

[pmli]

So far, this adds LowRankOperator and LowRankUpdatedOperator to operators.constructions. The LowRankUpdatedOperator uses Sherman-Morrison-Woodbury formula in apply_inverse and apply_inverse_adjoint.

The motivation for L M^-1 R^H form of the low-rank operator is from Riccati equations appearing in variants of balanced truncation (LQG balanced truncation with nonzero D, balanced stochastic truncation, bounded real balanced truncation with nonzero D, ...). In particular, I'm considering removing the S part in riccati.solve_ricc_lrcf and instead expect it to be given as a low-rank update in A.

What remains is making LincombOperator and/or algorithms.lincomb aware of these operators, which I'm not sure how to do. One approach might be adding at least three rules to lincomb:

• if all operators are LowRankOperators, return a single LowRankOperator
• if one of the operators is a LowRankOperator, return a LowRankUpdatedOperator
• if there is a LowRankUpdatedOperator and another operator, return a new LowRankUpdatedOperator

Another approach might be removing LowRankUpdatedOperator and making LincombOperator detect LowRankOperators and use SMW formula.

@sdrave Thoughts?

Closes #502.

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

[codecov]

Codecov Report

Merging #743 into master will increase coverage by 0.23%.
The diff coverage is 96.37%.

Impacted Files	Coverage Δ
src/pymor/algorithms/lincomb.py	74.52% <100%> (+16.35%)	⬆️
src/pymor/algorithms/to_matrix.py	99.09% <90%> (-0.91%)	⬇️
src/pymor/operators/constructions.py	86.44% <94.11%> (+1.04%)	⬆️
src/pymor/tools/random.py	100% <0%> (+6.66%)	⬆️

[pmli]

I removed the option for core to be the identity matrix by default to avoid issues of whether it is inverted or not in algorithms.lincomb. The disadvantage is that this makes apply methods less efficient because of additional multiplications with an identity matrix.

@sdrave Do you have an idea how to better test lincomb for low-rank operators?

[sdrave]

I think you should also apply some random vectors and see if op does the right thing.

[pmli]

@sdrave Thanks, I made the changes.

[sdrave]

The changes look good. However, I still think there should be better tests for assemble_lincomb. Basically, I would apply some random vectors to the assembled operator and compare to the result of manually applying the vector to the summands of the linear combination.

[pmli]

Aha, ok, I thought you meant only the operators themselves.
I've also extended algorithms.to_matrix.