Fix bug in colsupport for Hilbert (#66)

* Fix bug in colsupport for Hilbert

* add comments

* add comments

* Fix within interval stieltjes evaluation

* Use ChebyshevT

* Update test_stieltjes.jl

* Update test_interlace.jl

* Update Project.toml

Author: dlfivefifty

Project.toml

@@ -28,7 +28,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 ArrayLayouts = "0.7"
 BandedMatrices = "0.16.5"
 BlockArrays = "0.16"
-BlockBandedMatrices = "0.10"
+BlockBandedMatrices = "0.10.9"
 ContinuumArrays = "0.9.1"
 DomainSets = "0.5"
 FFTW = "1.1"
@@ -39,8 +39,8 @@ HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.11.1"
 InfiniteLinearAlgebra = "0.5.8"
 IntervalSets = "0.5"
-LazyArrays = "0.21.15"
-LazyBandedMatrices = "0.6.5"
+LazyArrays = "0.21.16"
+LazyBandedMatrices = "0.6.6"
 QuasiArrays = "0.7"
 SpecialFunctions = "1.0"
 julia = "1.6"

src/stieltjes.jl

@@ -57,7 +57,7 @@ associated(::ChebyshevT{T}) where T = ChebyshevU{T}()
 associated(::ChebyshevU{T}) where T = ChebyshevU{T}()
 
 
-const StieltjesPoint{T,V,D} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{BroadcastQuasiMatrix{T,typeof(-),Tuple{T,QuasiAdjoint{V,Inclusion{V,D}}}}}}
+const StieltjesPoint{T,W<:Number,V,D} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{BroadcastQuasiMatrix{T,typeof(-),Tuple{W,QuasiAdjoint{V,Inclusion{V,D}}}}}}
 const ConvKernel{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(-),Tuple{D1,QuasiAdjoint{T,D2}}}
 const Hilbert{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{ConvKernel{T,Inclusion{T,D1},Inclusion{T,D2}}}}
 const LogKernel{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(log),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{T,Inclusion{T,D1},Inclusion{T,D2}}}}}}
@@ -203,8 +203,8 @@ end
     P = wP.P
     w = orthogonalityweight(P)
     X = jacobimatrix(P)
-    z, x = parent(S).args[1].args
-    z in axes(P,1) && transpose((inv.(x .- x') * wP)[z,:])
+    z, xc = parent(S).args[1].args
+    z in axes(P,1) && return transpose(view(inv.(xc' .- xc) * wP,z,:))
     transpose((X'-z*I) \ [-sum(w)*_p0(P); zeros(∞)])
 end
 
@@ -213,8 +213,8 @@ sqrtx2(x::Real) = sign(x)*sqrt(x^2-1)
 
 @simplify function *(S::StieltjesPoint, wP::Weighted{<:Any,<:ChebyshevU})
     T = promote_type(eltype(S), eltype(wP))
-    z, x = parent(S).args[1].args
-    z in axes(wP,1) && transpose((inv.(x .- x') * wP)[z,:])
+    z, xc = parent(S).args[1].args
+    z in axes(wP,1) && return  (convert(T,π)*ChebyshevT()[z,2:end])'
     ξ = inv(z + sqrtx2(z))
     transpose(convert(T,π) * ξ.^oneto(∞))
 end
@@ -231,7 +231,7 @@ mutable struct HilbertVandermonde{T,MM} <: AbstractCachedMatrix{T}
     colsupport::Vector{Int}
 end
 
-HilbertVandermonde(M, data::Matrix) = HilbertVandermonde(M, data, size(data), Int[])
+HilbertVandermonde(M, data::Matrix) = HilbertVandermonde(M, data, size(data), fill(size(data,1), size(data,2)))
 size(H::HilbertVandermonde) = (ℵ₀,ℵ₀)
 function colsupport(H::HilbertVandermonde, j)
     resizedata!(H, H.datasize[1], maximum(j))
@@ -256,8 +256,9 @@ end
 @simplify function *(H::Hilbert{<:Any,<:Any,<:ChebyshevInterval}, W::Weighted{<:Any,<:ChebyshevU})
     x = axes(H,1)
     T̃ = chebyshevt(x)
-    ψ_1 = T̃ \ inv.(x .+ sqrtx2.(x))
+    ψ_1 = T̃ \ inv.(x .+ sqrtx2.(x)) # same ψ_1 = x .- sqrt(x^2 - 1) but with relative accuracy as x -> ∞
     data = convert(eltype(H),π) * Matrix(reshape(paddeddata(ψ_1),:,1))
+    # Operator has columns π * ψ_1^k 
     T̃ * HilbertVandermonde(Clenshaw(T̃ * ψ_1, T̃), data)
 end
 

test/test_interlace.jl

@@ -1,4 +1,4 @@
-using ClassicalOrthogonalPolynomials, BlockArrays, LazyBandedMatrices, Test
+using ClassicalOrthogonalPolynomials, BlockArrays, LazyBandedMatrices, FillArrays, Test
 import ClassicalOrthogonalPolynomials: PiecewiseInterlace
 
 @testset "Piecewise" begin
@@ -36,12 +36,37 @@ import ClassicalOrthogonalPolynomials: PiecewiseInterlace
     @testset "two-interval Hilbert" begin
         T1,T2 = chebyshevt((-2)..(-1)), chebyshevt(0..2)
         U1,U2 = chebyshevu((-2)..(-1)), chebyshevu(0..2)
-        W = PiecewiseInterlace(Weighted(U1), Weighted(U2)) 
+        W = PiecewiseInterlace(Weighted(U1), Weighted(U2))
         T = PiecewiseInterlace(T1, T2)
+        U = PiecewiseInterlace(U1, U2)
         x = axes(W,1)
         H = T \ inv.(x .- x') * W;
 
+        @test maximum(BlockArrays.blockcolsupport(H,Block(5))) ≤ Block(50)
+
         c = W \ broadcast(x -> exp(x)* (0 ≤ x ≤ 2 ? sqrt(2-x)*sqrt(x) : sqrt(-1-x)*sqrt(x+2)), x)
         @test T[0.5,1:200]'*(H*c)[1:200] ≈ -6.064426633490422
+
+        @testset "inversion" begin
+            H̃ = BlockHcat(Eye((axes(H,1),))[:,Block(1)], H)
+            @test BlockArrays.blockcolsupport(H̃,1) == Block.(1:1)
+            @test BlockArrays.blockcolsupport(H̃,2) == Block.(1:22)
+
+            UT = U \ T
+            D = U \ Derivative(x) * T
+            V = x -> x^4 - 10x^2
+            V_cfs = T \ V.(x)
+            Vp_cfs_U = D * V_cfs
+
+            N = 100
+            Vp_cfs_N = UT[Block.(1:N),Block.(1:N)] \ Vp_cfs_U[Block.(1:N)]
+
+            cμ = H̃[Block.(1:N), Block.(1:N)] \ Vp_cfs_N;
+            c1,c2 = cμ[Block(1)]
+            μ = W[:,Block.(1:N-1)] * cμ[Block.(2:N)]/2;
+            
+            # H * μ == Vp(x) + c1 on first interval
+            # H * μ == Vp(x) + c2 on second interval
+        end
     end
 end