Add convenience functions for symmetric eigen BandedMatrices (#240)

jishnub · web-flow · commit fcfe37ce5610 · 2023-04-28T11:29:32.000+05:30
* Add convenience functions for symmetric eigen BandedMatrices

* update docstring to mention n

* update docstring

* Add Ultraspherical test

* Add SymmetricEigensystem

* eigvals for SymmetricEigensystem

* Don't restrict domainspace

* access domains for PiecewiseSpace using components(d)

* Add SkewSymmetricEigensystem

* update docstrings

* move docstring

* Specify QuotientSpace type in SymmetricEigensystem constructor

* remove symmetric_bandmatrices_eigen

* collect matrices in skew-symmetric system
diff --git a/Project.toml b/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"
@@ -20,7 +20,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-ApproxFunBase = "0.8.21"
+ApproxFunBase = "0.8.22"
 ApproxFunBaseTest = "0.1"
 Aqua = "0.5"
 BandedMatrices = "0.16, 0.17"
diff --git a/src/ApproxFunOrthogonalPolynomials.jl b/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, normalizedspace
+                    LeftEndPoint, RightEndPoint, normalizedspace, promotedomainspace
 
 import DomainSets: Domain, indomain, UnionDomain, FullSpace, Point,
             Interval, ChebyshevInterval, boundary, rightendpoint, leftendpoint,
@@ -82,7 +82,7 @@ import SpecialFunctions: erfcx, dawson,
 
 using StaticArrays: SVector
 
-import LinearAlgebra: isdiag
+import LinearAlgebra: isdiag, eigvals
 
 points(d::IntervalOrSegmentDomain{T},n::Integer) where {T} =
     _maybefromcanonical(d, chebyshevpoints(float(real(eltype(T))), n))  # eltype to handle point
@@ -140,6 +140,7 @@ include("Spaces/Spaces.jl")
 include("roots.jl")
 include("specialfunctions.jl")
 include("fastops.jl")
+include("symeigen.jl")
 include("show.jl")
 
 end
diff --git a/src/Spaces/Jacobi/Jacobi.jl b/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)
diff --git a/src/symeigen.jl b/src/symeigen.jl
@@ -0,0 +1,118 @@
+export bandmatrices_eigen, SymmetricEigensystem, SkewSymmetricEigensystem
+
+abstract type EigenSystem end
+
+"""
+    SymmetricEigensystem(L::Operator, B::Operator, QuotientSpaceType::Type = QuotientSpace)
+
+Represent the eigensystem `L v = λ v` subject to `B v = 0`, where `L` is self-adjoint with respect to the standard `L2`
+inner product given the boundary conditions `B`.
+
+!!! note
+    No tests are performed to assert that the operator `L` is self-adjoint, and it's the user's responsibility
+    to ensure that the operators are compliant.
+
+The optional argument `QuotientSpaceType` specifies the type of space to be used to denote the quotient space in the basis
+recombination process. In most cases, the default choice of `QuotientSpace` is a good one. In specific instances where `B`
+is rank-deficient (e.g. it contains a column of zeros,
+which typically happens if one of the basis elements already satiafies the boundary conditions),
+one may need to choose this to be a `PathologicalQuotientSpace`.
+
+!!! note
+    No checks on the rank of `B` are carried out, and it's up to the user to specify the correct type.
+"""
+SymmetricEigensystem
+
+"""
+    SkewSymmetricEigensystem(L::Operator, B::Operator, QuotientSpaceType::Type = QuotientSpace)
+
+Represent the eigensystem `L v = λ v` subject to `B v = 0`, where `L` is skew-symmetric with respect to the standard `L2`
+inner product given the boundary conditions `B`.
+
+!!! note
+    No tests are performed to assert that the operator `L` is skew-symmetric, and it's the user's responsibility
+    to ensure that the operators are compliant.
+
+The optional argument `QuotientSpaceType` specifies the type of space to be used to denote the quotient space in the basis
+recombination process. In most cases, the default choice of `QuotientSpace` is a good one. In specific instances where `B`
+is rank-deficient (e.g. it contains a column of zeros,
+which typically happens if one of the basis elements already satiafies the boundary conditions),
+one may need to choose this to be a `PathologicalQuotientSpace`.
+
+!!! note
+    No checks on the rank of `B` are carried out, and it's up to the user to specify the correct type.
+"""
+SkewSymmetricEigensystem
+
+for SET in (:SymmetricEigensystem, :SkewSymmetricEigensystem)
+    @eval begin
+        struct $SET{LT,QT} <: EigenSystem
+            L :: LT
+            Q :: QT
+
+            function $SET(L, B, ::Type{QST} = QuotientSpace) where {QST}
+                L2, B2 = promotedomainspace((L, B))
+                if isambiguous(domainspace(L))
+                    throw(ArgumentError("could not detect spaces, please specify the domain spaces for the operators"))
+                end
+
+                QS = QST(B2)
+                S = domainspace(L2)
+                Q = Conversion(QS, S)
+                new{typeof(L2),typeof(Q)}(L2, Q)
+            end
+        end
+    end
+end
+
+function basis_recombination(SE::EigenSystem)
+    L, Q = SE.L, SE.Q
+    S = domainspace(L)
+    D1 = Conversion_normalizedspace(S)
+    D2 = Conversion_normalizedspace(S, Val(:backward))
+    R = D1*Q;
+    C = Conversion(S, rangespace(L))
+    P = cache(PartialInverseOperator(C, (0, bandwidth(L, 1) + bandwidth(R, 1) + bandwidth(C, 2))));
+    A = R'D1*P*L*D2*R
+    BB = R'R;
+
+    return A, BB
+end
+
+"""
+    bandmatrices_eigen(S::Union{SymmetricEigensystem, SkewSymmetricEigensystem}, n::Integer)
+
+Recast the symmetric/skew-symmetric eigenvalue problem `L v = λ v` subject to `B v = 0` to the generalized
+eigenvalue problem `SA v = λ SB v`, where `SA` and `SB` are banded operators, and
+return the `n × n` matrix representations of `SA` and `SB`.
+If `S isa SymmetricEigensystem`, the returned matrices will be `Symmetric`.
+
+!!! note
+    No tests are performed to assert that the system is symmetric/skew-symmetric, and it's the user's responsibility
+    to ensure that the operators are compliant.
+"""
+function bandmatrices_eigen(S::SymmetricEigensystem, n::Integer)
+    A, B = _bandmatrices_eigen(S, n)
+    SA = Symmetric(A, :L)
+    SB = Symmetric(B, :L)
+    return SA, SB
+end
+
+function bandmatrices_eigen(S::EigenSystem, n::Integer)
+    A, B = _bandmatrices_eigen(S, n)
+    A2 = tril(A, bandwidth(A,1))
+    B2 = tril(B, bandwidth(B,1))
+    Matrix(A2), Matrix(B2)
+end
+
+function _bandmatrices_eigen(S::EigenSystem, n::Integer)
+    AA, BB = basis_recombination(S)
+    A = AA[1:n,1:n]
+    B = BB[1:n,1:n]
+    return A, B
+end
+
+function eigvals(S::EigenSystem, n::Integer)
+    SA, SB = bandmatrices_eigen(S, n)
+    eigvals(SA, SB)
+end
diff --git a/test/EigTest.jl b/test/EigTest.jl
@@ -2,6 +2,7 @@ using ApproxFunOrthogonalPolynomials
 using ApproxFunBase
 using ApproxFunBase: bandwidth
 using BandedMatrices
