Low-rank RADI solver for Riccati equations

[lbalicki]

This is an implementation of the low-rank RADI algorithm mentioned in #692. One issue I encountered was that the other solvers can be applied to equations of the type AXE^T+EXA^T-(EXC^T+S) R^{-1}(EXC^T+S)^T+BB^T=0, whereas the RADI implementation can not deal with the matrices S and R. Since LQGBTReductor does not need S and R as well for solving Riccati equations, this should not be too big of an issue. But I'm not sure if my workaround in riccati.py is appropriate.

[renefritze]

Please rebase on current master to remove the failing MPI jobs in CI

[pmli]

@lbalicki Notice that the more general Riccati equation with S and R (where R is positive definite) can be rewritten as a simpler Riccati equation with new A, B B^T, and C^T C (A - S R^-1 C, B B^T - S R^-1 S^T, and C^T R^-1 C if I'm not mistaken). The difficulty with nonzero S is that the new A has a low-rank update. I'm working on a "low-rank updated operator" which would use the Sherman-Morrison-Woodbury formula in apply_inverse. For now, I think you can remove the elif config.HAVE_PYMESS: part in riccati.py (note that pymess also doesn't support general S and R, see bindings.pymess.solve_ricc_lrcf).

[pmli]

It seems lrradi.py has Windows style new lines (I see ^M at the end of almost every line in the editor I use). Could you replace them with Unix style?
Below are some other general comments. I'll take a look at the implementation in more detail later.

[lbalicki]

I adapted the implementation to the requested changes. For now RADI is dealing with matrices S and R exactly like bindings.pymess.solve_ricc_lrcf.

[pmli]

The implementation looks very good. It seems you found and fixed a typo in how the paper used SMW. I will ask Jens next week for a review.

I just have very minor comments about the code below:

• in a few docstrings, there is an empty line missing before Returns
• np.eye is used much more often in pyMOR compared to np.identity
• there are a few places with expressions like 1 / a * b, which could be replaced with b / a to reduce the number of operations and increase accuracy a bit
• the expression RF += V * a can be converted to RF.axpy(a, V), which should avoid forming a new VectorArray for V * a
• np.block can be replaced with np.vstack in a few places to reduce the number of brackets

[pmli]

I talked with @drittelhacker a bit and he didn't point to anything that should be changed.

[codecov]

Codecov Report

Merging #736 into master will increase coverage by 0.35%.
The diff coverage is 97.6%.

Impacted Files	Coverage Δ
src/pymor/algorithms/riccati.py	90.47% <100%> (+0.31%)	⬆️
src/pymor/algorithms/lrradi.py	97.56% <97.56%> (ø)
src/pymor/vectorarrays/numpy.py	83.79% <0%> (+0.75%)	⬆️
src/pymor/vectorarrays/constructions.py	100% <0%> (+85.71%)	⬆️

[pmli]

@lbalicki Please rebase and squash the commits changing lrradi.py to simplify the diffs related to line endings. Probably the easiest way is

git checkout radi
git rebase master
git reset master

and then git add+git commit for each of the changed files.

[lbalicki]

@pmli I now put all my changes into one commit. Or would it be better to have an individual commit for each changed file?

[pmli]

@lbalicki For this situation, I would say lrradi.py and bibliography.rst can be in the same commit, followed by algorithms/riccati.py and pymortests/riccati.py in separate commits.

[renefritze]

CI is blocked by jupyterhub/repo2docker#755
Hotfix incoming: #744

[pmli]

@lbalicki Another rebase on master should make the tests pass.

[pmli]

Great, thanks!