Associated OPs

Project.toml

@@ -14,6 +14,7 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 InfiniteLinearAlgebra = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
@@ -26,20 +27,21 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 [compat]
 ArrayLayouts = "0.6.2"
 BandedMatrices = "0.16.5"
-BlockArrays = "0.14.1, 0.15"
+BlockArrays = "0.15"
 BlockBandedMatrices = "0.10"
-ContinuumArrays = "0.6.4"
+ContinuumArrays = "0.7"
 DomainSets = "0.4, 0.5"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"
 FastTransforms = "0.11, 0.12"
 FillArrays = "0.11"
+HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.10.3"
 InfiniteLinearAlgebra = "0.5.3"
 IntervalSets = "0.5"
 LazyArrays = "0.21"
-LazyBandedMatrices = "0.5.1"
-QuasiArrays = "0.4.6"
+LazyBandedMatrices = "0.5.5"
+QuasiArrays = "0.5"
 SpecialFunctions = "0.10, 1"
 julia = "1.5"
 

src/ClassicalOrthogonalPolynomials.jl

@@ -2,7 +2,7 @@ module ClassicalOrthogonalPolynomials
 using ContinuumArrays, QuasiArrays, LazyArrays, FillArrays, BandedMatrices, BlockArrays,
     IntervalSets, DomainSets, ArrayLayouts, SpecialFunctions,
     InfiniteLinearAlgebra, InfiniteArrays, LinearAlgebra, FastGaussQuadrature, FastTransforms, FFTW,
-    LazyBandedMatrices
+    LazyBandedMatrices, HypergeometricFunctions
 
 import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
                 IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy,
@@ -126,6 +126,7 @@ Note that `X` is the transpose of the usual definition of the Jacobi matrix.
 """
 jacobimatrix(P) = error("Override for $(typeof(P))")
 
+
 """
     recurrencecoefficients(P)
 
@@ -261,10 +262,52 @@ _vec(a::InfiniteArrays.ReshapedArray) = _vec(parent(a))
 _vec(a::Adjoint{<:Any,<:AbstractVector}) = a'
 
 include("clenshaw.jl")
+include("ratios.jl")
+include("normalized.jl")
+include("lanczos.jl")
+
+function _tritrunc(_, X, n)
+    c,a,b = subdiagonaldata(X),diagonaldata(X),supdiagonaldata(X)
+    Tridiagonal(c[OneTo(n-1)],a[OneTo(n)],b[OneTo(n-1)])
+end
+
+function _tritrunc(::SymTridiagonalLayout, X, n)
+    a,b = diagonaldata(X),supdiagonaldata(X)
+    SymTridiagonal(a[OneTo(n)],b[OneTo(n-1)])
+end
+
+_tritrunc(X, n) = _tritrunc(MemoryLayout(X), X, n)
+
+jacobimatrix(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{Inclusion,OneTo}}) = 
+    _tritrunc(jacobimatrix(parent(V)), maximum(parentindices(V)[2]))
+
+grid(P::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,AbstractUnitRange}}) = 
+    eigvals(symtridiagonalize(jacobimatrix(P)))
 
+function golubwelsch(X)
+    D, V = eigen(symtridiagonalize(X))  # Eigenvalue decomposition
+    D, V[1,:].^2
+end
+
+function golubwelsch(V::SubQuasiArray)
+    x,w = golubwelsch(jacobimatrix(V))
+    w .*= sum(orthogonalityweight(parent(V)))
+    x,w
+end
+
+function factorize(L::SubQuasiArray{T,2,<:Normalized,<:Tuple{Inclusion,OneTo}}) where T
+    x,w = golubwelsch(L)
+    TransformFactorization(x, L[x,:]'*Diagonal(w))
+end
 
-factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{<:Inclusion,<:OneTo}}) where T =
-    TransformFactorization(grid(L), nothing, qr(L[grid(L),:])) # Use QR so type-stable
+
+function factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,OneTo}}) where T
+    x,w = golubwelsch(L)
+    Q = Normalized(parent(L))[parentindices(L)...]
+    D = L \ Q
+    F = factorize(Q)
+    TransformFactorization(F.grid, D*F.plan)
+end
 
 function factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{<:Inclusion,<:AbstractUnitRange}}) where T
     _,jr = parentindices(L)
@@ -290,15 +333,13 @@ function \(wA::WeightedOrthogonalPolynomial, wB::WeightedOrthogonalPolynomial)
     A\B
 end
 
