-
Notifications
You must be signed in to change notification settings - Fork 97
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
rand_la
refactoring
#1753
base: main
Are you sure you want to change the base?
rand_la
refactoring
#1753
Conversation
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 |
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. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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).
59376b0
to
24b67f5
Compare
I dislike having the |
38eec03
to
a2641e3
Compare
@artpelling, sorry, was on vacation. Will do proper review by the end of the week. |
I dislike it as well.
No. The deeper problem is that |
5a3bcf6
to
d21123f
Compare
3fabc0d
to
2b44de4
Compare
@sdrave I've just rebased again. How do we proceed? |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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?
It must be related to the recent changes. Running this for 1807d2d works.. |
@artpelling, just pushed a fix which then gives the same results as 1807d2d if I checked correctly. However, both yield estimated errors of |
src/pymor/algorithms/rand_la.py
Outdated
norms = np.sqrt( | ||
np.abs(norms**2) - spla.norm(self._projection_coeffs[:num_testvecs, :complement_basis_size], axis=1)**2 | ||
) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
These lines look a bit fishy.. @sdrave What are you computing here?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The norms of the images of the test vectors projected to the orthogonal complement of the current range basis. Why fishy?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
OK it seems like you copied the error with the brackets from my code. Sorry about that. Fixing it, will get rid of the nan
s. I noticed two things:
- the error estimator does not go below about
3e-7
, which is why for the tolerance given, it never stops. - the old version uses
atol=0
andrtol=0
ingram_schmidt
opposed to the default values in the new version.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@artpelling, what are your thoughts?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ok, will give it a shot.
to reproduce behavior of adaptive_rrf in main
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 therange
of anOperator
. IMO a class is warranted here, because it comes with anestimate_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
andadaptive_rrf
are now combined in onefind_range
method where I can passbasis_size
andtol
arguments. A basis of at least sizebasis_size
will be constructed with error smaller thantol
if tol is some float. For now, I have rewritten both functions to use the class and added aDeprecated
decorator(I think thats not how its done 🥴). I would be happy to remove them completely as well.Other features:
chooseableblock_size
for adaptive range finderthe blocksize can be doubled in each adaptive range finding step by settingincrease_block
toTrue
. This reduces the number of error estimations. Thebasis_size
is then deflated with a binary search to give the smallest basis for given tolerance.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.TODO: