From 068784e4e863d2ae1ae30509a5f4e7574485fd6d Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 21 Jan 2025 22:37:23 +0000 Subject: [PATCH 1/7] Introduce GramMatrixPolynomial --- test/test_gram.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 test/test_gram.jl diff --git a/test/test_gram.jl b/test/test_gram.jl new file mode 100644 index 0000000..82039da --- /dev/null +++ b/test/test_gram.jl @@ -0,0 +1,10 @@ +using ClassicalOrthogonalPolynomials, FastTransforms + +P = Legendre() +x = axes(P,1) +w = @.(1-x^2) + +μ = P'w +X = jacobimatrix(P) +n = 20 +@test GramMatrix(μ[1:2n], X[1:2n,1:2n]) ≈ (P' * (w .* P))[1:n,1:n] \ No newline at end of file From 7da0326f08b39e9370607a22c3d2dffa70055cc4 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 23 Jan 2025 17:02:27 +0000 Subject: [PATCH 2/7] Fixes to ConvertedOP --- Project.toml | 2 +- src/choleskyQR.jl | 19 +++++++++++++++++-- src/normalized.jl | 1 - 3 files changed, 18 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index eeb9112..429ac13 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ClassicalOrthogonalPolynomials" uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd" authors = ["Sheehan Olver "] -version = "0.14.2" +version = "0.14.3" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" diff --git a/src/choleskyQR.jl b/src/choleskyQR.jl index 362e464..15d4f9c 100644 --- a/src/choleskyQR.jl +++ b/src/choleskyQR.jl @@ -8,15 +8,21 @@ struct ConvertedOrthogonalPolynomial{T, WW<:AbstractQuasiVector{T}, XX, UU, PP} P::PP end -_p0(Q::ConvertedOrthogonalPolynomial) = _p0(Q.P) +_p0(Q::ConvertedOrthogonalPolynomial) = _p0(Q.P)/Q.U[1,1] axes(Q::ConvertedOrthogonalPolynomial) = axes(Q.P) + + +struct ConvertedOPLayout <: AbstractOPLayout end MemoryLayout(::Type{<:ConvertedOrthogonalPolynomial}) = ConvertedOPLayout() +equals_layout(::ConvertedOPLayout, ::ConvertedOPLayout, P, Q) = orthogonalityweight(P) == orthogonalityweight(Q) + + jacobimatrix(Q::ConvertedOrthogonalPolynomial) = Q.X orthogonalityweight(Q::ConvertedOrthogonalPolynomial) = Q.weight -# transform to P * U if needed for differentiation, etc. +# transform to P * inv(U) if needed for differentiation, etc. arguments(::ApplyLayout{typeof(*)}, Q::ConvertedOrthogonalPolynomial) = Q.P, ApplyArray(inv, Q.U) OrthogonalPolynomial(w::AbstractQuasiVector) = OrthogonalPolynomial(w, orthogonalpolynomial(singularities(w))) @@ -33,6 +39,12 @@ OrthogonalPolynomial(w::Function, P::AbstractQuasiMatrix) = OrthogonalPolynomial orthogonalpolynomial(w::Function, P::AbstractQuasiMatrix) = orthogonalpolynomial(w.(axes(P,1)), P) +simplifiable(L::Ldiv{ConvertedOPLayout,ConvertedOPLayout}) = simplifiable(Ldiv{ApplyLayout{typeof(*)}, ApplyLayout{typeof(*)}}(L.A, L.B)) +copy(L::Ldiv{ConvertedOPLayout,ConvertedOPLayout}) = copy(Ldiv{ApplyLayout{typeof(*)}, ApplyLayout{typeof(*)}}(L.A, L.B)) +copy(L::Ldiv{ConvertedOPLayout,Lay}) where Lay<:AbstractBasisLayout = copy(Ldiv{ApplyLayout{typeof(*)}, Lay}(L.A, L.B)) +copy(L::Ldiv{ConvertedOPLayout,Lay}) where Lay<:AbstractNormalizedOPLayout = copy(Ldiv{ApplyLayout{typeof(*)}, Lay}(L.A, L.B)) + + """ cholesky_jacobimatrix(w, P) @@ -274,3 +286,6 @@ function _fillqrbanddata!(J::QRJacobiData{:R,T}, inds::UnitRange{Int}) where T ev[k] = UX[k,k-1]/U[k+1,k+1]*(-U[k-1,k+1]/U[k-1,k-1]+U[k-1,k]*U[k,k+1]/(U[k-1,k-1]*U[k,k]))+UX[k,k]/U[k+1,k+1]*(-U[k,k+1]/U[k,k])+UX[k,k+1]/U[k+1,k+1] # this is dot(view(UX,k,k-1:k+1), U[k-1:k+1,k-1:k+1] \ ek) end end + + +diff_layout(::ConvertedOPLayout, P::AbstractQuasiMatrix, dims...) = diff_layout(ApplyLayout{typeof(*)}(), P, dims...) \ No newline at end of file diff --git a/src/normalized.jl b/src/normalized.jl index b0cf762..efc89e8 100644 --- a/src/normalized.jl +++ b/src/normalized.jl @@ -35,7 +35,6 @@ isnormalized(_) = false represents OPs that are of the form P * R where P is another family of OPs and R is upper-triangular. """ abstract type AbstractNormalizedOPLayout <: AbstractOPLayout end -struct ConvertedOPLayout <: AbstractNormalizedOPLayout end struct NormalizedOPLayout{LAY<:AbstractBasisLayout} <: AbstractNormalizedOPLayout end MemoryLayout(::Type{<:Normalized{<:Any, OPs}}) where OPs = NormalizedOPLayout{typeof(MemoryLayout(OPs))}() From fd138ca0aaeadf7d67749c2193b6cf6bcd1311bc Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Fri, 24 Jan 2025 08:02:43 +0000 Subject: [PATCH 3/7] Update choleskyQR.jl --- src/choleskyQR.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/choleskyQR.jl b/src/choleskyQR.jl index 15d4f9c..aef930d 100644 --- a/src/choleskyQR.jl +++ b/src/choleskyQR.jl @@ -15,7 +15,11 @@ axes(Q::ConvertedOrthogonalPolynomial) = axes(Q.P) struct ConvertedOPLayout <: AbstractOPLayout end MemoryLayout(::Type{<:ConvertedOrthogonalPolynomial}) = ConvertedOPLayout() -equals_layout(::ConvertedOPLayout, ::ConvertedOPLayout, P, Q) = orthogonalityweight(P) == orthogonalityweight(Q) +equals_layout(::ConvertedOPLayout, ::ConvertedOPLayout, P, Q) = orthogonalityweight(P) == orthogonalityweight(Q) +equals_layout(::ConvertedOPLayout, ::AbstractOPLayout, P, Q) = false # fix +equals_layout(::AbstractOPLayout, ::ConvertedOPLayout, P, Q) = false # fix + + jacobimatrix(Q::ConvertedOrthogonalPolynomial) = Q.X From 7ae14682473ccd64cac88be76d944014c0c60bb5 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 26 Jan 2025 08:02:28 +0000 Subject: [PATCH 4/7] make ConvertedOP normalized again --- Project.toml | 2 +- src/choleskyQR.jl | 7 ++----- test/test_choleskyQR.jl | 1 - 3 files changed, 3 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 429ac13..643681b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ClassicalOrthogonalPolynomials" uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd" authors = ["Sheehan Olver "] -version = "0.14.3" +version = "0.15" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" diff --git a/src/choleskyQR.jl b/src/choleskyQR.jl index aef930d..fc0e6e5 100644 --- a/src/choleskyQR.jl +++ b/src/choleskyQR.jl @@ -1,5 +1,5 @@ """ -Represent an Orthogonal polynomial which has a conversion operator from P, that is, Q = P * inv(U). +Represents orthonormal polynomials defined via a conversion operator from P, that is, Q = P * inv(U). """ struct ConvertedOrthogonalPolynomial{T, WW<:AbstractQuasiVector{T}, XX, UU, PP} <: OrthonormalPolynomial{T} weight::WW @@ -13,11 +13,8 @@ _p0(Q::ConvertedOrthogonalPolynomial) = _p0(Q.P)/Q.U[1,1] axes(Q::ConvertedOrthogonalPolynomial) = axes(Q.P) -struct ConvertedOPLayout <: AbstractOPLayout end +struct ConvertedOPLayout <: AbstractNormalizedOPLayout end MemoryLayout(::Type{<:ConvertedOrthogonalPolynomial}) = ConvertedOPLayout() -equals_layout(::ConvertedOPLayout, ::ConvertedOPLayout, P, Q) = orthogonalityweight(P) == orthogonalityweight(Q) -equals_layout(::ConvertedOPLayout, ::AbstractOPLayout, P, Q) = false # fix -equals_layout(::AbstractOPLayout, ::ConvertedOPLayout, P, Q) = false # fix diff --git a/test/test_choleskyQR.jl b/test/test_choleskyQR.jl index 839a0a7..73e8cfa 100644 --- a/test/test_choleskyQR.jl +++ b/test/test_choleskyQR.jl @@ -230,7 +230,6 @@ import LazyArrays: AbstractCachedMatrix, resizedata! @test Q[0.1,1] ≈ 1/sqrt(2) @test Q[0.1,1:10] ≈ Q̃[0.1,1:10] - # AWESOME, thanks TSGUT!! @test Q[0.1,10_000] ≈ Q̃[0.1,10_000] R = P \ Q From eb10d6671b220a4f699f058b04b1386988e2170b Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 26 Jan 2025 08:02:57 +0000 Subject: [PATCH 5/7] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 643681b..429ac13 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ClassicalOrthogonalPolynomials" uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd" authors = ["Sheehan Olver "] -version = "0.15" +version = "0.14.3" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" From 3634436623317f4f5a359e0cd9fffc33aea4d234 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 26 Jan 2025 08:04:45 +0000 Subject: [PATCH 6/7] Update choleskyQR.jl --- src/choleskyQR.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/choleskyQR.jl b/src/choleskyQR.jl index fc0e6e5..8f6b527 100644 --- a/src/choleskyQR.jl +++ b/src/choleskyQR.jl @@ -40,12 +40,6 @@ OrthogonalPolynomial(w::Function, P::AbstractQuasiMatrix) = OrthogonalPolynomial orthogonalpolynomial(w::Function, P::AbstractQuasiMatrix) = orthogonalpolynomial(w.(axes(P,1)), P) -simplifiable(L::Ldiv{ConvertedOPLayout,ConvertedOPLayout}) = simplifiable(Ldiv{ApplyLayout{typeof(*)}, ApplyLayout{typeof(*)}}(L.A, L.B)) -copy(L::Ldiv{ConvertedOPLayout,ConvertedOPLayout}) = copy(Ldiv{ApplyLayout{typeof(*)}, ApplyLayout{typeof(*)}}(L.A, L.B)) -copy(L::Ldiv{ConvertedOPLayout,Lay}) where Lay<:AbstractBasisLayout = copy(Ldiv{ApplyLayout{typeof(*)}, Lay}(L.A, L.B)) -copy(L::Ldiv{ConvertedOPLayout,Lay}) where Lay<:AbstractNormalizedOPLayout = copy(Ldiv{ApplyLayout{typeof(*)}, Lay}(L.A, L.B)) - - """ cholesky_jacobimatrix(w, P) @@ -287,6 +281,3 @@ function _fillqrbanddata!(J::QRJacobiData{:R,T}, inds::UnitRange{Int}) where T ev[k] = UX[k,k-1]/U[k+1,k+1]*(-U[k-1,k+1]/U[k-1,k-1]+U[k-1,k]*U[k,k+1]/(U[k-1,k-1]*U[k,k]))+UX[k,k]/U[k+1,k+1]*(-U[k,k+1]/U[k,k])+UX[k,k+1]/U[k+1,k+1] # this is dot(view(UX,k,k-1:k+1), U[k-1:k+1,k-1:k+1] \ ek) end end - - -diff_layout(::ConvertedOPLayout, P::AbstractQuasiMatrix, dims...) = diff_layout(ApplyLayout{typeof(*)}(), P, dims...) \ No newline at end of file From 7cda647c7d1a627f1238d6e90262c5371be0419d Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 26 Jan 2025 12:28:45 +0000 Subject: [PATCH 7/7] Update test_choleskyQR.jl --- test/test_choleskyQR.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/test/test_choleskyQR.jl b/test/test_choleskyQR.jl index 73e8cfa..905c82a 100644 --- a/test/test_choleskyQR.jl +++ b/test/test_choleskyQR.jl @@ -1,5 +1,5 @@ using Test, ClassicalOrthogonalPolynomials, BandedMatrices, LinearAlgebra, LazyArrays, ContinuumArrays, LazyBandedMatrices, InfiniteLinearAlgebra -import ClassicalOrthogonalPolynomials: cholesky_jacobimatrix, qr_jacobimatrix, orthogonalpolynomial +import ClassicalOrthogonalPolynomials: cholesky_jacobimatrix, qr_jacobimatrix, orthogonalpolynomial, _p0 import LazyArrays: AbstractCachedMatrix, resizedata! @testset "CholeskyQR" begin @@ -228,7 +228,7 @@ import LazyArrays: AbstractCachedMatrix, resizedata! @test Q == Q̃ @test Q̃ == Q - @test Q[0.1,1] ≈ 1/sqrt(2) + @test Q[0.1,1] ≈ _p0(Q) ≈ 1/sqrt(2) @test Q[0.1,1:10] ≈ Q̃[0.1,1:10] @test Q[0.1,10_000] ≈ Q̃[0.1,10_000] @@ -246,6 +246,9 @@ import LazyArrays: AbstractCachedMatrix, resizedata! @testset "Chebyshev" begin U = ChebyshevU() Q = orthogonalpolynomial(x -> (1+x^2)*sqrt(1-x^2), U) + x = axes(U,1) + + @test Q[0.1,1] ≈ _p0(Q) ≈ 1/sqrt(sum(expand(Weighted(U),x -> (1+x^2)*sqrt(1-x^2)))) @test bandwidths(Q\U) == (0,2) Q̃ = OrthogonalPolynomial(x -> (1+x^2)*sqrt(1-x^2), U)