-include("ratios.jl")
-include("normalized.jl")
-include("lanczos.jl")
-include("hermite.jl")
-include("jacobi.jl")
-include("chebyshev.jl")
-include("ultraspherical.jl")
-include("laguerre.jl")
-include("fourier.jl")
+
+include("classical/hermite.jl")
+include("classical/jacobi.jl")
+include("classical/chebyshev.jl")
+include("classical/ultraspherical.jl")
+include("classical/laguerre.jl")
+include("classical/fourier.jl")
 include("stieltjes.jl")
 
 

src/jacobi.jl → src/classical/jacobi.jl

@@ -163,8 +163,10 @@ copy(Q::HalfWeighted) = Q
 
 convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:a,T,<:Jacobi}) where T = JacobiWeight(Q.P.a,zero(T)) .* Q.P
 convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:b,T,<:Jacobi}) where T = JacobiWeight(zero(T),Q.P.b) .* Q.P
-convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:a,T,<:Normalized}) where T = JacobiWeight(Q.P.P.a,zero(T)) .* Q.P
-convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:b,T,<:Normalized}) where T = JacobiWeight(zero(T),Q.P.P.b) .* Q.P
+function convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{lr,T,<:Normalized}) where {T,lr}
+    w,_ = arguments(convert(WeightedOrthogonalPolynomial, HalfWeighted{lr}(Q.P.P)))
+    w .* Q.P
+end
 
 getindex(Q::HalfWeighted, x::Union{Number,AbstractVector}, jr::Union{Number,AbstractVector}) = convert(WeightedOrthogonalPolynomial, Q)[x,jr]
 
@@ -199,7 +201,7 @@ summary(io::IO, P::Jacobi) = print(io, "Jacobi($(P.a), $(P.b))")
 # transforms
 ###
 
-function grid(Pn::SubQuasiArray{T,2,<:AbstractJacobi,<:Tuple{<:Inclusion,<:AbstractUnitRange}}) where T
+function grid(Pn::SubQuasiArray{T,2,<:AbstractJacobi,<:Tuple{Inclusion,AbstractUnitRange}}) where T
     kr,jr = parentindices(Pn)
     ChebyshevGrid{1,T}(maximum(jr))
 end

src/ultraspherical.jl → src/classical/ultraspherical.jl

@@ -10,7 +10,7 @@ struct UltrasphericalWeight{T,Λ} <: AbstractJacobiWeight{T}
 end
 
 UltrasphericalWeight{T}(λ) where T = UltrasphericalWeight{T,typeof(λ)}(λ)
-UltrasphericalWeight(λ) = UltrasphericalWeight{typeof(λ),typeof(λ)}(λ)
+UltrasphericalWeight(λ) = UltrasphericalWeight{float(typeof(λ)),typeof(λ)}(λ)
 
 ==(a::UltrasphericalWeight, b::UltrasphericalWeight) = a.λ == b.λ
 
@@ -19,6 +19,7 @@ function getindex(w::UltrasphericalWeight, x::Number)
     (1-x^2)^(w.λ-one(w.λ)/2)
 end
 
+sum(w::UltrasphericalWeight{T}) where T = sqrt(convert(T,π))*exp(loggamma(one(T)/2 + w.λ)-loggamma(1+w.λ))
 
 
 struct Ultraspherical{T,Λ} <: AbstractJacobi{T}

src/lanczos.jl

@@ -68,6 +68,8 @@ struct LanczosConversion{T} <: LayoutMatrix{T}
 end
 
 size(::LanczosConversion) = (∞,∞)
+copy(R::LanczosConversion) = R
+
 bandwidths(::LanczosConversion) = (0,∞)
 colsupport(L::LanczosConversion, j) = 1:maximum(j)
 rowsupport(L::LanczosConversion, j) = minimum(j):∞

src/stieltjes.jl