+using LinearAlgebra
 using SpecialFunctions
 using Test
 
@@ -11,28 +12,24 @@ using Test
         # -𝒟² u = λu,  u'(±1) = 0.
         #
         d = Segment(-1..1)
-        S = Ultraspherical(0.5, d)
-        NS = NormalizedPolynomialSpace(S)
-        L = -Derivative(S, 2)
-        C = Conversion(domainspace(L), rangespace(L))
-        B = Neumann(S)
-        QS = QuotientSpace(B)
-        Q = Conversion(QS, S)
-        D1 = Conversion(S, NS)
-        D2 = Conversion(NS, S)
-        R = D1*Q
-        P = cache(PartialInverseOperator(C, (0, bandwidth(L, 1) + bandwidth(R, 1) + bandwidth(C, 2))))
-        A = R'D1*P*L*D2*R
-        B = R'R
-
-        # Currently a hack to avoid a LAPACK calling bug with a pencil (A, B)
-        # with smaller bandwidth in A.
-        n = 50
-        SA = Symmetric(Matrix(A[1:n,1:n]), :L)
-        SB = Symmetric(Matrix(B[1:n,1:n]), :L)
-        λ = eigvals(SA, SB)
-
-        @test λ[1:round(Int, 2n/5)] ≈ (π^2/4).*(0:round(Int, 2n/5)-1).^2
+        for S in (Ultraspherical(0.5, d), Legendre(d))
+            L = -Derivative(S, 2)
+            B = Neumann(S)
+
+            n = 50
+            Seig = SymmetricEigensystem(L, B)
+            SA, SB = bandmatrices_eigen(Seig, n)
+
+            # A hack to avoid a BandedMatrices bug on Julia v1.6, with a pencil (A, B)
+            # with smaller bandwidth in A.
+            λ = if VERSION >= v"1.8"
+                eigvals(SA, SB)
+            else
+                eigvals(Symmetric(Matrix(SA)), Symmetric(Matrix(SB)))
+            end
+
+            @test λ[1:round(Int, 2n/5)] ≈ (π^2/4).*(0:round(Int, 2n/5)-1).^2
+        end
     end
 
     @testset "Schrödinger with piecewise-linear potential with Dirichlet boundary conditions" begin
