Fix cholesky interface with singularities #168
Conversation
Codecov ReportAll modified and coverable lines are covered by tests ✅
Additional details and impacted files@@ Coverage Diff @@
## main #168 +/- ##
==========================================
+ Coverage 92.69% 92.70% +0.01%
==========================================
Files 17 17
Lines 1861 1864 +3
==========================================
+ Hits 1725 1728 +3
Misses 136 136 ☔ View full report in Codecov by Sentry. |
| @@ -43,7 +45,8 @@ cholesky_jacobimatrix(w::Function, P) = cholesky_jacobimatrix(w.(axes(P,1)), P) | |||
|
|
|||
| function cholesky_jacobimatrix(w::AbstractQuasiVector, P) | |||
| Q = normalized(P) | |||
| W = Symmetric(Q \ (w .* Q)) # Compute weight multiplication via Clenshaw | |||
| w_P = orthogonalityweight(P) | |||
| W = Symmetric(Q \ ((w ./ w_P) .* Q)) # Compute weight multiplication via Clenshaw | |||
There was a problem hiding this comment.
Am I missing something - doesn't this break the package?
The convention is that the supplied w is the modification, meaning that as the doc string says we "return the Jacobi matrix X associated to a quasi-matrix of polynomials orthogonal with respect to w(x) w_p(x) where w_p(x) is the weight of the polynomials in P".
If you divide by the orthogonality weight first then instead of the modification, you are supplying the target weight. This will break SemiclassicalOPs among other things.
I am surprised the tests didn't catch this but I guess it's because they are all Legendre based so the w_P term is a constant and thus irrelevant.
There was a problem hiding this comment.
We can change this convention of course but it's a breaking change and needs adjustments downstream.
There was a problem hiding this comment.
The convention should have been the same as Lanczos polynomial. Now it is.
if tests pass it’s by definition not breaking. Downstream packages were therefore using the routine wrong.
There was a problem hiding this comment.
Ok, technically true due to test design, but at the very least the docstring is wrong now and we should fix it downstream since it's our downstream. I'll just fix this in the speed improvements PR then if you think it isn't urgent
There was a problem hiding this comment.
The old version was just broken, eg
Is consistent with my change.
so yes just fix the docstring and down stream packages
No description provided.