@@ -1,78 +1,187 @@
+####
+# Associated
+####
+
+
+"""
+    AssociatedWeighted(P)
+
+We normalise so that `orthogonalityweight(::Associated)` is a probability measure.
+"""
+struct AssociatedWeight{T,OPs<:AbstractQuasiMatrix{T}} <: Weight{T}
+    P::OPs
+end
+axes(w::AssociatedWeight) = (axes(w.P,1),)
+
+sum(::AssociatedWeight{T}) where T = one(T)
+
+"""
+    Associated(P)
+
+constructs the associated orthogonal polynomials for P, which have the Jacobi matrix
+
+    jacobimatrix(P)[2:end,2:end]
+
+and constant first term. Or alternatively
+
+    w = orthogonalityweight(P)
+    A = recurrencecoefficients(P)[1]
+    Associated(P) == (w/(sum(w)*A[1]))'*((P[:,2:end]' - P[:,2:end]) ./ (x' - x))
+
+where `x = axes(P,1)`.
+"""
+
+struct Associated{T, OPs<:AbstractQuasiMatrix{T}} <: OrthogonalPolynomial{T}
+    P::OPs
+end
+
+associated(P) = Associated(P)
+
+axes(Q::Associated) = axes(Q.P)
+==(A::Associated, B::Associated) = A.P == B.P
+
+orthogonalityweight(Q::Associated) = AssociatedWeight(Q.P)
+
+function associated_jacobimatrix(X::Tridiagonal)
+    c,a,b = subdiagonaldata(X),diagonaldata(X),supdiagonaldata(X)
+    Tridiagonal(c[2:end], a[2:end], b[2:end])
+end
+
+function associated_jacobimatrix(X::SymTridiagonal)
+    a,b = diagonaldata(X),supdiagonaldata(X)
+    SymTridiagonal(a[2:end], b[2:end])
+end
+jacobimatrix(a::Associated) = associated_jacobimatrix(jacobimatrix(a.P))
+
+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 ConvKernel{T,D} = BroadcastQuasiMatrix{T,typeof(-),Tuple{D,QuasiAdjoint{T,D}}}
 const Hilbert{T,D} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{ConvKernel{T,Inclusion{T,D}}}}
 const LogKernel{T,D} = BroadcastQuasiMatrix{T,typeof(log),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{T,Inclusion{T,D}}}}}}
 const PowKernel{T,D,F<:Real} = BroadcastQuasiMatrix{T,typeof(^),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{T,Inclusion{T,D}}}},F}}
 
 
-@simplify function *(S::StieltjesPoint{<:Any,<:Any,<:ChebyshevInterval}, wT::WeightedBasis{<:Any,<:ChebyshevTWeight,<:ChebyshevT})
-    w,T = wT.args
-    J = jacobimatrix(T)
-    z, x = parent(S).args[1].args
-    transpose((J'-z*I) \ [-π; zeros(∞)])
+@simplify function *(H::Hilbert, w::ChebyshevTWeight)
+    T = promote_type(eltype(H), eltype(w))
+    zeros(T, axes(w,1))
+end
+
+@simplify function *(H::Hilbert, w::ChebyshevUWeight)
+    T = promote_type(eltype(H), eltype(w))
+    fill(convert(T,π), axes(w,1))
 end
 
-@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval}, wT::WeightedBasis{<:Any,<:ChebyshevTWeight,<:ChebyshevT}) 
+@simplify function *(H::Hilbert, w::LegendreWeight)
+    T = promote_type(eltype(H), eltype(w))
+    x = axes(w,1)
+    log.(x .+ 1) .- log.(1 .- x)
+end
+
+@simplify function *(H::Hilbert, wT::Weighted{<:Any,<:ChebyshevT}) 
     T = promote_type(eltype(H), eltype(wT))
-    ApplyQuasiArray(*, ChebyshevU{T}(), _BandedMatrix(Fill(-convert(T,π),1,∞), ℵ₀, -1, 1))
+    ChebyshevU{T}() * _BandedMatrix(Fill(-convert(T,π),1,∞), ℵ₀, -1, 1)
 end
 
-@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval}, wU::WeightedBasis{<:Any,<:ChebyshevUWeight,<:ChebyshevU}) 
+@simplify function *(H::Hilbert, wU::Weighted{<:Any,<:ChebyshevU}) 
     T = promote_type(eltype(H), eltype(wU))
-    ApplyQuasiArray(*, ChebyshevT{T}(), _BandedMatrix(Fill(convert(T,π),1,∞), ℵ₀, 1, -1))
+    ChebyshevT{T}() * _BandedMatrix(Fill(convert(T,π),1,∞), ℵ₀, 1, -1)
+end
+
+
+@simplify function *(H::Hilbert, wP::Weighted{<:Any,<:OrthogonalPolynomial}) 
+    P = wP.P
+    w = orthogonalityweight(P)
+    A = recurrencecoefficients(P)[1]
+    (-A[1]*sum(w))*[zero(axes(P,1)) associated(P)] + (H*w) .* P
 end
 
+@simplify *(H::Hilbert, P::Legendre) = H * Weighted(P)
+
 ### 
 # LogKernel
 ###
 
-@simplify function *(L::LogKernel{<:Any,<:ChebyshevInterval}, wT::WeightedBasis{<:Any,<:ChebyshevTWeight,<:ChebyshevT}) 
+@simplify function *(L::LogKernel, wT::Weighted{<:Any,<:ChebyshevT}) 
     T = promote_type(eltype(L), eltype(wT))