@@ -42,25 +39,14 @@ using Test
         # where V = 100|x|.
         #
         d = Segment(-1..0)∪Segment(0..1)
-        S = PiecewiseSpace(Ultraspherical.(0.5, d.domains))
-        NS = PiecewiseSpace(NormalizedUltraspherical.(0.5, d.domains))
+        S = PiecewiseSpace(Ultraspherical.(0.5, components(d)))
         V = 100Fun(abs, S)
         L = -Derivative(S, 2) + V
-        C = Conversion(domainspace(L), rangespace(L))
         B = [Dirichlet(S); continuity(S, 0:1)]
-        QS = QuotientSpace(B)
-        Q = Conversion(QS, S)
-        D1 = Conversion(S, NS)
-        D2 = Conversion(NS, S)
-        R = D1*Q
-        P = cache(PartialInverseOperator(C, (0, bandwidth(L, 1) + bandwidth(R, 1) + bandwidth(C, 2))))
-        A = R'D1*P*L*D2*R
-        B = R'R
 
+        Seig = SymmetricEigensystem(L, B)
         n = 100
-        SA = Symmetric(A[1:n,1:n], :L)
-        SB = Symmetric(B[1:n,1:n], :L)
-        λ = eigvals(SA, SB)
+        λ = eigvals(Seig, n)
 
         @test λ[1] ≈ parse(BigFloat, "2.19503852085715200848808942880214615154684642693583513254593767079468401198338e+01")
     end
@@ -71,26 +57,17 @@ using Test
         #
         # where V = 1000[χ_[-1,-1/2](x) + χ_[1/2,1](x)].
         #
-        d = Segment(-1..(-0.5))∪Segment(-0.5..0.5)∪Segment(0.5..1)
-        S = PiecewiseSpace(Ultraspherical.(0.5, d.domains))
-        NS = PiecewiseSpace(NormalizedUltraspherical.(0.5, d.domains))
+        d = Segment(-1..(-0.5)) ∪ Segment(-0.5..0.5) ∪ Segment(0.5..1)
+        S = PiecewiseSpace(Ultraspherical.(0.5, components(d)))
+
         V = Fun(x->abs(x) ≥ 1/2 ? 1000 : 0, S)
         L = -Derivative(S, 2) + V
-        C = Conversion(domainspace(L), rangespace(L))
         B = [Dirichlet(S); continuity(S, 0:1)]
-        QS = QuotientSpace(B)
-        Q = Conversion(QS, S)
-        D1 = Conversion(S, NS)
-        D2 = Conversion(NS, S)
-        R = D1*Q
-        P = cache(PartialInverseOperator(C, (0, bandwidth(L, 1) + bandwidth(R, 1) + bandwidth(C, 2))))
-        A = R'D1*P*L*D2*R
-        B = R'R
+
+        Seig = SymmetricEigensystem(L, B)
 
         n = 150
-        SA = Symmetric(A[1:n,1:n], :L)
-        SB = Symmetric(B[1:n,1:n], :L)
-        λ = eigvals(SA, SB)
+        λ = eigvals(Seig, n)
         # From Lee--Greengard (1997).
         λtrue = [2.95446;5.90736;8.85702;11.80147].^2
         @test norm((λ[1:4] - λtrue)./λ[1:4]) < 1e-5
@@ -102,36 +79,27 @@ using Test
         #
         # where V = x + 100δ(x-0.25).
         #
-        d = Segment(-1..0.25)∪Segment(0.25..1)
-        S = PiecewiseSpace(Ultraspherical.(0.5, d.domains))
-        NS = PiecewiseSpace(NormalizedUltraspherical.(0.5, d.domains))
+        d = Segment(-1..0.25) ∪ Segment(0.25..1)
+        S = PiecewiseSpace(Ultraspherical.(0.5, components(d)))
         V = Fun(identity, S)
         L = -Derivative(S, 2) + V
