Add convenience functions for symmetric eigen BandedMatrices

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.6.26"
+version = "0.6.27"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"

src/ApproxFunOrthogonalPolynomials.jl

@@ -44,7 +44,7 @@ import ApproxFunBase: Fun, SubSpace, WeightSpace, NoSpace, HeavisideSpace,
                     ℓ⁰, recα, recβ, recγ, ℵ₀, ∞, RectDomain,
                     assert_integer, supportsinplacetransform, ContinuousSpace, SpecialEvalPtType,
                     isleftendpoint, isrightendpoint, evaluation_point, boundaryfn, ldiffbc, rdiffbc,
-                    LeftEndPoint, RightEndPoint
+                    LeftEndPoint, RightEndPoint, promotedomainspace
 
 import DomainSets: Domain, indomain, UnionDomain, FullSpace, Point,
             Interval, ChebyshevInterval, boundary, rightendpoint, leftendpoint,
@@ -139,6 +139,7 @@ include("Spaces/Spaces.jl")
 include("roots.jl")
 include("specialfunctions.jl")
 include("fastops.jl")
+include("symeigen.jl")
 include("show.jl")
 
 end

src/Spaces/Jacobi/Jacobi.jl

@@ -14,13 +14,18 @@ end
 Jacobi(b::T,a::T,d::Domain) where {T<:Number} =
     Jacobi{typeof(d),promote_type(T,real(prectype(d)))}(b, a, d)
 Legendre(domain = ChebyshevInterval()) = Jacobi(0,0,Domain(domain)::Domain)
+Legendre(s::PolynomialSpace) = Legendre(Jacobi(s))
+Legendre(s::Jacobi) = s.a == s.b == 0 ? s : throw(ArgumentError("can't convert $s to Legendre"))
 Jacobi(b::Number,a::Number,d=ChebyshevInterval()) = Jacobi(promote(b, a)..., Domain(d)::Domain)
 Jacobi(A::Ultraspherical) = Jacobi(order(A)-0.5,order(A)-0.5,domain(A))
 Jacobi(A::Chebyshev) = Jacobi(-0.5,-0.5,domain(A))
+Jacobi(A::Jacobi) = A
 
 const NormalizedJacobi{D<:Domain,R,T} = NormalizedPolynomialSpace{Jacobi{D,R,T},D,R}
 NormalizedJacobi(s...) = NormalizedPolynomialSpace(Jacobi(s...))
+NormalizedJacobi(s::Space) = NormalizedPolynomialSpace(Jacobi(_stripnorm(s)))
 NormalizedLegendre(d...) = NormalizedPolynomialSpace(Legendre(d...))
+NormalizedLegendre(s::Space) = NormalizedPolynomialSpace(Legendre(_stripnorm(s)))
 
 function normalization(::Type{T}, sp::Jacobi, k::Int) where T
     x = FastTransforms.Anαβ(k, sp.a, sp.b)

src/symeigen.jl

@@ -0,0 +1,45 @@
+export symmetric_bandmatrices_eigen
+
+function symmetrize_operator_eigen(L, B)
+    L2, B2 = promotedomainspace((L, B))
+    if isambiguous(domainspace(L))
+        throw(ArgumentError("could not detect spaces, please specify the domain spaces for the operators"))
+    end
+    d = domain(L2)
+    if !(domainspace(L2) == Legendre(d) || domainspace(L2) == Ultraspherical(0.5, d))
+        throw(ArgumentError("domainspace of the operators must be $(Legendre(d)) or $(Ultraspherical(0.5,d)) "*
+            "for the symmetric banded conversion, received $(domainspace(L2))"))
+    end
+    S = domainspace(L2)
+    NS = NormalizedPolynomialSpace(S)
+    D1 = Conversion(S, NS)
+    D2 = Conversion(NS, S);
+    QS = QuotientSpace(B2)
+    Q = Conversion(QS, S)
+    R = D1*Q;
+    C = Conversion(domainspace(L2), rangespace(L2))
+    P = cache(PartialInverseOperator(C, (0, bandwidth(L2, 1) + bandwidth(R, 1) + bandwidth(C, 2))));
+    A = R'D1*P*L2*D2*R
+    BB = R'R;
+
+    return A, BB
+end
+
+"""
+    symmetric_bandmatrices_eigen(L::Operator, B::Operator, n::Integer)
+
+Recast the self-adjoint eigenvalue problem `L v = λ v` subject to `B v = 0` to the generalized
+eigenvalue problem `SA v = λ SB v`, where `SA` and `SB` are `Symmetric(::BandedMatrix)`es, and
+return `SA` and `SB`. The `domainspace` of the operators must be one of `Legedre(::Domain)` or
+`Ultraspherical(0.5, ::Domain)`.
+
+!!! note
+    No tests are performed to assert that the system is self-adjoint, and it's the user's responsibility
+    to ensure that this holds.
+"""
+function symmetric_bandmatrices_eigen(L, B, n)
+    AA, BB = symmetrize_operator_eigen(L, B)
+    SA = Symmetric(AA[1:n,1:n], :L)
+    SB = Symmetric(BB[1:n,1:n], :L)
+    return SA, SB
+end

test/EigTest.jl

@@ -221,4 +221,12 @@ using Test
         λ = eigvals(Matrix([B;L][1:n,1:n]), Matrix([B-B;C][1:n,1:n]))
         @test eltype(λ) == Complex{Float64}
     end
+
+    @testset "convenience function for symmetric problems" begin
+        L = -Derivative(Legendre(0..1))^2
+        B = Dirichlet()
+        SA, SB = ApproxFunOrthogonalPolynomials.symmetric_bandmatrices_eigen(L, B, 20)
+        λ = eigvals(SA, SB)[1:4]
+        @test λ ≈ (1:4).^2 .* pi^2 rtol=1e-8
+    end
 end

test/JacobiTest.jl

@@ -183,17 +183,29 @@ using Static
         end
 
         @testset "conversion between spaces" begin
-            for u in (Ultraspherical(1), Chebyshev())
+            for u in (Ultraspherical(1), Chebyshev(), Jacobi(1,1), Legendre())
                 @test NormalizedPolynomialSpace(Jacobi(u)) ==
                     NormalizedJacobi(NormalizedPolynomialSpace(u))
             end
             for j in (Legendre(), Jacobi(1,1))
                 @test NormalizedPolynomialSpace(Ultraspherical(j)) ==
                     NormalizedUltraspherical(NormalizedPolynomialSpace(j))
             end
+            @test NormalizedLegendre(NormalizedLegendre()) == NormalizedLegendre()
+            @test NormalizedJacobi(NormalizedJacobi(1,1)) == NormalizedJacobi(1,1)
+            @test Legendre(Legendre()) == Legendre()
+            @test Legendre(Ultraspherical(0.5)) == Legendre()
+            @test_throws ArgumentError Legendre(Chebyshev())
         end
 
         @test ApproxFunOrthogonalPolynomials.normalization(ComplexF64, Jacobi(-0.5, -0.5), 0) ≈ pi
+
+        @testset "domainspace promotion bug" begin
+            X = Multiplication(Fun(Legendre()), NormalizedLegendre()) → Jacobi(2,2)
+            Y = X : Legendre()
+            @test domainspace(Y) == Legendre()
+            @test Y * Fun(Legendre()) ≈ Fun(x->x^2, Legendre())
+        end
     end
 
     @testset "inplace transform" begin