-    ApplyQuasiArray(*, ChebyshevT{T}(), Diagonal(Vcat(-π*log(2*one(T)),-convert(T,π)./(1:∞))))
+    ChebyshevT{T}() * Diagonal(Vcat(-π*log(2*one(T)),-convert(T,π)./(1:∞)))
 end
 
 
+
 ### 
 # PowKernel
 ###
 
-@simplify function *(K::PowKernel{<:Any,<:ChebyshevInterval}, wT::WeightedBasis{<:Any,<:JacobiWeight,<:Jacobi}) 
+@simplify function *(K::PowKernel, wT::Weighted{<:Any,<:Jacobi}) 
     T = promote_type(eltype(K), eltype(wT))
     cnv,α = K.args
     x,y = K.args[1].args[1].args
     @assert x' == y
     β = (-α-1)/2
-    ApplyQuasiArray(*, ChebyshevT{T}(), Diagonal(1:∞))
+    error("Not implemented")
+    # ChebyshevT{T}() * Diagonal(1:∞)
 end
 
 
 ####
 # StieltjesPoint
 ####
 
-@simplify function *(S::StieltjesPoint, wT::SubQuasiArray{<:Any,2,<:WeightedBasis,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
+stieltjesmoment_jacobi_normalization(n::Int,α::Real,β::Real) = 2^(α+β)*gamma(n+α+1)*gamma(n+β+1)/gamma(2n+α+β+2)
+
+@simplify function *(S::StieltjesPoint, w::AbstractJacobiWeight)
+    α,β = w.a,w.b
+    z,_ = parent(S).args[1].args
+    (x = 2/(1-z);stieltjesmoment_jacobi_normalization(0,α,β)*HypergeometricFunctions.mxa_₂F₁(1,α+1,α+β+2,x))
+end
+
+@simplify function *(S::StieltjesPoint, wP::Weighted)
+    P = wP.P
+    w = orthogonalityweight(P)
+    X = jacobimatrix(P)
+    z, x = parent(S).args[1].args
+    transpose((X'-z*I) \ [-sum(w)*_p0(P); zeros(∞)])
+end
+
+
+@simplify function *(S::StieltjesPoint, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
     P = parent(wT)
     z, x = parent(S).args[1].args
     z̃ = inbounds_getindex(parentindices(wT)[1], z)
     x̃ = axes(P,1)
     (inv.(z̃ .- x̃') * P)[:,parentindices(wT)[2]]
 end
 
-@simplify function *(H::Hilbert, wT::SubQuasiArray{<:Any,2,<:WeightedBasis,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}) 
+@simplify function *(H::Hilbert, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}) 
     P = parent(wT)
     x = axes(P,1)
     apply(*, inv.(x .- x'), P)[parentindices(wT)...]
 end
 
 
-@simplify function *(L::LogKernel, wT::SubQuasiArray{<:Any,2,<:WeightedBasis,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}}) 
+@simplify function *(L::LogKernel, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}}) 
     V = promote_type(eltype(L), eltype(wT))
-    P = parent(wT)
+    wP = parent(wT)
     kr, jr = parentindices(wT)
-    @assert P isa WeightedBasis{<:Any,<:ChebyshevWeight,<:Chebyshev}
-    x = axes(P,1)
-    w,T = P.args
-    D = T \ apply(*, log.(abs.(x .- x')), P)
+    x = axes(wP,1)
+    w,T = arguments(ApplyLayout{typeof(*)}(), wP)
+    @assert w isa ChebyshevTWeight
+    @assert T isa ChebyshevT
+    D = T \ (log.(abs.(x .- x')) * wP)
     c = inv(2*kr.A)
     T[kr,:] * Diagonal(Vcat(2*convert(V,π)*c*log(c), 2c*D.diag.args[2]))
-end
+end
+
+
+
+### generic fallback
+for Op in (:Hilbert, :StieltjesPoint, :LogKernel, :PowKernel)
+    @eval @simplify function *(H::$Op, wP::WeightedBasis{<:Any,<:Weight,<:Any}) 
+        w,P = wP.args
+        Q = OrthogonalPolynomial(w)
+        (H * Weighted(Q)) * (Q \ P)
+    end
+end

test/test_jacobi.jl

@@ -242,7 +242,7 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
         w = JacobiWeight(1.0,1.0)
         wS = w .* S
 
-        W = Diagonal(w)
+        W = QuasiDiagonal(w)
         @test W[0.1,0.2] ≈ 0.0
     end