-        C = Conversion(domainspace(L), rangespace(L))
+
         B4 = zeros(Operator{ApproxFunBase.prectype(S)}, 1, 2)
         B4[1, 1] = -Evaluation(component(S, 1), rightendpoint, 1) - 100*0.5*Evaluation(component(S, 1), rightendpoint)
         B4[1, 2] = Evaluation(component(S, 2), leftendpoint, 1)  - 100*0.5*Evaluation(component(S, 2), leftendpoint)
         B4 = ApproxFunBase.InterlaceOperator(B4, PiecewiseSpace, ApproxFunBase.ArraySpace)
         B = [Evaluation(S, -1); Evaluation(S, 1) + Evaluation(S, 1, 1); continuity(S, 0); B4]
-        QS = QuotientSpace(B)
-        Q = Conversion(QS, S)
-        D1 = Conversion(S, NS)
-        D2 = Conversion(NS, S)
-        R = D1*Q
-        P = cache(PartialInverseOperator(C, (0, bandwidth(L, 1) + bandwidth(R, 1) + bandwidth(C, 2))))
-        A = R'D1*P*L*D2*R
-        B = R'R
 
         n = 100
-        SA = Symmetric(A[1:n,1:n], :L)
-        SB = Symmetric(B[1:n,1:n], :L)
+        Seig = SymmetricEigensystem(L, B)
+        SA, SB = bandmatrices_eigen(Seig, n)
+        λ, Q = eigen(SA, SB);
 
+        QS = QuotientSpace(B)
         k = 3
-
-        λ, Q = eigen(SA, SB);
         u_QS = Fun(QS, Q[:, k])
         u_S = Fun(u_QS, S)
-        u = Fun(u_S, PiecewiseSpace(Chebyshev.(d.domains)))
+        u = Fun(u_S, PiecewiseSpace(Chebyshev.(components(d))))
         u /= sign(u'(-1))
         u1, u2 = components(u)
 
@@ -149,22 +117,13 @@ using Test
         #
         d = Segment(big(-1.0)..big(1.0))
         S = Ultraspherical(big(0.5), d)
-        NS = NormalizedPolynomialSpace(S)
         L = -Derivative(S, 2)
         C = Conversion(domainspace(L), rangespace(L))
         B = Dirichlet(S)
-        QS = QuotientSpace(B)
-        Q = Conversion(QS, S)
-        D1 = Conversion(S, NS)
-        D2 = Conversion(NS, S)
-        R = D1*Q
-        P = cache(PartialInverseOperator(C, (0, bandwidth(L, 1) + bandwidth(R, 1) + bandwidth(C, 2))))
-        A = R'D1*P*L*D2*R
-        B = R'R
 
         n = 300
-        SA = Symmetric(A[1:n,1:n], :L)
-        SB = Symmetric(B[1:n,1:n], :L)
+        Seig = SymmetricEigensystem(L, B)
+        SA, SB = bandmatrices_eigen(Seig, n)
         BSA = BandedMatrix(SA)
         BSB = BandedMatrix(SB)
         begin
@@ -189,21 +148,13 @@ using Test
         #
         d = Segment(-1..1)
         S = Ultraspherical(0.5, d)
-        NS = NormalizedPolynomialSpace(S)
         Lsk = Derivative(S)
-        Csk = Conversion(domainspace(Lsk), rangespace(Lsk))
         B = Evaluation(S, -1) + Evaluation(S, 1)
-        QS = PathologicalQuotientSpace(B)
-        Q = Conversion(QS, S)
-        D1 = Conversion(S, NS)
-        D2 = Conversion(NS, S)
-        R = D1*Q
-        Psk = cache(PartialInverseOperator(Csk, (0, 4 + bandwidth(Lsk, 1) + bandwidth(R, 1) + bandwidth(Csk, 2))))
-        Ask = R'D1*Psk*Lsk*D2*R
-        Bsk = R'R
+        Seig = SkewSymmetricEigensystem(Lsk, B, PathologicalQuotientSpace)
 
         n = 100
-        λim = imag(sort!(eigvals(Matrix(tril(Ask[1:n,1:n], 1)), Matrix(tril(Bsk[1:n,1:n], 2))), by = abs))
+        λ = eigvals(Seig, n)
+        λim = imag(sort!(λ, by = abs))
 
         @test abs.(λim[1:2:round(Int, 2n/5)]) ≈ π.*(0.5:round(Int, 2n/5)/2)
         @test abs.(λim[2:2:round(Int, 2n/5)]) ≈ π.*(0.5:round(Int, 2n/5)/2)
diff --git a/test/JacobiTest.jl b/test/JacobiTest.jl
