rand_la refactoring

[artpelling]

Main changes:

Based off @sdrave's randomness PR #1736 , I have restructured the rand_la module.

The main contribution is a RandomizedRangeFinder class that can be used to calculate an approximate basis for the range of an Operator. IMO a class is warranted here, because it comes with an estimate_error method which implements the error estimator from @andreasbuhr which I found to be more accurate than the classical one from Sec 4.4 in Halko et. al. (I need to conduct more tests, however).

Moreover, the constructed bases are cached such that one can always come back later and enlarge the range approximation without having to recalculate from scratch (I would like such a functionality with the data-driven reductors such as ERA and think this functionality should be handled by the range approximator and not the reductor).

In addition to this, the functionalities of rrf and adaptive_rrf are now combined in one find_range method where I can pass basis_size and tol arguments. A basis of at least size basis_size will be constructed with error smaller than tol if tol is some float. For now, I have rewritten both functions to use the class and added a Deprecated decorator (I think thats not how its done 🥴). I would be happy to remove them completely as well.

Other features:

• chooseable block_size for adaptive range finder
• the blocksize can be doubled in each adaptive range finding step by setting increase_block to True. This reduces the number of error estimations. The basis_size is then deflated with a binary search to give the smallest basis for given tolerance.
• the adaptive_rrf did not implement subspace iterations. It is now possible to do this with the class but it is lacking theoretic foundation for now.
• nicer logging

TODO:

• refactor SVD and GHEP methods
• prove error bound for subspace iterations
• better names for attributes and methods
• docstrings
• unit tests

[artpelling]

I would like to get some feedback on the current state! Any suggestions regarding code style, variable names, efficiency are greatly welcome. Also anything regarding the error bound with subspace iterations would be helpful.

I also suspect that the caching does not work correctly. When I add a print statement in _c_est, it gets printed multiple times.

[renefritze]

I've changed the PR target to the branch of #1736 to make the diff smaller. We'll change that back to main once that other PR has landed.

[sdrave]

Overall I like the refactoring. I left several comments. Also, before merging I would like to check, if the error bound is still true in case of power iterations (should be fine, I think).

[artpelling]

I dislike having the complex (fka iscomplex) arguments in the randomized methods. IMO this should depend on the operator. Does anyone know a nice way to check the operator for that?

[sdrave]

@artpelling, sorry, was on vacation. Will do proper review by the end of the week.

[sdrave]

I dislike having the complex (fka iscomplex) arguments in the randomized methods.

I dislike it as well.

IMO this should depend on the operator. Does anyone know a nice way to check the operator for that?

No. The deeper problem is that VectorSpaces in pyMOR actually do note care about their underlying field. You can take a real VectorArray and multiply it by 1j, and suddenly it's complex. Even for backends that don't support complex numbers. If you search the issues/PRs, you will find lengthy discussions why we did this. The consequence is, that neither op.source or op.range will give you any information about the realness of the Operator. We could add an attribute for that or even a RuleTable-based algorithms. But then you have to decide whether or not Operators are 'real' or 'unknown' by default. Both would have its issues, but I am open to debate this. However, as long as this is the only place where we run into this problem, it might be better just to live with the current state.

[artpelling]

@sdrave I've just rebased again. How do we proceed?

[sdrave]

I did a small regression test using

import numpy as np
import scipy.linalg as spla

from pymor.algorithms.rand_la import adaptive_rrf
from pymor.operators.numpy import NumpyMatrixOperator


def random_svd_matrix(m, n, cond):
    rng = np.random.default_rng(0)
    U = rng.normal(size=(n, m))
    U = spla.qr(U, mode='economic')[0]
    S = np.diag(np.geomspace(1/cond, 1, m))
    V = rng.normal(size=(m, m))
    V = spla.qr(V)[0]
    return V @ S @ U.T


A = random_svd_matrix(100, 100, 1e20)
range = adaptive_rrf(NumpyMatrixOperator(A), tol=1e-9).to_numpy().T
B = range @ range.T @ A
print(spla.svdvals(A - B)[0])

While the results seem to be correct for main, with the changes the algorithm fails completely. Do you want to take a look?

[artpelling]

A = random_svd_matrix(100, 100, 1e20)
range = adaptive_rrf(NumpyMatrixOperator(A), tol=1e-9).to_numpy().T
B = range @ range.T @ A
print(spla.svdvals(A - B)[0])

While the results seem to be correct for main, with the changes the algorithm fails completely. Do you want to take a look?

It must be related to the recent changes. Running this for 1807d2d works..

[sdrave]

@artpelling, just pushed a fix which then gives the same results as 1807d2d if I checked correctly. However, both yield estimated errors of nan for basis dimension 40 onwards. main does not have such issues. If I print maxnorm / np.sqrt(2. * lambda_min) * erfinv(testfail**(1. / num_testvecs)) I get nice estimates for all iterations. Also, the estimates seem to be different, already for the first iterations.

[artpelling]

These lines look a bit fishy.. @sdrave What are you computing here?

[sdrave]

The norms of the images of the test vectors projected to the orthogonal complement of the current range basis. Why fishy?

[artpelling]

Ah sorry I misread the brackets, I thought you were squaring the np.sqrt. Nonetheless, I think we should add an absolute value here to avoid negative numbers inside of the root.

[artpelling]

OK it seems like you copied the error with the brackets from my code. Sorry about that. Fixing it, will get rid of the nans. I noticed two things:

1. the error estimator does not go below about 3e-7, which is why for the tolerance given, it never stops.
2. the old version uses atol=0 and rtol=0 in gram_schmidt opposed to the default values in the new version.

[sdrave]

With the changes I just pushed, I now get the exact same convergence history as for main -- until the estimated error reaches half machine precision. I believe that the square root of squared differences expression here is the culprit. So the question is how to proceed here.

We introduced this new approach for estimating the error in order to be able to estimate the error for smaller basis sizes without having to recompute the projection onto the orthogonal complement. Personally, I am not too much interested in that feature. As long as the number of test vectors is not increased, these estimates could simply be returned as the convergence history of the algorithm. It would also be possible to later increase the number of test vectors or the basis size. Only going back would be the problem.

BTW, I don't think that taking the absolute value of the difference makes sense. When the difference is zero, then obviously the estimate isn't reliable anymore and an error should be thrown.

[sdrave]

@artpelling, what are your thoughts?

[sdrave]

@artpelling, any ideas how to proceed? In case you don't care that much, I could have another shot at refactoring the module. However, I would probably only allow the number of basis/test vectors to increase to avoid the numercial issues described above.

[artpelling]

I have to admit that this entire PR seems quite far away at the moment. It would take me some time to get back into it (which I don't have in the next couple of weeks). Still, I would be extremely happy if the module could be refactored because I actually need to use it.

So I am very happy for you (or anyone) to take the lead on this. In May, I should have some resources to talk more about the refactoring in case there are some issues.

[sdrave]

Ok, will give it a shot.