Skip to content

Fix cholesky interface with singularities #168

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.

Sign up for GitHub

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

Merged
merged 3 commits into from
Jan 21, 2024
Merged

Fix cholesky interface with singularities #168

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
83749d2
Fix cholesky interface with singularities
dlfivefifty Jan 20, 2024
24491ca
Update test_choleskyQR.jl
dlfivefifty Jan 21, 2024
3ea4841
v0.12.3
dlfivefifty Jan 21, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "ClassicalOrthogonalPolynomials"
uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
authors = ["Sheehan Olver <[email protected]>"]
version = "0.12.2"
version = "0.12.3"

[deps]
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
Expand Down
8 changes: 2 additions & 6 deletions examples/marchenkopasturops.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,12 +4,10 @@ using ClassicalOrthogonalPolynomials, Plots
r = 0.5
lmin, lmax = (1-sqrt(r))^2, (1+sqrt(r))^2
U = chebyshevu(lmin..lmax)
x = axes(U,1)
w = @. 1/(2π) * sqrt((lmax-x)*(x-lmin))/(x*r)

# Q is a quasimatrix such that Q[x,k+1] is equivalent to
# qₖ(x), the k-th orthogonal polynomial wrt to w
Q = LanczosPolynomial(w, U)
Q = OrthogonalPolynomial(x -> 1/(2π) * sqrt((lmax-x)*(x-lmin))/(x*r), U)

# The Jacobi matrix associated with Q, as an ∞×∞ SymTridiagonal
J = jacobimatrix(Q)
Expand All @@ -22,7 +20,5 @@ a,b = 5,10
c,d = sqrt(a/(a+b) * (1-1/(a+b))), sqrt(1/(a+b) * (1-a/(a+b)))
lmin,lmax = (c-d)^2,(c+d)^2
U = chebyshevu(lmin..lmax)
x = axes(U,1)
w = @. (a+b) * sqrt((x-lmin)*(lmax-x))/(2π*x*(1-x))
Q = LanczosPolynomial(w, U)
Q = OrthogonalPolynomial(x -> (a+b) * sqrt((x-lmin)*(lmax-x))/(2π*x*(1-x)), U)
plot(Q[:,1:7])
7 changes: 5 additions & 2 deletions src/choleskyQR.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,9 +26,11 @@ function OrthogonalPolynomial(w::AbstractQuasiVector, P::AbstractQuasiMatrix)
ConvertedOrthogonalPolynomial(w, X, parent(X.dv).U, Q)
end

orthogonalpolynomial(w::AbstractQuasiVector) = OrthogonalPolynomial(w)
orthogonalpolynomial(wP...) = OrthogonalPolynomial(wP...)
orthogonalpolynomial(w::SubQuasiArray) = orthogonalpolynomial(parent(w))[parentindices(w)[1],:]

OrthogonalPolynomial(w::Function, P::AbstractQuasiMatrix) = OrthogonalPolynomial(w.(axes(P,1)), P)
orthogonalpolynomial(w::Function, P::AbstractQuasiMatrix) = orthogonalpolynomial(w.(axes(P,1)), P)


"""
Expand All @@ -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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can change this convention of course but it's a breaking change and needs adjustments downstream.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The old version was just broken, eg

ClassicalOrthogonalPolynomials.jl/src/choleskyQR.jl

Line 22 in 5256e8e

OrthogonalPolynomial(w::AbstractQuasiVector) = OrthogonalPolynomial(w, orthogonalpolynomial(singularities(w)))

Is consistent with my change.

so yes just fix the docstring and down stream packages

return cholesky_jacobimatrix(W, Q)
end
function cholesky_jacobimatrix(W::AbstractMatrix, Q)
Expand Down
15 changes: 12 additions & 3 deletions test/test_choleskyQR.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
using Test, ClassicalOrthogonalPolynomials, BandedMatrices, LinearAlgebra, LazyArrays, ContinuumArrays, LazyBandedMatrices, InfiniteLinearAlgebra
import ClassicalOrthogonalPolynomials: cholesky_jacobimatrix, qr_jacobimatrix
import ClassicalOrthogonalPolynomials: cholesky_jacobimatrix, qr_jacobimatrix, orthogonalpolynomial
import LazyArrays: AbstractCachedMatrix, resizedata!


Expand All @@ -10,7 +10,7 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
P = Normalized(legendre(0..1))
x = axes(P,1)
J = jacobimatrix(P)
wf(x) = x^2*(1-x)
wf = x -> x^2*(1-x)
# compute Jacobi matrix via cholesky
Jchol = cholesky_jacobimatrix(wf, P)
# compute Jacobi matrix via classical recurrence
Expand Down Expand Up @@ -224,7 +224,7 @@ import LazyArrays: AbstractCachedMatrix, resizedata!

@test Q == Q
@test P ≠ Q
@test Q ≠ P
@test Q ≠ P
@test Q == Q̃
@test Q̃ == Q

Expand All @@ -243,4 +243,13 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
# need to fix InfiniteLinearAlgebra to add AdaptiveBandedLayout
@test_broken R[1:10,1:10] isa BandedMatrix
end

@testset "Chebyshev" begin
U = ChebyshevU()
Q = orthogonalpolynomial(x -> (1+x^2)*sqrt(1-x^2), U)
@test bandwidths(Q\U) == (0,2)

Q̃ = OrthogonalPolynomial(x -> (1+x^2)*sqrt(1-x^2), U)
@test jacobimatrix(Q)[1:10,1:10] == jacobimatrix(Q̃)[1:10,1:10]
end
end