WIP: Replace Lanzos with Cholesky

Project.toml

@@ -14,17 +14,17 @@ InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
+SingularIntegrals = "d7440221-8b5e-42fc-909c-0567823f424a"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-ArrayLayouts = "0.8"
+ArrayLayouts = "1"
 BandedMatrices = "0.17"
-ClassicalOrthogonalPolynomials = "0.6, 0.7"
 ContinuumArrays = "0.10, 0.11, 0.12"
-FillArrays = "0.13"
+FillArrays = "1"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.12.9"
-LazyArrays = "0.22.14"
+LazyArrays = "1"
 QuasiArrays = "0.9"
 SpecialFunctions = "1.0, 2"
 julia = "1.7"

src/SemiclassicalOrthogonalPolynomials.jl

@@ -7,12 +7,13 @@ import Base: getindex, axes, size, \, /, *, +, -, summary, ==, copy, sum, unsafe
 
 import ArrayLayouts: MemoryLayout, ldiv, diagonaldata, subdiagonaldata, supdiagonaldata
 import BandedMatrices: bandwidths, AbstractBandedMatrix, BandedLayout, _BandedMatrix
-import LazyArrays: resizedata!, paddeddata, CachedVector, CachedMatrix, CachedAbstractVector, LazyMatrix, LazyVector, arguments, ApplyLayout, colsupport, AbstractCachedVector,
+import LazyArrays: resizedata!, paddeddata, CachedVector, CachedMatrix, CachedAbstractVector, LazyMatrix, LazyVector, arguments, ApplyLayout, colsupport, AbstractCachedVector, ApplyArray,
                     AccumulateAbstractVector, LazyVector, AbstractCachedMatrix, BroadcastLayout
 import ClassicalOrthogonalPolynomials: OrthogonalPolynomial, recurrencecoefficients, jacobimatrix, normalize, _p0, UnitInterval, orthogonalityweight, NormalizedOPLayout,
                                         Bidiagonal, Tridiagonal, SymTridiagonal, symtridiagonalize, normalizationconstant, LanczosPolynomial,
-                                        OrthogonalPolynomialRatio, Weighted, WeightLayout, UnionDomain, oneto, Hilbert, WeightedBasis, HalfWeighted,
-                                        Associated, golubwelsch, associated, AbstractOPLayout, weight
+                                        OrthogonalPolynomialRatio, Weighted, WeightLayout, UnionDomain, oneto, WeightedBasis, HalfWeighted,
+                                        golubwelsch, AbstractOPLayout, weight, cholesky_jacobimatrix, qr_jacobimatrix
+import SingularIntegrals: Hilbert, Associated
 import InfiniteArrays: OneToInf, InfUnitRange
 import ContinuumArrays: basis, Weight, @simplify, AbstractBasisLayout, BasisLayout, MappedBasisLayout, grid, plotgrid, _equals, ExpansionLayout
 import FillArrays: SquareEye
@@ -93,20 +94,13 @@ RaisedOP(Q, y::Number) = RaisedOP(Q, OrthogonalPolynomialRatio(Q,y))
 function jacobimatrix(P::RaisedOP{T}) where T
     ℓ = P.ℓ
     X = jacobimatrix(P.Q)
-    a,b = diagonaldata(X), supdiagonaldata(X)
+    a,b = view(X,band(0)), view(X,band(1)) # a,b = diagonaldata(X), supdiagonaldata(X)
     # non-normalized lower diag of Jacobi
     v = Vcat(zero(T),b .* ℓ)
     c = BroadcastVector((ℓ,a,b,sa,v) -> ℓ*a + b - b*ℓ^2 - sa*ℓ + ℓ*v, ℓ, a, b, a[2:∞], v)
     Tridiagonal(c, BroadcastVector((ℓ,a,b,v) -> a - b * ℓ + v, ℓ,a,b,v), b)
 end
 
-
-
-
-
-
-
-
 """
   the bands of the Jacobi matrix
 """
