Skip to content

Commit e75aacc

Browse files
TSGutTSGut
committed
start implementing Cholesky-Jacobi approach
1 parent 94d5843 commit e75aacc

File tree

4 files changed

+167
-0
lines changed

4 files changed

+167
-0
lines changed

‎src/ClassicalOrthogonalPolynomials.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -361,6 +361,9 @@ include("classical/ultraspherical.jl")
361361
include("classical/laguerre.jl")
362362
include("classical/fourier.jl")
363363
include("stieltjes.jl")
364+
include("decompOPs.jl")
365+
366+
export cholesky_jacobimatrix
364367

365368

366369
end # module

‎src/decompOPs.jl

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
# This currently takes the weight multiplication operator as input.
2+
# I will probably change this to take the weight function instead.
3+
function cholesky_jacobimatrix(W)
4+
bands = CholeskyJacobiBands(W) # the cached array only needs to store two bands bc of symmetry
5+
return SymTridiagonal(bands[1,:],bands[2,:])
6+
end
7+
8+
# The generated Jacobi operators are symmetric tridiagonal, so we store their data as two bands
9+
mutable struct CholeskyJacobiBands{T} <: AbstractCachedMatrix{T}
10+
data
11+
U
12+
X
13+
datasize::Int
14+
array
15+
end
16+
17+
# SymTridiagonal currently doesn't parse as Symmetric, so here's a Q&D workaround for conversion
18+
symmjacobim(J::SymTridiagonal) = Symmetric(BandedMatrix(0=>J.dv, 1=>J.ev))
19+
20+
# Computes the initial data for the Jacobi operator bands
21+
function CholeskyJacobiBands(W)
22+
T = eltype(W)
23+
U = cholesky(W).U
24+
X = symmjacobim(jacobimatrix(Normalized(Legendre()[affine(zero(T)..one(T),Inclusion(-one(T)..one(T))),:])))
25+
dat = zeros(T,2,10)
26+
for k in 1:10
27+
dat[1,k] = (U * (X * (U \ [zeros(k-1); 1; zeros(∞)])))[k]
28+
dat[2,k] = (U * (X * (U \ [zeros(k); 1; zeros(∞)])))[k]
29+
end
30+
return CholeskyJacobiBands{T}(dat, U, X, 10, dat)
31+
end
32+
33+
size(::CholeskyJacobiBands) = (2,ℵ₀) # Stored as two infinite bands
34+
35+
# Resize and filling functions for cached implementation
36+
function resizedata!(K::CholeskyJacobiBands, nm::Integer)
37+
νμ = K.datasize
38+
if nm > νμ
39+
olddata = copy(K.data)
40+
K.data = similar(K.data, nm, nm)
41+
K.data[axes(olddata)...] = olddata
42+
inds = νμ:nm
43+
cache_filldata!(K, inds)
44+
K.datasize = nm
45+
end
46+
K
47+
end
48+
function cache_filldata!(J::CholeskyJacobiBands, inds)
49+
for k in inds
50+
J.data[1,k] = (J.U * (J.X * (J.U \ [zeros(k-1); 1; zeros(∞)])))[k]
51+
J.data[2,k] = (J.U * (J.X * (J.U \ [zeros(k); 1; zeros(∞)])))[k]
52+
end
53+
end
54+
55+
# Need to check which of these are actually useful in practice with tests
56+
function getindex(K::CholeskyJacobiBands, k::Integer, j::Integer)
57+
resizedata!(K, max(k,j))
58+
K.data[k, j]
59+
end
60+
function getindex(K::CholeskyJacobiBands, kr::Integer, jr::UnitRange{Integer})
61+
resizedata!(K, maximum(jr))
62+
K.data[kr, jr]
63+
end
64+
function getindex(K::CholeskyJacobiBands, I::Vararg{Int,2})
65+
resizedata!(K,maximum(I))
66+
getindex(K.data,I[1],I[2])
67+
end

‎test/runtests.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ include("test_normalized.jl")
4040
include("test_lanczos.jl")
4141
include("test_stieltjes.jl")
4242
include("test_interlace.jl")
43+
include("test_decompOPs.jl")
4344

