Fix cholesky interface with singularities (#168)

* Fix cholesky interface with singularities

* Update test_choleskyQR.jl

* v0.12.3

Author: dlfivefifty

Project.toml

@@ -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"

examples/marchenkopasturops.jl

@@ -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)
@@ -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])

src/choleskyQR.jl

@@ -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)
 
 
 """
@@ -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
     return cholesky_jacobimatrix(W, Q)
 end
 function cholesky_jacobimatrix(W::AbstractMatrix, Q)

test/test_choleskyQR.jl

@@ -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!
 
 
@@ -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
@@ -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
 
@@ -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