[algorithms.samdp] add samdp algorithm

[lbalicki]

I implemented the the subspace accelerated dominant pole algorithm (https://ieeexplore.ieee.org/document/1717547) which computes the dominant poles of the transfer function of a LTI system. It can be used in the future to add modal MOR methods to pyMOR.

[pmli]

Hi @lbalicki, thanks for the PR. Adding a test in src/pymortests/algorithms/ would be good (something simple, maybe a random or a specially constructed example).

[pmli]

I gave mostly minor cosmetic comments. I'll take a closer look at the algorithm later.

[lbalicki]

I made the suggested changes and in order to be more consistent with the other algorithms I decided to change the B and C parameters to VectorArrays.

[pmli]

Here are some mostly minor issues (found using pylint).

[pmli]

Some more minor remarks below.

I tried running the code with the building model from SLICOT:

from pymor.basic import *
from pymor.algorithms.samdp import samdp
build = LTIModel.from_mat_file('build.mat')
dp = samdp(build.A, build.E, build.B.as_range_array(), build.C.as_source_array(), 2)

but I get NaNs appearing.

[lbalicki]

Some more minor remarks below.

I tried running the code with the building model from SLICOT:

from pymor.basic import *
from pymor.algorithms.samdp import samdp
build = LTIModel.from_mat_file('build.mat')
dp = samdp(build.A, build.E, build.B.as_range_array(), build.C.as_source_array(), 2)

but I get NaNs appearing.

I had set the default initial shift to zero, which seems to cause problems in some cases. Now, if no initial shifts are provided, a random complex shift is used. For me the model from SLICOT works this way.

[pmli]

Some improvements concerning NumPy and converting from Matlab to Python style (e.g., it is more natural for residues[i], i.e. residues[i, :, :], to be the ith residue instead of residues[:, :, i]).

[pmli]

Thank you for the changes, @lbalicki. I've added some more comments below.

I've tried the code with almost all SLICOT benchmarks and it works very well except for the CDplayer benchmark. For that example, samdp gave me a result only when nwanted is 1. @lbalicki, do you get the same? Do you see if changing some options could improve things? I guess it could simply be a badly conditioned problem (LR-ADI also has issues with it, as far as I remember). So, either way, I would be for merging after the changes below.

[lbalicki]

It seems like there was still an issue when the eigenvectors of the dominant poles were complex. With the changes I have made, the CDplayer model from SLICOT works even with large values for nwanted.

[pmli]

Ok, great, CDplayer is working nicely, and a few other benchmarks also perform better (e.g., random). I have a few more comments now:

1. I noticed that multiple runs can give different results, at least in terms of convergence. Therefore, it would be good to give the user the option to pick the seed for np.random.uniform, or maybe just give a suggestion in the documentation to try setting a different initial shift if there is no convergence.
2. There are two more methods here (twosided_rqi and select_max_eig) that might be better as "private" methods (i.e., renaming to _twosided_rqi and _select_max_eig) if users are not expected to use them directly.
3. I tried samdp on filter2D and rail1357 benchmarks from the Oberwolfach collection. For larger nwanted (e.g., 20), it returns poles with non-zero imaginary parts (although, relatively small), even though both models have only real poles. Do you see a way to remedy this?

[lbalicki]

Ok, great, CDplayer is working nicely, and a few other benchmarks also perform better (e.g., random). I have a few more comments now:

1. I noticed that multiple runs can give different results, at least in terms of convergence. Therefore, it would be good to give the user the option to pick the seed for `np.random.uniform`, or maybe just give a suggestion in the documentation to try setting a different initial shift if there is no convergence.

2. There are two more methods here (`twosided_rqi` and `select_max_eig`) that might be better as "private" methods (i.e., renaming to `_twosided_rqi` and `_select_max_eig`) if users are not expected to use them directly.

3. I tried `samdp` on [filter2D](https://morwiki.mpi-magdeburg.mpg.de/morwiki/index.php/Tunable_Optical_Filter) and [rail1357](https://morwiki.mpi-magdeburg.mpg.de/morwiki/index.php/Steel_Profile) benchmarks from the Oberwolfach collection. For larger `nwanted` (e.g., 20), it returns poles with non-zero imaginary parts (although, relatively small), even though both models have only real poles. Do you see a way to remedy this?

Regarding 1. and 2. I just added a seed=0 option and I changed the name of the functions.

For 3. I changed the default tolerances to tol=1e-11 and imagtol=1e-7 which makes the poles of both examples have zero imaginary part. In general the user might have to adjust these tolerances to get accurate results. I also tried to use relative error tolerances in some parts of the algorithm, but that did not lead to better results for these models either.

[pmli]

For 3. I changed the default tolerances to tol=1e-11 and imagtol=1e-7 which makes the poles of both examples have zero imaginary part. In general the user might have to adjust these tolerances to get accurate results. I also tried to use relative error tolerances in some parts of the algorithm, but that did not lead to better results for these models either.

Ok, interesting. I guess the algorithm is really sensitive to the choices of the tolerances. I tried the new tolerances for the 'random' example in the SLICOT benchmark and it looks like samdp can return less than nwanted poles. E.g., if nwanted is 10, for tol=1e-11 it returns 4, while for tol=1e-10 it returns 20. @lbalicki Do you get the same? Is this a feature or a bug?

Also, it would be good to return the test for which='LM'. I see that one test would fail. Maybe the assert condition can be relaxed or the tolerances adjusted.

[lbalicki]

It seems like the problem with the tol option is that it can happen that the _twosided_rqi function is not able to converge (i.e., the residual will not drop below tol no matter how many iterations are performed). In this case, restarts are performed until the maximum amount of restarts is reached and samdp can not return the desired amount of poles.

My suggestion would be to print a warning if _twosided_rqi does not converge to the desired tolerance, but continue with the computed poles anyway. This works well for the 'random' SLICOT benchmark. Does this make sense?

For the tests I would maybe change assert np.average(val) < np.average(dom_val) to assert np.median(val) < np.median(dom_val) instead of relaxing the condition.

[pmli]

My suggestion would be to print a warning if _twosided_rqi does not converge to the desired tolerance, but continue with the computed poles anyway. This works well for the 'random' SLICOT benchmark. Does this make sense?

It makes sense to me.

[pmli]

Now the rail1357 example (for nwanted equal to 20) looks wrong. It seems samdp does not return poles in conjugate pairs.

[lbalicki]

I made a few changes which make the algorithm more stable. Also, dealing with complex poles works better now.

[pmli]

Now the rail1357 example (for nwanted equal to 20) looks wrong. It seems samdp does not return poles in conjugate pairs.

It looks like there are now similar issues for the transmission line model from SLICOT. For nwanted=50, it returns only real poles, and with imagtol=1e-8 it doesn't return conjugate pairs.

[lbalicki]

Comparing with the poles computed by scipy.linalg.eig it looks like it is correct that the poles are real for nwanted=50 in that example. Nonetheless, I fixed the issue with the conjugate pairs and made some other improvements as well.

[pmli]

These are the poles I get for tline model, colored by dominance (log10 of norm(residual)/norm(pole.real)):

This is the result of current samdp for nwanted=50:

and using an older version (39e0a03):

The older version looks better to me...

The results were generated using tline_samdp.zip

[codecov]

Codecov Report

Merging #834 into master will not change coverage.
The diff coverage is n/a.

[lbalicki]

It seems that this happened because I switched to relative tolerances for the residuals. After some testing it seems that absolute tolerances work better for almost all examples.

One example that does not work better is the filter2D model. I get a complex pole for tol=1e-10 and nwanted=50 with the newest version. The model works fine for tol=1e-15 but then other models fail (e.g. the building model from SLICOT). Overall, it seems to be difficult to balance the tol option such that all examples yield good results and regarding this, the implementation here behaves similarly. Is this behaviour acceptable?

[pmli]

I can confirm that samdp works well for all the SLICOT benchmarks I tried.

From what I've heard from @drittelhacker, sensitivity to tolerances is somewhat acceptable with this algorithm.

[drittelhacker]

I can confirm that samdp works well for all the SLICOT benchmarks I tried.

From what I've heard from @drittelhacker, sensitivity to tolerances is somewhat acceptable with this algorithm.

I would not exactly say "acceptable" but "to be expected"

[sdrave]

Looks very good to me, although I cannot really comment on the algorithmic details. @lbalicki, you should add an entry to AUTHORS.md ..

[pmli]

It's merged. Thank you, @lbalicki.