Marchenko–Pastur example, fix plotting of Q[:,1:n] (#109)

* Marchenko–Pastur example, fix plotting of Q[:,1:4]

* restore plotgrid

* Update normalized.jl

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.6.6"
+version = "0.6.7"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

examples/marchenkopasturops.jl

@@ -0,0 +1,28 @@
+using ClassicalOrthogonalPolynomials, Plots
+
+# MP law
+r = 0.5
+lmin, lmax = (1-sqrt(r))^2,  (1+sqrt(r))^2
+U = chebyshevu(lmin..lmax)
+x = axes(U,1)
+w = @. 1/(2π) * sqrt((lmax-x)*(x-lmin))/(x*r)
+
+# Q is a quasimatrix such that Q[x,k+1] is equivalent to
+# qₖ(x), the k-th orthogonal polynomial wrt to w
+Q = LanczosPolynomial(w, U)
+
+# The Jacobi matrix associated with Q, as an ∞×∞ SymTridiagonal
+J = jacobimatrix(Q)
+
+# plot q₀,…,q₆
+plot(Q[:,1:7])
+
+# Wachter law
+a,b = 5,10
+c,d = sqrt(a/(a+b) * (1-1/(a+b))), sqrt(1/(a+b) * (1-a/(a+b)))
+lmin,lmax = (c-d)^2,(c+d)^2
+U = chebyshevu(lmin..lmax)
+x = axes(U,1)
+w = @. (a+b) * sqrt((x-lmin)*(lmax-x)/(2π*x*(1-x)))
+Q = LanczosPolynomial(w, U)
+plot(Q[:,1:7])

src/ClassicalOrthogonalPolynomials.jl

@@ -245,9 +245,19 @@ _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,Any}}) =
+grid(P::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,OneTo}}) =
     eigvals(symtridiagonalize(jacobimatrix(P)))
 
+function grid(Pn::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,Any}})
+    kr,jr = parentindices(Pn)
+    grid(parent(Pn)[:,oneto(maximum(jr))])
+end
+
+function plotgrid(Pn::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,Any}}) where T
+    kr,jr = parentindices(Pn)
+    grid(parent(Pn)[:,oneto(40maximum(jr))])
+end
+
 function golubwelsch(X)
     D, V = eigen(symtridiagonalize(X))  # Eigenvalue decomposition
     D, V[1,:].^2

src/classical/jacobi.jl

@@ -153,7 +153,7 @@ summary(io::IO, P::Jacobi) = print(io, "Jacobi($(P.a), $(P.b))")
 # transforms
 ###
 
-function grid(Pn::SubQuasiArray{T,2,<:AbstractJacobi,<:Tuple{Inclusion,Any}}) where T
+function grid(Pn::SubQuasiArray{T,2,<:AbstractJacobi,<:Tuple{Inclusion,OneTo}}) where T
     kr,jr = parentindices(Pn)
     ChebyshevGrid{1,T}(maximum(jr))
 end

src/normalized.jl

@@ -83,7 +83,7 @@ QuasiArrays.ApplyQuasiArray(Q::Normalized) = ApplyQuasiArray(*, arguments(ApplyL
 
 ArrayLayouts.mul(Q::Normalized, C::AbstractArray) = ApplyQuasiArray(*, Q, C)
 
-grid(Q::SubQuasiArray{<:Any,2,<:Normalized,<:Tuple{Inclusion,Any}}) = grid(view(parent(Q).P, parentindices(Q)...))
+grid(Q::SubQuasiArray{<:Any,2,<:Normalized,<:Tuple{Inclusion,OneTo}}) = grid(view(parent(Q).P, parentindices(Q)...))
 plotgrid(Q::SubQuasiArray{<:Any,2,<:Normalized,<:Tuple{Inclusion,Any}}) = plotgrid(view(parent(Q).P, parentindices(Q)...))
 
 # transform_ldiv(Q::Normalized, C::AbstractQuasiArray) = Q.scaling .\ (Q.P \ C)

test/test_lanczos.jl

@@ -263,4 +263,36 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedLayout, ort
         R = U \ Q;
         @test R[1:5,1:5] isa Matrix{Float64}
     end
+
+    @testset "Marchenko–Pastur" begin
+        # MP law
+        r = 0.5
+        lmin, lmax = (1-sqrt(r))^2,  (1+sqrt(r))^2
+        U = chebyshevu(lmin..lmax)
+        x = axes(U,1)
+        w = @. 1/(2π) * sqrt((lmax-x)*(x-lmin))/(x*r)
+
+        # Q is a quasimatrix such that Q[x,k+1] is equivalent to
+        # qₖ(x), the k-th orthogonal polynomial wrt to w
+        Q = LanczosPolynomial(w, U)
+
+        # The Jacobi matrix associated with Q, as an ∞×∞ SymTridiagonal
+        J = jacobimatrix(Q)
+
+        @test J[1:3,1:3] ≈ SymTridiagonal([1,1.5,1.5], fill(1/sqrt(2),2))
+
+        # plot q₀,…,q₆
+        @test plotgrid(Q[:,1:7]) ≈ eigvals(Symmetric(J[1:280,1:280]))
+    end
+
+    @testset "Wachter law" begin
+        a,b = 5,10
+        c,d = sqrt(a/(a+b) * (1-1/(a+b))), sqrt(1/(a+b) * (1-a/(a+b)))
+        lmin,lmax = (c-d)^2,(c+d)^2
+        U = chebyshevu(lmin..lmax)
+        x = axes(U,1)
+        w = @. (a+b) * sqrt((x-lmin)*(lmax-x)/(2π*x*(1-x)))
+        Q = LanczosPolynomial(w, U)
+        @test Q[0.5,1:3] ≈ [0.9384176649676137,1.2140849874743076,0.5387898199391473]
+    end
 end