Add tridiagonal Reverse Cholesky (#167)

* Generalise QL

* tridiagonal reverse cholesky

* Update test_infql.jl

* v0.7.5

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.7.4"
+version = "0.7.5"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -16,17 +16,22 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 
 [compat]
-Aqua = "0.6"
-ArrayLayouts = "1.0.12"
+Aqua = "0.8"
+ArrayLayouts = "1.9.2"
 BandedMatrices = "0.17.19, 1"
 BlockArrays = "0.16.14"
 BlockBandedMatrices = "0.12"
 FillArrays = "1"
 InfiniteArrays = "0.13"
 LazyArrays = "1.3"
-LazyBandedMatrices = "0.8.7, 0.9"
-MatrixFactorizations = "2.1"
+LazyBandedMatrices = "0.9"
+LinearAlgebra = "1"
+MatrixFactorizations = "2.2"
+Random = "1"
 SemiseparableMatrices = "0.3"
+SpecialFunctions = "2"
+StaticArrays = "1"
+Test = "1"
 julia = "1.9"
 
 [extras]

src/InfiniteLinearAlgebra.jl

@@ -14,15 +14,15 @@ import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, trian
                         _bidiag_forwardsub!, mulreduce, RangeCumsum, _factorize, transposelayout, ldiv!, lmul!, mul, CNoPivot
 import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns, BandedLayout,
                         _default_banded_broadcast, banded_similar
-import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row
+import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row, AbstractFillMatrix, AbstractFillVector
 import InfiniteArrays: OneToInf, InfUnitRange, Infinity, PosInfinity, InfiniteCardinal, InfStepRange, AbstractInfUnitRange, InfAxes, InfRanges
 import LinearAlgebra: matprod, qr, AbstractTriangular, AbstractQ, adjoint, transpose, AdjOrTrans, copymutable
 import LazyArrays: applybroadcaststyle, CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, LazyArrayStyle,
                     resizedata!, MemoryLayout, AbstractLazyLayout,
                     factorize, sub_materialize, LazyLayout, LazyArrayStyle, layout_getindex,
                     applylayout, ApplyLayout, PaddedLayout, CachedLayout, AbstractCachedVector, AbstractCachedMatrix, cacheddata, zero!, MulAddStyle, ApplyArray,
                     LazyArray, LazyMatrix, LazyVector, paddeddata, arguments, resizedata!, simplifiable, simplify, LazyLayouts
-import MatrixFactorizations: ul, ul!, _ul, ql, ql!, _ql, QLPackedQ, getL, getR, getQ, getU, reflector!, reflectorApply!, QL, QR, QRPackedQ,
+import MatrixFactorizations: ul, ul!, ul_layout, ql, ql!, ql_layout, reversecholesky_layout, QLPackedQ, getL, getR, getQ, getU, reflector!, reflectorApply!, QL, QR, QRPackedQ,
                             QRPackedQLayout, AdjQRPackedQLayout, QLPackedQLayout, AdjQLPackedQLayout, LayoutQ, copymutable_size
 
 import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange, blockcolsupport, BlockLayout, AbstractBlockLayout, BlockSlice
@@ -121,6 +121,7 @@ include("banded/hessenbergq.jl")
 include("banded/infbanded.jl")
 include("blockbanded/blockbanded.jl")
 include("banded/infqltoeplitz.jl")
+include("banded/infreversecholeskytoeplitz.jl")
 include("infql.jl")
 include("infqr.jl")
 include("inful.jl")

src/banded/infbanded.jl

@@ -34,7 +34,7 @@ sub_materialize(_, V::SubArray{<:Any,1,<:AbstractMatrix,Tuple{InfBandCartesianIn
     _inf_banded_sub_materialize(MemoryLayout(parent(V)), V)
 
 const TriToeplitz{T} = Tridiagonal{T,Fill{T,1,Tuple{OneToInf{Int}}}}
-const ConstRowMatrix{T} = ApplyMatrix{T,typeof(*),<:Tuple{<:AbstractVector,<:AbstractFill{<:Any,2,Tuple{OneTo{Int},OneToInf{Int}}}}}
+const ConstRowMatrix{T} = ApplyMatrix{T,typeof(*),<:Tuple{<:AbstractVector,<:AbstractFillMatrix{<:Any,Tuple{OneTo{Int},OneToInf{Int}}}}}
 const PertConstRowMatrix{T} = Hcat{T,<:Tuple{Array{T},<:ConstRowMatrix{T}}}
 const InfToeplitz{T} = BandedMatrix{T,<:ConstRowMatrix{T},OneToInf{Int}}
 const PertToeplitz{T} = BandedMatrix{T,<:PertConstRowMatrix{T},OneToInf{Int}}
@@ -337,10 +337,26 @@ for Typ in (:ConstRows, :PertConstRows)
     end
 end
 
+"""
+    TridiagonalToeplitzLayout
+
+represents a matrix which is tridiagonal and toeplitz. Must support
+`subdiagonalconstant`, `diagonalconstant`, `supdiagonalconstant`.
+"""
 struct TridiagonalToeplitzLayout <: AbstractLazyBandedLayout end
 const BandedToeplitzLayout = BandedColumns{ConstRows}
 const PertToeplitzLayout = BandedColumns{PertConstRows}
 const PertTriangularToeplitzLayout{UPLO,UNIT} = TriangularLayout{UPLO,UNIT,BandedColumns{PertConstRows}}
+struct BidiagonalToeplitzLayout <: AbstractLazyBandedLayout end
+struct PertBidiagonalToeplitzLayout <: AbstractLazyBandedLayout end
+struct PertTridiagonalToeplitzLayout <: AbstractLazyBandedLayout end
+
+const InfToeplitzLayouts = Union{TridiagonalToeplitzLayout, BandedToeplitzLayout, BidiagonalToeplitzLayout,
+                                 PertToeplitzLayout, PertTriangularToeplitzLayout, PertBidiagonalToeplitzLayout, PertTridiagonalToeplitzLayout}
+
+subdiagonalconstant(A) = getindex_value(subdiagonaldata(A))
+diagonalconstant(A) = getindex_value(diagonaldata(A))
+supdiagonalconstant(A) = getindex_value(supdiagonaldata(A))
 
 
 _BandedMatrix(::BandedToeplitzLayout, A::AbstractMatrix) =
@@ -428,18 +444,13 @@ _bandedfill_mul(M::MulAdd, ::Tuple{Any,InfAxes}, ::Tuple{InfAxes,InfAxes}) = App
 _bandedfill_mul(M::MulAdd, ::Tuple{InfAxes,InfAxes}, ::Tuple{InfAxes,Any}) = ApplyArray(*, M.A, M.B)
 _bandedfill_mul(M::MulAdd, ::Tuple{Any,InfAxes}, ::Tuple{InfAxes,Any}) = ApplyArray(*, M.A, M.B)
 
-mulreduce(M::Mul{BandedToeplitzLayout, BandedToeplitzLayout}) = ApplyArray(M)
-mulreduce(M::Mul{BandedToeplitzLayout}) = ApplyArray(M)
-mulreduce(M::Mul{BandedToeplitzLayout,<:PaddedLayout}) = MulAdd(M)
-mulreduce(M::Mul{<:Any, BandedToeplitzLayout}) = ApplyArray(M)
-mulreduce(M::Mul{<:BandedColumns{<:AbstractFillLayout}, PertToeplitzLayout}) = ApplyArray(M)
-mulreduce(M::Mul{<:PertToeplitzLayout, <:BandedColumns{<:AbstractFillLayout}}) = ApplyArray(M)
-mulreduce(M::Mul{<:BandedColumns{<:AbstractFillLayout}, BandedToeplitzLayout}) = ApplyArray(M)
-mulreduce(M::Mul{BandedToeplitzLayout, <:BandedColumns{<:AbstractFillLayout}}) = ApplyArray(M)
-mulreduce(M::Mul{<:AbstractQLayout, BandedToeplitzLayout}) = ApplyArray(M)
-mulreduce(M::Mul{<:AbstractQLayout, PertToeplitzLayout}) = ApplyArray(M)
-mulreduce(M::Mul{<:DiagonalLayout, BandedToeplitzLayout}) = Lmul(M)
-mulreduce(M::Mul{BandedToeplitzLayout, <:DiagonalLayout}) = Rmul(M)
+mulreduce(M::Mul{<:InfToeplitzLayouts, <:InfToeplitzLayouts}) = ApplyArray(M)
+mulreduce(M::Mul{<:InfToeplitzLayouts}) = ApplyArray(M)
+mulreduce(M::Mul{<:InfToeplitzLayouts,<:PaddedLayout}) = MulAdd(M)
+mulreduce(M::Mul{<:Any, <:InfToeplitzLayouts}) = ApplyArray(M)
+mulreduce(M::Mul{<:AbstractQLayout, <:InfToeplitzLayouts}) = ApplyArray(M)
+mulreduce(M::Mul{<:DiagonalLayout, <:InfToeplitzLayouts}) = Lmul(M)
+mulreduce(M::Mul{<:InfToeplitzLayouts, <:DiagonalLayout}) = Rmul(M)
 
 
 function _bidiag_forwardsub!(M::Ldiv{<:Any,<:PaddedLayout})
@@ -491,8 +502,6 @@ for Typ in (:(LinearAlgebra.Tridiagonal{<:Any,<:InfFill}),
     end
 end
 
-struct BidiagonalToeplitzLayout <: AbstractLazyBandedLayout end
-
 for Typ in (:(LinearAlgebra.Bidiagonal{<:Any,<:InfFill}),
             :(LazyBandedMatrices.Bidiagonal{<:Any,<:InfFill,<:InfFill}))
     @eval begin
@@ -501,6 +510,9 @@ for Typ in (:(LinearAlgebra.Bidiagonal{<:Any,<:InfFill}),
     end
 end
 
+*(A::LinearAlgebra.Bidiagonal{<:Any,<:InfFill}, B::LinearAlgebra.Bidiagonal{<:Any,<:InfFill}) =
+    mul(A, B)
+
 # fall back for Ldiv
 triangularlayout(::Type{<:TriangularLayout{UPLO,'N'}}, ::TridiagonalToeplitzLayout) where UPLO = BidiagonalToeplitzLayout()
 materialize!(L::MatLdivVec{BidiagonalToeplitzLayout,Lay}) where Lay = materialize!(Ldiv{BidiagonalLayout{FillLayout,FillLayout},Lay}(L.A, L.B))
@@ -518,3 +530,10 @@ copy(A::Transpose{T,<:BandedMatrix{T,<:Any,OneToInf{Int}}}) where T = transpose(
 
 Base.typed_hcat(::Type{T}, A::BandedMatrix{<:Any,<:Any,OneToInf{Int}}, B::AbstractVecOrMat...) where T = Hcat{T}(A, B...)
 
+
+
+###
+# SymTriPertToeplitz
+###
+
+MemoryLayout(::Type{<:SymTriPertToeplitz}) = PertTridiagonalToeplitzLayout()

src/banded/infqltoeplitz.jl

@@ -41,8 +41,8 @@ end
 
 
 
-function ql(Op::TriToeplitz{T}; kwds...) where {T<:Real}
-    Z, A, B = Op.dl.value, Op.d.value, Op.du.value
+function ql_layout(::TridiagonalToeplitzLayout, ::NTuple{2,OneToInf{Int}}, Op::AbstractMatrix{T}, args...; kwds...) where {T<:Real}
+    Z, A, B = subdiagonalconstant(Op), diagonalconstant(Op), supdiagonalconstant(Op)
     d, e = tail_de([Z, A, B]; kwds...) # fixed point of QL but with two Qs, one that changes sign
     X = [Z A B; zero(T) d e]
     F = ql_X!(X)

src/banded/infreversecholeskytoeplitz.jl

@@ -0,0 +1,44 @@
+function reversecholesky_layout(::TridiagonalToeplitzLayout, ::NTuple{2,OneToInf{Int}}, A, ::NoPivot; kwds...)
+    # [a b     = [α β       * [α
+    #  b a ⋱        ⋱ ⋱]     β α
+    #    ⋱ ⋱]
+    #  [a       b'      [a    b'        [1  b'inv(U)'   ]   [a - b'inv(U)'inv(U)*b      [1
+    #   b    A   ] =     b U*U']    =       U           ]                           I]   inv(U)b    U']
+
+    # since inv(U) = [inv(α) …] we have inv(U)*(b*e₁) = b/α
+    # we also assume Toeplitz structure so that α^2 = a - b'inv(U)'inv(U)*b = a - b^2/α^2
+    # Thus we have a Quadratic equation:
+    #
+    #   α^4 - a*α^2 + b^2 = 0
+    #
+    # i.e. α^2 = (a ± sqrt(a^2 - 4b^2))/2.
+    # We also have αβ = b.
+
+    a = diagonalconstant(A)
+    b = supdiagonalconstant(A)
+    @assert b == subdiagonalconstant(A)
+
+    α² = (a + sqrt(a^2 - 4b^2))/2
+    # (a - sqrt(a^2 - 4b^2))/2 other branch gives non-invertible U
+    
+    α = sqrt(α²)
+    β = b/α
+
+    ReverseCholesky(Bidiagonal(Fill(α,∞), Fill(β,∞), :U), 'U', 0)
+end
+
+function reversecholesky_layout(::PertTridiagonalToeplitzLayout, ::NTuple{2,OneToInf{Int}}, A, ::NoPivot; kwds...)
+    a = diagonaldata(A)
+    aₙ,a∞ = arguments(vcat, a)
+    b = supdiagonaldata(A)
+    bₙ,b∞ = arguments(vcat, b)
+    U∞, = reversecholesky(SymTridiagonal(a∞, b∞))
+
+    n = max(length(aₙ), length(bₙ)+1)
+    Aₙ = SymTridiagonal([aₙ; float(a∞[1:(n-length(aₙ))])], [bₙ; float(b∞[1:(n-length(bₙ)-1)])])
+    α = U∞[1,1]
+    b = getindex_value(b∞)
+    Aₙ[end,end] -= b^2/α^2
+    Uₙ, = reversecholesky(Aₙ)
+    ReverseCholesky(Bidiagonal([Uₙ.dv; U∞.dv], [Uₙ.ev; U∞.ev], :U), 'U', 0)
+end

src/infql.jl

@@ -281,7 +281,7 @@ function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N',BandedColumns{Per
     b
 end
 
-_ql(layout, ::NTuple{2,OneToInf{Int}}, A, args...; kwds...) = error("Not implemented")
+ql_layout(layout, ::NTuple{2,OneToInf{Int}}, A, args...; kwds...) = error("Not implemented")
 
 _data_tail(::PaddedLayout, a) = paddeddata(a), zero(eltype(a))
 _data_tail(::AbstractFillLayout, a) = Vector{eltype(a)}(), getindex_value(a)
@@ -293,7 +293,7 @@ function _data_tail(::ApplyLayout{typeof(vcat)}, a)
 end
 _data_tail(a) = _data_tail(MemoryLayout(a), a)
 
-function _ql(::SymTridiagonalLayout, ::NTuple{2,OneToInf{Int}}, A, args...; kwds...)
+function ql_layout(::Union{PertTridiagonalToeplitzLayout,SymTridiagonalLayout}, ::NTuple{2,OneToInf{Int}}, A, args...; kwds...)
     T = eltype(A)
     d,d∞ = _data_tail(A.dv)
     ev,ev∞ = _data_tail(A.ev)
@@ -312,7 +312,7 @@ end
 
 # TODO: This should be redesigned as ql(BandedMatrix(A))
 # But we need to support dispatch on axes
-function _ql(::TridiagonalLayout, ::NTuple{2,OneToInf{Int}}, A, args...; kwds...)
+function ql_layout(::TridiagonalLayout, ::NTuple{2,OneToInf{Int}}, A, args...; kwds...)
     T = eltype(A)
     d,d∞ = _data_tail(A.d)
     dl,dl∞ = _data_tail(A.dl)

src/inful.jl

@@ -34,7 +34,7 @@ function _ultailL1(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
     C*(V*Diagonal(inv.(λs))/V)
 end
 
-function _ul(::TridiagonalToeplitzLayout, J::AbstractMatrix, ::Val{false}; check::Bool = true)
+function ul_layout(::TridiagonalToeplitzLayout, J::AbstractMatrix, ::Val{false}; check::Bool = true)
     C = getindex_value(subdiagonaldata(J))
     A = getindex_value(diagonaldata(J))
     B = getindex_value(supdiagonaldata(J))
@@ -43,15 +43,15 @@ function _ul(::TridiagonalToeplitzLayout, J::AbstractMatrix, ::Val{false}; check
     UL(Tridiagonal(Fill(convert(typeof(L),C),∞), Fill(L,∞), Fill(U,∞)), OneToInf(), 0)
 end
 
-function _ul(::TridiagonalToeplitzLayout, J::AbstractMatrix, ::Val{true}; check::Bool = true)
+function ul_layout(::TridiagonalToeplitzLayout, J::AbstractMatrix, ::Val{true}; check::Bool = true)
     C = getindex_value(subdiagonaldata(J))
     A = getindex_value(diagonaldata(J))
     B = getindex_value(supdiagonaldata(J))
     A^2 ≥ 4B*C || error("Pivotting not implemented")
     ul(J, Val(false))
 end
 
-function _ul(::BlockTridiagonalToeplitzLayout, J::AbstractMatrix, ::Val{false}; check::Bool = true)
+function ul_layout(::BlockTridiagonalToeplitzLayout, J::AbstractMatrix, ::Val{false}; check::Bool = true)
     C = getindex_value(subdiagonaldata(blocks(J)))
     A = getindex_value(diagonaldata(blocks(J)))
     B = getindex_value(supdiagonaldata(blocks(J)))

test/runtests.jl

@@ -15,12 +15,7 @@ import LazyBandedMatrices: BroadcastBandedBlockBandedLayout, BroadcastBandedLayo
 
 using Aqua
 @testset "Project quality" begin
-    Aqua.test_all(InfiniteLinearAlgebra, ambiguities=false, unbound_args=false, piracy=false,
-        # Project.toml formatting issue on v1.6
-        # Pkg issue: https://github.com/JuliaLang/Pkg.jl/issues/3481
-        # Aqua workaround: https://github.com/JuliaTesting/Aqua.jl/issues/105#issuecomment-1551405866
-        # we only check the formatting on more recent versions
-        project_toml_formatting = VERSION>=v"1.7")
+    Aqua.test_all(InfiniteLinearAlgebra, ambiguities=false, unbound_args=false, piracies=false)
 end
 
 @testset "chop" begin
@@ -477,3 +472,4 @@ include("test_infql.jl")
 include("test_infqr.jl")
 include("test_inful.jl")
 include("test_infcholesky.jl")
+include("test_infreversecholesky.jl")

test/test_infql.jl

@@ -220,7 +220,7 @@ using ArrayLayouts: TriangularLayout, UnknownLayout
             Aplain = LinearAlgebra.Tridiagonal([[1., 2.]; Fill(1.,∞)], [[1.,2.]; Fill(3.,∞)], [[1., 2.]; Fill(1.,∞)])
             Qsym, Lsym = ql(Asym, 1e-10)
             Qplain, Lplain = ql(Aplain, 1e-10)
-            
+
             @test size(Qsym) == (ℵ₀, ℵ₀)
             @test size(Lsym) == (ℵ₀, ℵ₀)
             @test size(Qplain) == (ℵ₀, ℵ₀)
@@ -250,7 +250,7 @@ using ArrayLayouts: TriangularLayout, UnknownLayout
             A = im * LinearAlgebra.Tridiagonal([[1., 2.]; Fill(1.,∞)], [[1.,2.]; Fill(3.,∞)], [[1., 2.]; Fill(1.,∞)])
             Abanded = _BandedMatrix(conj.(Hcat(Vcat(1.,A.du),A.d,A.dl)'), ℵ₀, 1, 1)
             F = ql(Abanded)
-            @test (F.Q[1:51,1:51]*F.L[1:51,1:51])[1:50,1:50] ≈ A[1:50,1:50] 
+            @test (F.Q[1:51,1:51]*F.L[1:51,1:51])[1:50,1:50] ≈ A[1:50,1:50]
             @test MemoryLayout(F.L.data) == LazyBandedLayout()
             @test bandwidths(F.L) == (2,0)
         end
@@ -273,9 +273,15 @@ using ArrayLayouts: TriangularLayout, UnknownLayout
                              Vcat([[-1. 1.; 1. 1.]], Fill([-1. 1.; 1. 1.], ∞)),
                              Vcat([[0. 0.; 1. 0.]], Fill([0. 0.; 1. 0.], ∞)))
 
-        A[1,1] = 2 
+        A[1,1] = 2
 
         x = -0.95
         @test ql(A-x*I).L[1,1] isa Float64
     end
+
+    @testset "SymTridiagonal QL" begin
+        A = SymTridiagonal(Fill(3, ∞), Fill(1, ∞))
+        Q,L = ql(A)
+        @test (Q*L)[1:10,1:10] ≈ A[1:10,1:10]
+    end
 end

test/test_infreversecholesky.jl

@@ -0,0 +1,16 @@
+using InfiniteLinearAlgebra, MatrixFactorizations, ArrayLayouts, LinearAlgebra, Test
+
+
+
+@testset "infreversecholesky" begin
+    @testset "Tri Toeplitz" begin
+        A = SymTridiagonal(Fill(3, ∞), Fill(1, ∞))
+        U, = reversecholesky(A)
+        @test (U*U')[1:10,1:10] ≈ A[1:10,1:10]
+    end
+
+    @testset "Pert Tri Toeplitz" begin
+        A = SymTridiagonal([[4,5, 6]; Fill(3, ∞)], [[2,3]; Fill(1, ∞)])
+        @test reversecholesky(A).U[1:100,1:100] ≈ reversecholesky(A[1:1000,1:1000]).U[1:100,1:100]
+    end
+end