4445
@testset "Auto-diff" begin
4546
U = Ultraspherical(1)

‎test/test_decompOPs.jl

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
using Test, ClassicalOrthogonalPolynomials, BandedMatrices, LinearAlgebra, LazyArrays, ContinuumArrays, LazyBandedMatrices
2+
import ClassicalOrthogonalPolynomials: symmjacobim
3+
4+
@testset begin "Test the Q&D conversion to BandedMatrix format"
5+
# Legendre
6+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
7+
x = axes(P,1)
8+
J = jacobimatrix(P)
9+
Jx = symmjacobim(J)
10+
@test J[1:100,1:100] == Jx[1:100,1:100]
11+
# Jacobi
12+
P = Normalized(Jacobi(1,2)[affine(0..1,Inclusion(-1..1)),:])
13+
x = axes(P,1)
14+
J = jacobimatrix(P)
15+
Jx = symmjacobim(J)
16+
@test J[1:100,1:100] == Jx[1:100,1:100]
17+
end
18+
19+
@testset begin "Basic properties of Cholesky based Jacobi operators"
20+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
21+
x = axes(P,1)
22+
J = jacobimatrix(P)
23+
Jx = symmjacobim(J)
24+
w = (I - Jx^2)
25+
# banded cholesky for symmetric-tagged W
26+
@test cholesky(w).U isa UpperTriangular
27+
# compute Jacobi matrix via cholesky
28+
Jchol = cholesky_jacobimatrix(w)
29+
@test Jchol isa LazyBandedMatrices.SymTridiagonal
30+
end
31+
32+
@testset begin "Comparison with Lanczos and Classical for w(x) = x^2*(1-x)"
33+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
34+
x = axes(P,1)
35+
J = jacobimatrix(P)
36+
Jx = symmjacobim(J)
37+
w = (Jx^2 - Jx^3)
38+
# compute Jacobi matrix via cholesky
39+
Jchol = cholesky_jacobimatrix(w)
40+
# compute Jacobi matrix via classical recurrence
41+
Q = Normalized(Jacobi(1,2)[affine(0..1,Inclusion(-1..1)),:])
42+
Jclass = jacobimatrix(Q)
43+
# compute Jacobi matrix via Lanczos
44+
wf = x.^2 .* (1 .- x)
45+
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
46+
# Comparison with Lanczos
47+
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
48+
# Comparison with Classical
49+
@test Jchol[1:500,1:500] ≈ Jclass[1:500,1:500]
50+
end
51+
52+
@testset begin "Comparison with Lanczos for w(x) = (1-x^2)"
53+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
54+
x = axes(P,1)
55+
J = jacobimatrix(P)
56+
Jx = symmjacobim(J)
57+
w = (I - Jx^2)
58+
# compute Jacobi matrix via cholesky
59+
Jchol = cholesky_jacobimatrix(w)
60+
# compute Jacobi matrix via Lanczos
61+
wf = (1 .- x.^2)
62+
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
63+
# Comparison
64+
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
65+
end
66+
67+
@testset begin "Comparison with Lanczos for w(x) = (1-x^4)"
68+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
69+
x = axes(P,1)
70+
J = jacobimatrix(P)
71+
Jx = symmjacobim(J)
72+
w = (I - Jx^4)
73+
# compute Jacobi matrix via cholesky
74+
Jchol = cholesky_jacobimatrix(w)
75+
# compute Jacobi matrix via Lanczos
76+
wf = (1 .- x.^4)
77+
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
78+
# Comparison
79+
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
80+
end
81+
82+
@testset begin "Comparison with Lanczos for w(x) = (1.014-x^3)"
83+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
84+
x = axes(P,1)
85+
J = jacobimatrix(P)
86+
Jx = symmjacobim(J)
87+
t = 1.014
88+
w = (t*I - Jx^3)
89+
# compute Jacobi matrix via cholesky
90+
Jchol = cholesky_jacobimatrix(w)
91+
# compute Jacobi matrix via Lanczos
92+
wf = (t .- x.^3)
93+
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
94+
# Comparison
95+
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
96+
end

0 commit comments

Comments
 (0)