@@ -141,7 +135,6 @@ SemiclassicalJacobiBand{dv}(a,b,ℓ) where dv = SemiclassicalJacobiBand{dv,promo
 copy(r::SemiclassicalJacobiBand) = r # immutable
 
 
-
 """
    SemiclassicalJacobi(t, a, b, c)
 
@@ -176,7 +169,7 @@ function semiclassical_jacobimatrix(t, a, b, c)
     if a == 0 && b == 0 && c == -1
         # for this special case we can generate the Jacobi operator explicitly
         N = (1:∞)
-        α = T.(neg1c_αcfs(t)) # cached coefficients
+        α = neg1c_αcfs(t) # cached coefficients
         A = Vcat((α[1]+1)/2 , -N./(N.*4 .- 2).*α .+ (N.+1)./(N.*4 .+ 2).*α[2:end].+1/2)
         C = -(N)./(N.*4 .- 2)
         B = Vcat((α[1]^2*3-α[1]*α[2]*2-1)/6 , -(N)./(N.*4 .+ 2).*α[2:end]./α)
@@ -186,10 +179,15 @@ function semiclassical_jacobimatrix(t, a, b, c)
     P = jacobi(b, a, UnitInterval{T}())
     iszero(c) && return symtridiagonalize(jacobimatrix(P))
     x = axes(P,1)
-    jacobimatrix(LanczosPolynomial(@.(x^a * (1-x)^b * (t-x)^c), P))
+    if iseven(c) # QR case
+        return qr_jacobimatrix(x->(t.-x).^(c÷2), Normalized(jacobi(a,b,UnitInterval{T}())))
+    elseif isodd(c)
+        return semiclassical_jacobimatrix(SemiclassicalJacobi(t, a, b, c-1), a, b, c)
+    else # if c is not an even integer, use Lanczos for now
+        return jacobimatrix(LanczosPolynomial(@.(x^a * (1-x)^b * (t-x)^c), P))
+    end
 end
 
-
 function symraised_jacobimatrix(Q, y)
     ℓ = OrthogonalPolynomialRatio(Q,y)
     X = jacobimatrix(Q)
@@ -198,9 +196,35 @@ function symraised_jacobimatrix(Q, y)
 end
 
 function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c)
-    if  iszero(a) && iszero(b) && c == -1 # (a,b,c) = (0,0,-1) special case
+    # special cases 
+    if iszero(a) && iszero(b) && c == -1 # (a,b,c) = (0,0,-1) special case
         semiclassical_jacobimatrix(Q.t, a, b, c)
-    elseif a == Q.a+1 && b == Q.b && c == Q.c  # raising by 1
+    elseif iszero(c) # classical Jacobi polynomial special case
+        jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c))
+    elseif a == Q.a && b == Q.b && c == Q.c # same basis
+        jacobimatrix(Q)
+    end
+
+    # QR decomposition approach for raising
+    Δa = a-Q.a
+    Δb = b-Q.b
+    Δc = c-Q.c
+    if iseven(a-Q.a) && iseven(b-Q.b) && iseven(c-Q.c) && Δa≥0 && Δb≥0 && Δc≥0
+        M = Diagonal(Ones(∞))
+        if !iszero(Δa)
+            M = ApplyArray(*,M,(Q.X)^((a-Q.a)÷2))
+        end
+        if !iszero(Δb)
+            M = ApplyArray(*,M,(I-Q.X)^((b-Q.b)÷2))
+        end
+        if !iszero(Δc)
+            M = ApplyArray(*,M,(Q.t*I-Q.X)^((c-Q.c)÷2))
+        end
+        return qr_jacobimatrix(M,Q,false)
+    end
+
+    # Fallback from decomposition methods
+    if a == Q.a+1 && b == Q.b && c == Q.c  # raising by 1
         symraised_jacobimatrix(Q, 0)
     elseif a == Q.a && b == Q.b+1 && c == Q.c
         symraised_jacobimatrix(Q, 1)
@@ -212,20 +236,18 @@ function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c)
         symlowered_jacobimatrix(Q, :a)
     elseif a == Q.a && b == Q.b-1 && c == Q.c
         symlowered_jacobimatrix(Q, :b)
-    elseif a > Q.a  # iterative raising
+    elseif a > Q.a  # iterative raising by 1
         semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a+1, Q.b, Q.c, Q), a, b, c)
     elseif b > Q.b
         semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b+1, Q.c, Q), a, b, c)
     elseif c > Q.c
         semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c+1, Q), a, b, c)
-    elseif b < Q.b  # iterative lowering
+    elseif b < Q.b  # iterative lowering by 1
         semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b-1, Q.c, Q), a, b, c)
     elseif c < Q.c
         semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c-1, Q), a, b, c)
     elseif a < Q.a 
         semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a-1, Q.b, Q.c, Q), a, b, c)
-    elseif a == Q.a && b == Q.b && c == Q.c # same basis
-        jacobimatrix(Q)
     else
         error("Not Implemented")
     end
@@ -291,6 +313,46 @@ function semijacobi_ldiv(P::SemiclassicalJacobi, Q)
     (P \ R) * _p0(R̃) * (R̃ \ Q)
 end
 
+# returns conversion operator from SemiclassicalJacobi P to SemiclassicalJacobi Q.
+function semijacobi_ldiv(Q::SemiclassicalJacobi,P::SemiclassicalJacobi)
+    @assert Q.t ≈ P.t
+    # For polynomial modifications, use Cholesky. Use Lanzos otherwise.
+    M = Diagonal(Ones(∞))
+    (Q == P) && return M
+    Δa = Q.a-P.a
+    Δb = Q.b-P.b
+    Δc = Q.c-P.c
+    if iseven(Δa) && iseven(Δb) && iseven(Δc)
+        if !iszero(Q.a-P.a)
+            M = ApplyArray(*,M,(P.X)^((Q.a-P.a)÷2))
+        end
+        if !iszero(Q.b-P.b)
+            M = ApplyArray(*,M,(I-P.X)^((Q.b-P.b)÷2))
+        end
+        if !iszero(Q.c-P.c)
+            M = ApplyArray(*,M,(Q.t*I-P.X)^((Q.c-P.c)÷2))
+        end
+        K = qr(M).R
+        return ApplyArray(*, Diagonal(sign.(view(K,band(0))).*Fill(abs.(1/K[1]),∞)), K) # match our normalization choice P_0(x) = 1
+    elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc)
+        if !iszero(Q.a-P.a)
+            M = ApplyArray(*,M,(P.X)^(Q.a-P.a))
+        end
+        if !iszero(Q.b-P.b)
+            M = ApplyArray(*,M,(I-P.X)^(Q.b-P.b))
+        end
+        if !iszero(Q.c-P.c)
+            M = ApplyArray(*,M,(Q.t*I-P.X)^(Q.c-P.c))
+        end
+        K = cholesky(Symmetric(M)).U
+        return ApplyArray(*, Diagonal(Fill(1/K[1],∞)), K) # match our normalization choice P_0(x) = 1
+    else # fallback option for non-integer weight modification
+        R = SemiclassicalJacobi(P.t, mod(P.a,-1), mod(P.b,-1), mod(P.c,-1))
+        R̃ = toclassical(R)
+        return (P \ R) * _p0(R̃) * (R̃ \ Q)
+    end
+end
+
 struct SemiclassicalJacobiLayout <: AbstractOPLayout end
 MemoryLayout(::Type{<:SemiclassicalJacobi}) = SemiclassicalJacobiLayout()
 
@@ -322,6 +384,7 @@ function copy(L::Ldiv{SemiclassicalJacobiLayout,SemiclassicalJacobiLayout})
     inv(M_Q) * L' * M_P
 end
 
+\(A::SemiclassicalJacobi, B::SemiclassicalJacobi) = semijacobi_ldiv(A, B)
 \(A::LanczosPolynomial, B::SemiclassicalJacobi) = semijacobi_ldiv(A, B)
 \(A::SemiclassicalJacobi, B::LanczosPolynomial) = semijacobi_ldiv(A, B)
 function \(w_A::WeightedSemiclassicalJacobi, w_B::WeightedSemiclassicalJacobi)
@@ -455,7 +518,27 @@ function SemiclassicalJacobiFamily{T}(data::Vector, t, a, b, c) where T
 end
 
 SemiclassicalJacobiFamily(t, a, b, c) = SemiclassicalJacobiFamily{float(promote_type(typeof(t),eltype(a),eltype(b),eltype(c)))}(t, a, b, c)
-SemiclassicalJacobiFamily{T}(t, a, b, c) where T = SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c))], t, a, b, c)
+function SemiclassicalJacobiFamily{T}(t, a, b, c) where T
+    ℓa = length(a)
+    ℓb = length(b)
+    ℓc = length(c)
+    # We need to start with a hierarchy containing two entries
+    if  ℓa > 1 && ℓb > 1 && ℓc > 1 
+        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, a[2], b[2], c[2])], t, a, b, c)
+    elseif ℓb > 1 && ℓc > 1 
+        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, first(a), b[2], c[2])], t, a, b, c)
+    elseif ℓa > 1 && ℓc > 1 
+        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, a[2], first(b), c[2])], t, a, b, c)
+    elseif ℓa > 1 && ℓb > 1 
+        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, a[2], b[2], first(c))], t, a, b, c)
+    elseif ℓc > 1 
+        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, first(a), first(b), c[2])], t, a, b, c)
+    elseif ℓa > 1 
+        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, a[2], first(b), first(c))], t, a, b, c)
+    else #if ℓb > 1 
+        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, first(a), b[2], first(c))], t, a, b, c)
+    end
+end
 
 Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Number, b::Number, c::Number) = SemiclassicalJacobi(t, a, b, c)
 Base.broadcasted(::Type{SemiclassicalJacobi{T}}, t::Number, a::Number, b::Number, c::Number) where T = SemiclassicalJacobi{T}(t, a, b, c)
@@ -471,7 +554,7 @@ _broadcast_getindex(a::Number,k) = a
 function LazyArrays.cache_filldata!(P::SemiclassicalJacobiFamily, inds::AbstractUnitRange)
     t,a,b,c = P.t,P.a,P.b,P.c
     for k in inds
-        P.data[k] = SemiclassicalJacobi(t, _broadcast_getindex(a,k), _broadcast_getindex(b,k), _broadcast_getindex(c,k), P.data[k-1])
+        P.data[k] = SemiclassicalJacobi(t, _broadcast_getindex(a,k), _broadcast_getindex(b,k), _broadcast_getindex(c,k), P.data[k-2])
     end
     P
 end

src/lowering.jl

@@ -1,5 +1,5 @@
 ###########
-# Generic methods for obtaining Jacobi matrices with one of the a, b and c parameters lowered by 1.
+# Generic fallback methods for obtaining Jacobi matrices with one of the a, b and c parameters lowered by 1.
 # passing symbols :a, :b or :c into lowindex determines which parameter is lowered
 ######
 function initialα_gen(P::SemiclassicalJacobi, lowindex::Symbol)
@@ -163,8 +163,8 @@ end
 
 ###########
 # Methods for the special case of computing SemiclassicalJacobi(t,0,0,-1)
-######
 # As this can be built from SemiclassicalJacobi(t,0,0,0) which is just shifted and normalized Legendre(), we have more explicit methods at our disposal due to explicitly known coefficients for the Legendre bases.
+######
 
 ###
 #   α coefficients implementation, cached, direct and via recurrences

src/twoband.jl

@@ -98,7 +98,6 @@ function jacobimatrix(R::TwoBandJacobi{T}) where T
     # M = (L / L[1,1])' # equal to R.Q \ R.P
 
     Tridiagonal(Interlace(L.dv/L[1,1], (ρ^2-1) * L.ev), Zeros{T}(∞), Interlace((1-ρ^2) * L.dv,L.ev/(-L[1,1])))
-    # Tridiagonal(Interlace(L.dv/L[1,1], (1-ρ^2) * L.ev), Zeros{T}(∞), Interlace((1-ρ^2) * L.dv, L.ev/L[1,1]))
 end
 
 

test/runtests.jl

@@ -117,7 +117,7 @@ end
         @testset "Mass matrix" begin
             @test (T'*(w_T .* T))[1:10,1:10] ≈ sum(w_T)I
             M = W'*(w_W .* W)
-            @test (sum(w_T)*inv(R)'L)[1:10,1:10] ≈ M[1:10,1:10]
+            @test (sum(w_T)*inv(R[1:10,1:10])'L[1:10,1:10]) ≈ M[1:10,1:10]
             @test T[0.1,1:10]' ≈ W[0.1,1:10]' * R[1:10,1:10]
         end
     end
@@ -450,3 +450,4 @@ end
 include("test_derivative.jl")
 include("test_twoband.jl")
 include("test_lowering.jl")
+