-
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
shifted Cholesky QR algorithm #2177
base: main
Are you sure you want to change the base?
Conversation
d772ddf
to
d3ce675
Compare
I thought a bit about how to add an where the (upper) Cholesky factor where |
I guess the important part in terms of the (RR)QR interface is that the first |
I've pushed a sketch. This is the best I could come up with so far. Not sure, if the copying and appending of |
d55f5d3
to
349e6c6
Compare
89bbf1f
to
d19d3c5
Compare
I am not sure, if the logger should produce a warning, if |
Adding more logs might be nice, I guess. It shouldn't have too much, but it will probably have less than |
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.
To make this easier to read, I'd probably:
- merge
_shifted_choslesky
into the main function - replace
_CholeskyShifter
by_compute_shift(A, X, product, product_norm, check_finite)
- replace
_ProductNorm
byproduct_norm = None
at the beginning of the main loop andif product_norm is None: product_norm = 1 if product is None else np.sqrt(np.abs(eigs(self.product, k=1)[0][0]))
30ad697
to
0fdb3de
Compare
I went ahead and put everything into the main function. |
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.
Just minor comments for now. Would be good to have tests (similar to pymortests/algorithms/gram_schmidt
).
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 like the code much better now!
Tests are still missing. These should also include testing the offset
argument.
Docs build is failing with
/builds/pymor/pymor/docs/source/autoapi/pymor/algorithms/cholesky_qr/index.rst:19: WARNING: Bullet list ends without a blank line; unexpected unindent.
Further, the algorithm does not seem to be tested for complex arrays:
.../pymor/src/pymor/algorithms/cholesky_qr.py:112: ComplexWarning: Casting complex values to real discards the imaginary part
Rinv = spla.lapack.dtrtri(Rx)[0].T
00:47 shifted_chol_qr: Iteration 2
.../src/pymor/algorithms/cholesky_qr.py:124: ComplexWarning: Casting complex values to real discards the imaginary part
spla.blas.dtrmm(1, Rx, Ri, overwrite_b=True)
00:47 shifted_chol_qr: Iteration 3
.../pymor/src/pymor/algorithms/cholesky_qr.py:140: ComplexWarning: Casting complex values to real discards the imaginary part
R[:offset, offset:] = Bi
.../pymor/src/pymor/algorithms/cholesky_qr.py:141: ComplexWarning: Casting complex values to real discards the imaginary part
R[offset:, offset:] = Ri
6beaab3
to
3dd7d9c
Compare
@artpelling, any ideas why the your tests seem to be hanging in CI. Have you tried locally? |
@artpelling, @pmli, tests are really hanging because |
Or return |
You mean if |
Yes, exactly. |
It seems that the tests are failing because of an infinite loop when a zero vector array is given. I guess shifted Cholesky can cheaply check if the input is zero and exit early. |
Ah thats why. I pushed a fix thank you! |
It looks like tests fail now because from pymor.algorithms.chol_qr import shifted_chol_qr
from pymor.vectorarrays.numpy import NumpyVectorSpace
A = NumpyVectorSpace(5).ones(2)
print(A)
Q, R = shifted_chol_qr(A)
print(Q)
print(R) produces
Fixing this doesn't seem simple to me. Also, |
Indeed, I looked into it but am not sure how to deal with this. The question is what we want the algorithm to do in this case. I might need some help here to develop a strategy.
This one could be fixed by prescribing a minimum shift. |
Based on the discussion in the previous community meeting, and since the "correct" approach is unknown (should use unpublished pivoted Cholesky, I guess), shifted Cholesky QR should probably raise an exception if the computed |
Regarding testing, I'd suggest to call |
I hope I fixed the tests now as discussed. @maxbindhak noticed that the computation of the 2-norm is the bottleneck of the method as it stands. In a discussion with @drittelhacker, we realized that we can speed up this computation by exploiting the Hermitian structure of the gramian (using |
I just fixed another error in the shift computation when @sdrave This PR is ready now from my side. |
|
||
if offset == len(A): | ||
return A, np.eye(len(A)) | ||
elif A.sup_norm().max() < 1e-16: |
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 don't think it's good practice to dismiss A
as zero based on some absolute norm. In general, we don't know anything about the scaling of A
. Is this needed for something?
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.
It serves to detect "zero" input, as well as inputs with zero length or dimension. I did this mainly to make the tests pass (took me quite some time).
@sdrave I've implemented all changes, thanks for the review! Unfortunately, one test now fails in |
Not sure why it was not failing before. If I run this locally, I get examples with assert np.all(almost_equal(U, onb.lincomb(p.apply2(onb, U).T), rtol=1e-13)) Raising the tolerance to 1e-11 fixes this. |
Thanks for investigating and finding the issue. Would it be ok to catch the
I was also puzzled by this. Only explanation would be that the old code didn't enter
Yes I had this, too. I will raise the tolerance. |
Let's first see if we can quickly fix the actual issue. If not, I'd say it's fine to proceed this way. |
Adds the Cholesky based QR algorithm from this paper with an oblique inner product. Shifts are applied if necessary (see Algorithm 4.1 in the paper).