WIP: Implement Cholesky-based OPs #63
Conversation
Codecov Report
@@ Coverage Diff @@
## main #63 +/- ##
==========================================
+ Coverage 89.11% 89.43% +0.32%
==========================================
Files 13 14 +1
Lines 1552 1590 +38
==========================================
+ Hits 1383 1422 +39
+ Misses 169 168 -1
Continue to review full report at Codecov.
|
|
Updated now to work as follows P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
x = axes(P,1)
J = jacobimatrix(P)
wf(x) = x*(1-x^2)*exp(-x^2)
# compute Jacobi matrix via cholesky
Jchol = cholesky_jacobimatrix(wf, P)
# compute Jacobi matrix via Lanczos
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
# Comparison with Lanczos
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]and with arbitrary OP basis, at least arbitrary Jacobi basis. When P is not Legendre, W is assumed to be the weight modification operator rather than simply the weight multiplication operator, i.e. it is generally the (new weight)/(old weight) multiplication operator. |
|
Shouldn't this just replace lanczos? |
|
I think it's too early to say? If it ends up being better, yes, but there's still a lot to figure out so I'm keeping them both around to compare. |
|
If I understand correctly, the inverse Cholesky factor of the polynomial weight modifier is the connection coefficient matrix? |
|
Yes, that's right @MikaelSlevinsky. There is also a slightly different version that would do this with a different decomposition but for now this is the one I'm going to mainly work on. Via Clenshaw approach approximately-banded multiplication operators it also works for non-polynomial weight modification but I'm still in the process of exploring exactly what the pros and cons of this approach are. |
|
I think this is really cool! You might have to spell out the cons of this connection to linear algebra -- I only see pros. For example,
|
|
Yes, it's really cool isn't it? I am very excited about this as an approach overall, I just don't want to promise too much before I've had the time to properly stress test it. We only very recently figured out this connection (as well as a QR analogue which I might implement at some point which in my tests so far is also performing quite well), following on from Sheehan's blogpost about the Lanczos approach a while back. There is some previous notice of this connection in the literature but it seems it never was used for anything and stayed a side-observation but it seems that decomposition methods can do a lot of good work for us. From what I can tell it's more stable than Lanczos as long as both Lanczos and Cholesky approach are started from the same initial basis. It also seems that like Lanczos, the closer the initial weight is to the final one the better accuracy we get. I'd say that's all fairly expected. I am working on the analysis part, this is sort of my in-practice sanity check until it's fully ready. It's can all be done adaptively and banded so ideally as Sheehan said this can replace Lanczos in most scenarios. I'm very interested in your first point actually, maybe we should talk more about what you have in mind? |
|
Hmm, I overlooked something and now it's not quite that simple. I'll describe the polynomial case then outline how to modify it for rational w(x). Polynomial w(x). Let It follows that if Rational w(x). Let The relationship between the inverse connection coefficients and the symmetric and banded multiplication matrices (In both cases, the modified Jacobi matrix |
|
I think I am already making use of that? The polynomial case as you describe it is more or less exactly how I first found this relationship and once one has that it's straightforward to generalise this for non-Legendre as well (although it's particularly pretty for Legendre). The case of rational weight is then somewhat analogous to the problem of starting with a non-Legendre basis. Ideally the weights would cancel such that u(x)/v(x) is a polynomial, e.g. mapping between different Jacobi basis already works: P = Normalized(Jacobi(1,1))
# modification by (1-x)^2 should map from Jacobi(1,1) to Jacobi(3,1)
wf(x) = (1-x)^2
# compute Jacobi(3,1) jacobimatrix via cholesky
Jchol = cholesky_jacobimatrix(wf, P)
# compare with canonical result
Q = Normalized(Jacobi(3,1))
J = jacobimatrix(Q)
@test J[1:500,1:500] ≈ Jchol[1:500,1:500] and in those cases we can just take modifier P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
x = axes(P,1)
J = jacobimatrix(P)
wf(x) = exp(x^2)/(1.4-x)^2
# compute Jacobi matrix via cholesky
Jchol = cholesky_jacobimatrix(wf, P)
# compute Jacobi matrix via Lanczos
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
# Comparison with Lanczos
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]So far so good, but this is of course a bad idea if the weight has singularities on the domain boundary as it will be terribly slow - for more well-behaved things it's perfectly fine and quite fast. This is why there is potential for use in semiclassical bases. I've been wondering if we can do better for weights with singularities but it isn't obvious how to. Computing |
|
FYI |
|
I guess the rational case can be done in a few steps. Assume that
Is the full connection problem not just If this works, I think it's better to do rational modifications like this rather than approximating the rational as a polynomial. It would still be nice to do the full rational case in one sweep. |
|
I guess solving |
|
No, not quite (if I understand you correctly) but I've been thinking about this for a bit now as well and I also think we can do it. Start in the basis with u(x) as its weight (if it's not a classical one, then just move there first via Lanczos or Cholesky - I think you were also thinking along these lines?). Compute the |
|
Ah, there you go, you saw the issue already. Yes, that is entirely doable. A possible problem from the adaptive perspective is that I think the |
|
@dlfivefifty any opinions on For finite matrices it's equivalent to Cholesky if we apply permutation matrices in the right places. Not really an option for infinite case. If we can do this in a fast and stable way then rational polynomial weights are within reach with this approach. Basic examples already work for me. |
|
For Chebyshev they’ll be toeplitz so coming back from ∞ is Wiener hopf legendre case needs some thought |
|
Processing numerator and denominator modifications separately is sensitive to common-root problems, where neither multiplication operator is positive-definite. That means either pre-processing to remove common roots or thinking of a combined modification. |
|
The full rational case can also be done with an infinite QL decomposition. Suppose
The advantage in working with both U and V in the same step is that common-root problems would likely only cause numerical zeros in some outermost bands of some banded matrices. The tricky part is finding the QL decomposition. Maybe it can be done by doubling the dimensions until convergence sets in? |
|
QL decompositions of infinite matrices are already implemented in InfiniteLinearAlgebra.jl for perturbations of Toeplitz operators, perhaps there is a way to make that work for asymptotically Toeplitz operators as well. Uniqueness of the decomposition may also be a concern. Just scaling it to convergence is probably an option as you say, at least my basic testing suggests that it would work, but it would be slow to use that for anything adaptive since the decomposition would have to be re-computed in full instead of extended. Might be worth the cost in the end for rational weights, almost certainly still better than approximating it as a polynomial, but I'll give it some thought if we can do more. |
|
Since |
Starting to properly work on this now @dlfivefifty. Might also be of interest to @MarcoFasondini. Here's an example of what works so far:
Since it's intended to be used for normalised OPs, I'm using symmetry to store the Jacobi operator's bands as a (2,ℵ₀) cached array. It doesn't actually at all need the shift that I use in the tests - those only are there at the moment because I want to eventually compare to SemiclassicalOPs.jl.
There are a few limitations to fully exploring the limits of this approach as it stands:
1) Adaptive cholesky seems to currently not want to play well with the Clenshaw type, meaning that we can't just compute the weight multiplication matrix via the standard way of writing something likeW = P \ (w.(x) .* P). I also cannot see an acceptably fast way to convert from Clenshaw to something more useful. Any ideas @dlfivefifty?Nevermind, this actually works already, it just didn't know that W is symmetric. Easily fixed via
Symmetric(W). That's great, I will update with more complicated examples soon.J = jacobimatrix(P)to a symmetric Banded matrix at the moment in some places. This works fine as a proof of concept for the approach but we inevitably will want to not have that unnecessary conversion in there.But the basic idea of the method seems to work. If we work on the above two points we can see where it starts to run into problems as well as see if it performs better or worse than Lanczos in accuracy / speed.