Adaptive RandomizedERAReductor

[artpelling]

This PR adds a new RandomizedERAReductor based on (https://doi.org/10.1137/20M1327616).

I've also corrected a mistake in the error bound of the original ERAReductor that originated from a typo in the original paper.

The RandomizedERAReductor can also be used in an adaptive way by using a leave-one-out error estimator that I will introduce in an upcoming paper. The estimator is not a true upper bound for the error but very fast to evaluate and performs much better than all other available estimators in practice (the established error bounds are grossly conservative).

[artpelling]

The status of the branch is still pre "Great Transposition", I will have to carefully migrate everything after my holiday.

[artpelling]

@sdrave Thanks for transposing. I've checked out your changes and they look good. I could reproduce the results from my paper, so numerically (at least with the LOO estimator) it is fine. I will add descriptive docstrings and tests soon!

[codecov]

Codecov Report

❌ Patch coverage is 34.06593% with 60 lines in your changes missing coverage. Please review.
✅ Project coverage is 88.08%. Comparing base (648a321) to head (e1ffac8).

Files with missing lines	Patch %	Lines
src/pymor/reductors/era.py	28.91%	59 Missing ⚠️
src/pymor/algorithms/rand_la.py	87.50%	1 Missing ⚠️

Additional details and impacted files

Flag	Coverage Δ
github_actions	78.51% <34.06%> (-0.10%)	⬇️
gitlab_ci	88.00% <34.06%> (-0.12%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing lines	Coverage Δ
src/pymor/algorithms/rand_la.py	89.17% <87.50%> (+0.97%)	⬆️
src/pymor/reductors/era.py	42.35% <28.91%> (-13.79%)	⬇️

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[artpelling]

@sdrave I am pretty happy with it now. Maybe there is a more elegant way using pyMOR defaults to set the rrf_opts?

[sdrave]

Thanks, @artpelling. I'm very busy this week. Will try to look into the PR by the end of next week.

[sdrave]

Hi @artpelling, sorry for the long waiting time. Could you move the changes in rand_la.py to a different PR and briefly explain the changes made?

[sdrave]

With #2556 being merged, could you rebase on main, @artpelling?

[sdrave]

Hi @artpelling, thanks again for the PR. I have taken a closer look now. There are two things I'd like to discuss:

• You are doing the actual SVD part of the randomized SVD in reduce instead of calling randomized_svd. If I understand correctly, the reasons for this are that:

  1. randomized_svd does not support error tolerances,
  2. randomized_svd cannot be called repeatedly for different tolerances/ranks
     without completely recomputing the approximate range of _H.

  I feel that both points would be better addressed by improving randomized_svd. At least the missing tolerance parameter(s) should be definitely added there. Regarding different ranks: do you really need the ability to increase the rank a posteriori? (Decreasing it could be done by simply caching the SVD and truncating it as needed.) If yes, I think we should consider something like a RandomizedSVD class, or the possiblity to pass/return a RandomizedRangeFinder to/from randomized_svd.

• Can you explain why you need num_left/num_right in the constructor for the randomized version? If this is only about recycling the range finder, you could also keep theses parameters in reduce and throw the range finder away in case their values change. Alternatively, one could also move these parameters to the constructor for the deterministic version?

I think that if we address both points, then both ERA versions would share a lot more code, making it easier to see where the differences lie. Also maintaining the code would be easier ... - What do you think?

[artpelling]

• You are doing the actual SVD part of the randomized SVD in reduce instead of calling randomized_svd. If I understand correctly, the reasons for this are that:

  1. randomized_svd does not support error tolerances,
  2. randomized_svd cannot be called repeatedly for different tolerances/ranks
     without completely recomputing the approximate range of _H.

Yes.

I feel that both points would be better addressed by improving randomized_svd. At least the missing tolerance parameter(s) should be definitely added there.

The design of randomized_svd, as many other methods, does not fit the reductor style with caching and reusing computations. I already refactored the rand_la module a lot to accommodate RandomizedERA, so I drew the line there.

Regarding different ranks: do you really need the ability to increase the rank a posteriori? (Decreasing it could be done by simply caching the SVD and truncating it as needed.)

Yes that is the key feature for me. Knowing the target model order is impossible for me. The appeal of the adaptive pipeline is that I can always come back later to refine.

Can you explain why you need num_left/num_right in the constructor for the randomized version? If this is only about recycling the range finder, you could also keep theses parameters in reduce and throw the range finder away in case their values change. Alternatively, one could also move these parameters to the constructor for the deterministic version?

I want it in the constructor for the reason you described. I dislike throwing away the range finder because I cannot revert this (imagine a setting where I play around with different tangential directions). For me, different number of tangential directions all need a dedicated reductor. For classical ERA, it is a nice feature that one can have caching for different tangential directions. Because this reductor cannot be applied to very large problems anyway, caching more also does not hurt in most cases. The inconsistency in API is warranted IMO (I though a lot about this as I remember). We can of course discuss this.

I think that if we address both points, then both ERA versions would share a lot more code, making it easier to see where the differences lie. Also maintaining the code would be easier ... - What do you think?

I agree that the architecture could be improved a bit more. However, I also see that this feature is open for a long time now and I think the RandomizedERA class is quite short already. A RandomizedSVD class is what I would have preferred when I started implementing. I am a bit worried that a refactor of randomized_svd will trigger a lot of other design discussions as it raises questions about consistency with qr_svd and orth.

[sdrave]

A RandomizedSVD class is what I would have preferred when I started implementing. I am a bit worried that a refactor of randomized_svd will trigger a lot of other design discussions as it raises questions about consistency with qr_svd and orth.

Let me see if I can come up with a reasonable RandomizedSVD class by the end of the week. I don't think we need a larger debate for doing so. I'd keep randomized_svd as frontend for the user, so there would be no breaking change.

BTW, have you considered using oversampling in your algorithm, @artpelling?

[sdrave]

@artpelling, I have implemented a RandomizedSVD class in #2566. Please take a look if that would be sufficient for your needs. I think the only notable difference to your code should be that it uses qr_svd instead of SciPy for the low-rank SVD. In case this really is a bottleneck for you, we could consider making the SVD alg configurable and adding a scipy_svd method to the svd_va module.

Before merging #2566 some documentation would be needed. I'd also like to think about adding more tolerance parameters to compute_svd. However, I'm not entirely sure what the exact definition of atol and rtol would be for a randomized algorithm. Apart from that, I think that #2566 is mostly good to go.

[artpelling]

implemented a RandomizedSVD class in #2566. Please take a look if that would be sufficient for your needs. I think the only notable difference to your code should be that it uses qr_svd instead of SciPy for the low-rank SVD. In case this really is a bottleneck for you, we could consider making the SVD alg configurable and adding a scipy_svd method to the svd_va module.

Thanks! It looks quite good. I would need to run some tests but I think I would need a SciPy SVD. It would be okay for me if pyMOR did not support that natively, as long as I can monkey patch RandomizedSVD myself easily.

Regarding the tolerances, I am not sure and they are intertwined with the error estimator in the range finder. I suppose it would make most sense to have atol and rtol refer to the expected values of the "normal" quantitites. We could also look here for some inspiration https://github.com/agillman20/librla/blob/main/distrib/python/librla.py

[sdrave]

Thanks! It looks quite good. I would need to run some tests but I think I would need a SciPy SVD. It would be okay for me if pyMOR did not support that natively, as long as I can monkey patch RandomizedSVD myself easily.

I just added a scipy_svd to the svd_va module in #2566. You can now select the used algorithm in RandomizedSVD.__init__.

Regarding the tolerances, I am not sure and they are intertwined with the error estimator in the range finder. I suppose it would make most sense to have atol and rtol refer to the expected values of the "normal" quantitites.

I agree that everything else might get messy. I have added the standard tolerance parameters.

From my perspective, #2566 is now feature complete. I will add documentation/more tests after you had a chance to check if the design is sufficient for you, @artpelling.

[artpelling]

@sdrave I've just rebased on main but have a paper deadline. I can start integrating the new class in August. Feel free to move ahead without me though :)

I should definitely be able to have this complete by the hard freeze 👍