WIP: Implement Cholesky-based OPs

src/ClassicalOrthogonalPolynomials.jl

@@ -361,6 +361,9 @@ include("classical/ultraspherical.jl")
 include("classical/laguerre.jl")
 include("classical/fourier.jl")
 include("stieltjes.jl")
+include("decompOPs.jl")
+
+export cholesky_jacobimatrix
 
 
 end # module

src/decompOPs.jl

@@ -0,0 +1,64 @@
+# This currently takes the weight multiplication operator as input.
+# I will probably change this to take the weight function instead.
+function cholesky_jacobimatrix(w, P::OrthogonalPolynomial)
+    W = Symmetric(P \ (w.(axes(P,1)) .* P)) # Compute weight multiplication via Clenshaw
+    bands = CholeskyJacobiBands(W, P) # the cached array only needs to store two bands bc of symmetry
+    return SymTridiagonal(bands[1,:],bands[2,:])
+end
+
+# The generated Jacobi operators are symmetric tridiagonal, so we store their data as two bands
+mutable struct CholeskyJacobiBands{T} <: AbstractCachedMatrix{T}
+    data::Matrix{T}
+    U::UpperTriangular
+    X::SymTridiagonal{T}
+    datasize::Int
+    array
+end
+
+# Computes the initial data for the Jacobi operator bands
+function CholeskyJacobiBands(W::Symmetric{T}, P::OrthogonalPolynomial) where T
+    @assert P isa Normalized
+    U = cholesky(W).U
+    X = jacobimatrix(P)
+    dat = zeros(T,2,10)
+    for k in 1:10
+        dat[1,k] = (U * (X * (U \ [zeros(k-1); 1; zeros(∞)])))[k]
+        dat[2,k] = (U * (X * (U \ [zeros(k); 1; zeros(∞)])))[k]
+    end
+    return CholeskyJacobiBands{T}(dat, U, X, 10, dat)
+end
+
+size(::CholeskyJacobiBands) = (2,ℵ₀) # Stored as two infinite bands
+
+# Resize and filling functions for cached implementation
+function resizedata!(K::CholeskyJacobiBands, nm::Integer)
+    νμ = K.datasize
+    if nm > νμ
+        olddata = copy(K.data)
+        K.data = similar(K.data, 2, nm)
+        K.data[axes(olddata)...] = olddata
+        inds = νμ:nm
+        cache_filldata!(K, inds)
+        K.datasize = nm
+    end
+    K
+end
+function cache_filldata!(J::CholeskyJacobiBands, inds::UnitRange{Int})
+    for k in inds
+        J.data[1,k] = (J.U * (J.X * (J.U \ [zeros(k-1); 1; zeros(∞)])))[k]
+        J.data[2,k] = (J.U * (J.X * (J.U \ [zeros(k); 1; zeros(∞)])))[k]
+    end
+end
+
+function getindex(K::CholeskyJacobiBands, k::Integer, j::Integer)
+    resizedata!(K, max(k,j))
+    K.data[k, j]
+end
+function getindex(K::CholeskyJacobiBands, kr::Integer, jr::UnitRange{Int})
+    resizedata!(K, maximum(jr))
+    K.data[kr, jr]
+end
+function getindex(K::CholeskyJacobiBands, I::Vararg{Int,2})
+    resizedata!(K,maximum(I))
+    getindex(K.data,I[1],I[2])
+end

test/runtests.jl

@@ -40,6 +40,7 @@ include("test_normalized.jl")
 include("test_lanczos.jl")
 include("test_stieltjes.jl")
 include("test_interlace.jl")
+include("test_decompOPs.jl")
 
 @testset "Auto-diff" begin
     U = Ultraspherical(1)

test/test_decompOPs.jl

@@ -0,0 +1,139 @@
+using Test, ClassicalOrthogonalPolynomials, BandedMatrices, LinearAlgebra, LazyArrays, ContinuumArrays, LazyBandedMatrices
+import ClassicalOrthogonalPolynomials: CholeskyJacobiBands
+import LazyArrays: AbstractCachedMatrix
+import LazyBandedMatrices: SymTridiagonal
+
+@testset "Basic properties" begin
+    # basis
+    P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+    x = axes(P,1)
+    J = jacobimatrix(P)
+    # example weight
+    w(x) = (1 - x^2)
+    W = Symmetric(P \ (w.(x) .* P))
+    # banded cholesky for symmetric-tagged W
+    @test cholesky(W).U isa UpperTriangular
+    # compute Jacobi matrix via cholesky
+    Jchol = cholesky_jacobimatrix(w,P)
+    @test Jchol isa SymTridiagonal
+    # CholeskyJacobiBands object
+    Cbands = CholeskyJacobiBands(W,P)
+    @test Cbands isa CholeskyJacobiBands
+    @test Cbands isa AbstractCachedMatrix
+    @test getindex(Cbands,1,100) == getindex(Cbands,1,1:100)[100]
+end
+
+@testset "Comparison with Lanczos and Classical" begin
+    @testset "Using Clenshaw for polynomial weights" begin
+        @testset "w(x) = x^2*(1-x)" begin
+            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+            x = axes(P,1)
+            J = jacobimatrix(P)
+            wf(x) = x^2*(1-x)
+            # compute Jacobi matrix via cholesky
+            Jchol = cholesky_jacobimatrix(wf, P)
+            # compute Jacobi matrix via classical recurrence
+            Q = Normalized(Jacobi(1,2)[affine(0..1,Inclusion(-1..1)),:])
+            Jclass = jacobimatrix(Q)
+            # 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]
+            # Comparison with Classical
+            @test Jchol[1:500,1:500] ≈ Jclass[1:500,1:500]
+        end
+
+        @testset "w(x) = (1-x^2)" begin
+            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+            x = axes(P,1)
+            J = jacobimatrix(P)
+            wf(x) = (1-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]
+        end
+
+        @testset "w(x) = (1-x^4)" begin
+            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+            x = axes(P,1)
+            J = jacobimatrix(P)
+            wf(x) = (1-x^4)
+            # 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]
+        end
+
+        @testset "w(x) = (1.014-x^3)" begin
+            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+            x = axes(P,1)
+            J = jacobimatrix(P)
+            wf(x) = 1.014-x^4
+            # 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]
+        end
+    end
+
+    @testset "Using Clenshaw with exponential weights" begin
+        @testset "w(x) = exp(x)" begin
+            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+            x = axes(P,1)
+            J = jacobimatrix(P)
+            wf(x) = exp(x)
+            # 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]
+        end
+
+        @testset "w(x) = (1-x)*exp(x)" begin
+            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+            x = axes(P,1)
+            J = jacobimatrix(P)
+            wf(x) = (1-x)*exp(x)
+            # 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]
+        end
+
+        @testset "w(x) = (1-x^2)*exp(x^2)" begin
+            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+            x = axes(P,1)
+            J = jacobimatrix(P)
+            wf(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]
+        end
+
+        @testset "w(x) = x*(1-x^2)*exp(-x^2)" begin
+            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]
+        end
+    end
+end