SymTridiagonal ql (#112)

* Support for SymTridiagonal in ∞-QL

* Update Project.toml

* Update test_infqr.jl

* add tests

* more tests

* Update test_infqr.jl

* Update test_infql.jl

Author: dlfivefifty

Project.toml

@@ -17,16 +17,16 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 
 [compat]
-ArrayLayouts = "0.8.9"
+ArrayLayouts = "0.8.11"
 BandedMatrices = "0.17"
 BlockArrays = "0.16.14"
 BlockBandedMatrices = "0.11.5"
 DSP = "0.7"
 FillArrays = "0.13"
 InfiniteArrays = "0.12"
 LazyArrays = "0.22"
-LazyBandedMatrices = "0.8"
-MatrixFactorizations = "0.9"
+LazyBandedMatrices = "0.8.1"
+MatrixFactorizations = "0.9.2"
 SemiseparableMatrices = "0.3"
 julia = "1.6"
 

examples/toeplitz.jl

@@ -11,7 +11,7 @@ using InfiniteLinearAlgebra, BandedMatrices, PyPlot
 # Basic routines for plotting
 ###
 
-function ℓ11(A,λ; kwds...) 
+function ℓ11(A,λ; kwds...)
     try 
         abs(ql(A-λ*I; kwds...).L[1,1]) 
     catch DomainError 

src/InfiniteLinearAlgebra.jl

@@ -17,10 +17,10 @@ import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row
 import InfiniteArrays: OneToInf, InfUnitRange, Infinity, PosInfinity, InfiniteCardinal, InfStepRange, AbstractInfUnitRange, InfAxes, InfRanges
 import LinearAlgebra: matprod, qr, AbstractTriangular, AbstractQ, adjoint, transpose, AdjOrTrans
 import LazyArrays: applybroadcaststyle, CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, LazyArrayStyle,
-                    resizedata!, MemoryLayout,
+                    resizedata!, MemoryLayout, most,
                     factorize, sub_materialize, LazyLayout, LazyArrayStyle, layout_getindex,
-                    applylayout, ApplyLayout, PaddedLayout, zero!, MulAddStyle,
-                    LazyArray, LazyMatrix, LazyVector, paddeddata
+                    applylayout, ApplyLayout, PaddedLayout, CachedLayout, cacheddata, zero!, MulAddStyle,
+                    LazyArray, LazyMatrix, LazyVector, paddeddata, arguments
 import MatrixFactorizations: ul, ul!, _ul, ql, ql!, _ql, QLPackedQ, getL, getR, getU, reflector!, reflectorApply!, QL, QR, QRPackedQ,
                             QRPackedQLayout, AdjQRPackedQLayout, QLPackedQLayout, AdjQLPackedQLayout, LayoutQ
 

src/banded/infqltoeplitz.jl

@@ -123,6 +123,12 @@ end
 mul(A::ProductQ, x::AbstractVector) = _productq_mul(A, x)
 
 
+
+mul(Q::ProductQ, X::AbstractMatrix) = ApplyArray(*, Q.Qs...) * X
+mul(X::AbstractMatrix, Q::ProductQ) = X * ApplyArray(*, Q.Qs...)
+
+
+
 # LQ where Q is a product of orthogonal operations
 struct QLProduct{T,QQ<:Tuple,LL} <: Factorization{T}
     Qs::QQ

src/infql.jl

@@ -52,15 +52,11 @@ function qltail(Z::Number, A::Number, B::Number)
 
     e = sqrt(n^2 - abs2(B))
     d = σ*e*Z/n
-
-    Q = 
+ 
     ql!([Z A B; 0 d e])
 end
 
 
-ql(A::SymTriPertToeplitz{T}; kwds...) where T = ql_hessenberg!(BandedMatrix(A, (bandwidth(A,1)+bandwidth(A,2),bandwidth(A,2))); kwds...)
-ql(A::SymTridiagonal{T}; kwds...) where T = ql_hessenberg!(BandedMatrix(A, (bandwidth(A,1)+bandwidth(A,2),bandwidth(A,2))); kwds...)
-ql(A::TriPertToeplitz{T}; kwds...) where T = ql_hessenberg!(BandedMatrix(A, (bandwidth(A,1)+bandwidth(A,2),bandwidth(A,2))); kwds...)
 ql_hessenberg(A::InfBandedMatrix{T}; kwds...) where T = ql_hessenberg!(BandedMatrix(A, (bandwidth(A,1)+bandwidth(A,2),bandwidth(A,2))); kwds...)
 
 toeptail(B::BandedMatrix{T}) where T = 
@@ -303,3 +299,78 @@ function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N',BandedColumns{Per
     end
     b
 end
+
+_ql(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)
+_data_tail(::CachedLayout, a) = cacheddata(a), getindex_value(a.array)
+function _data_tail(::ApplyLayout{typeof(vcat)}, a)
+    args = arguments(vcat, a)
+    dat,tl = _data_tail(last(args))
+    vcat(most(args)..., dat), tl
+end
+_data_tail(a) = _data_tail(MemoryLayout(a), a)
+
+function _ql(::SymTridiagonalLayout, ::NTuple{2,OneToInf{Int}}, A, args...; kwds...)
+    T = eltype(A)
+    d,d∞ = _data_tail(A.dv)
+    ev,ev∞ = _data_tail(A.ev)
+    
+    m = max(length(d), length(ev)+1)
+    dat = zeros(T, 3, m)
+    dat[1,2:1+length(ev)] .= ev
+    dat[1,2+length(ev):end] .= ev∞
+    dat[2,1:length(d)] .= d
+    dat[2,1+length(d):end] .= d∞
+    dat[3,1:length(ev)] .= ev
+    dat[3,1+length(ev):end] .= ev∞
+
+    ql(_BandedMatrix(Hcat(dat, [ev∞,d∞,ev∞] * Ones{T}(1,∞)), ℵ₀, 1, 1), args...; kwds...)
+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...)
+    T = eltype(A)
+    d,d∞ = _data_tail(A.d)
+    dl,dl∞ = _data_tail(A.dl)
+    du,du∞ = _data_tail(A.du)
+    
+    m = max(length(d), length(du)+1, length(dl))
+    dat = zeros(T, 3, m)
+    dat[1,2:1+length(du)] .= du
+    dat[1,2+length(du):end] .= du∞
+    dat[2,1:length(d)] .= d
+    dat[2,1+length(d):end] .= d∞
+    dat[3,1:length(dl)] .= dl
+    dat[3,1+length(dl):end] .= dl∞
+
+    ql(_BandedMatrix(Hcat(dat, [du∞,d∞,dl∞] * Ones{T}(1,∞)), ℵ₀, 1, 1), args...; kwds...)
+end
+
+
+###
+# L*Q special case
+###
+
+copy(M::Mul{TriangularLayout{'L', 'N', PertToeplitzLayout}, HessenbergQLayout{'L'}}) =
+    ApplyArray(*, M.A, M.B)
+
+copy(M::Mul{HessenbergQLayout{'L'}, TriangularLayout{'L', 'N', PertToeplitzLayout}}) =
+    ApplyArray(*, M.A, M.B)
+
+
+function LazyBandedMatrices._SymTridiagonal(::Tuple{TriangularLayout{'L', 'N', PertToeplitzLayout}, HessenbergQLayout{'L'}}, A)
+    T = eltype(A)
+    L,Q = arguments(*, A)
+    Ldat,L∞ = arguments(hcat, L.data.data)
+    Qdat, Q∞ = arguments(vcat, Q.q)
+
+    m = max(size(Ldat,2)+2, length(Qdat)+1)
+    dv = [A[k,k] for k=1:m]
+    ev = [A[k,k+1] for k=1:m-1]
+    SymTridiagonal([dv; Fill(dv[end],∞)], [ev; Fill(ev[end],∞)])
+end

src/infqr.jl

@@ -7,8 +7,9 @@ end
 
 function AdaptiveQRData(::Union{SymmetricLayout{<:AbstractBandedLayout},AbstractBandedLayout}, A::AbstractMatrix{T}) where T
     l,u = bandwidths(A)
-    data = BandedMatrix{T}(undef,(2l+u+1,0),(l,l+u)) # pad super
-    AdaptiveQRData(CachedArray(data,A), Vector{T}(), 0)
+    FT = float(T)
+    data = BandedMatrix{FT}(undef,(2l+u+1,0),(l,l+u)) # pad super
+    AdaptiveQRData(CachedArray(data,A), Vector{FT}(), 0)
 end
 
 function AdaptiveQRData(::AbstractAlmostBandedLayout, A::AbstractMatrix{T}) where T
@@ -123,8 +124,8 @@ function getindex(F::AdaptiveQRFactors, k::Int, j::Int)
     F.data.data[k,j]
 end
 
-colsupport(F::QRPackedQ{<:Any,<:AdaptiveQRFactors}, j) = colsupport(F.factors, j)
-rowsupport(F::QRPackedQ{<:Any,<:AdaptiveQRFactors}, j) = rowsupport(F.factors, j)
+colsupport(F::QRPackedQ{<:Any,<:AdaptiveQRFactors}, j) = 1:last(colsupport(F.factors, j))
+rowsupport(F::QRPackedQ{<:Any,<:AdaptiveQRFactors}, j) = first(rowsupport(F.factors, j)):size(F,2)
 
 blockcolsupport(F::QRPackedQ{<:Any,<:AdaptiveQRFactors}, j) = blockcolsupport(F.factors, j)
 
@@ -349,4 +350,7 @@ ldiv!(F::QR{<:Any,<:AdaptiveQRFactors}, b::LayoutVector; kwds...) = ldiv!(F.R, l
 
 
 factorize(A::BandedMatrix{<:Any,<:Any,<:OneToInf}) = qr(A)
-qr(A::SymTridiagonal{T,<:AbstractFill{T,1,Tuple{OneToInf{Int}}}}) where T = adaptiveqr(A)
+qr(A::SymTridiagonal{T,<:AbstractFill{T,1,Tuple{OneToInf{Int}}}}) where T = adaptiveqr(A)
+
+copy(M::Mul{<:QRPackedQLayout{<:AdaptiveLayout}}) = ApplyArray(*, M.A, M.B)
+copy(M::Mul{<:Any,<:QRPackedQLayout{<:AdaptiveLayout}}) = ApplyArray(*, M.A, M.B)

test/test_infql.jl

@@ -179,4 +179,22 @@ import BandedMatrices: _BandedMatrix
         @test Q[1:10,1:10] ≈ Qn[1:10,1:10] * diagm(0 => [1; -Ones(9)] )
         @test (Q'A)[1:10,1:10] ≈ diagm(0 => [1; -Ones(9)] ) * Ln[1:10,1:10] ≈ L[1:10,1:10]
     end
+
+    @testset "Tridiagonal QL" begin
+        for A in (LinearAlgebra.SymTridiagonal([[1,2]; Fill(3,∞)], [[1, 2]; Fill(1,∞)]),
+                LinearAlgebra.Tridiagonal([[1, 2]; Fill(1,∞)], [[1,2]; Fill(3,∞)], [[1, 2]; Fill(1,∞)]),
+                LazyBandedMatrices.SymTridiagonal([[1,2]; Fill(3,∞)], [[1, 2]; Zeros(∞)]),
+                LazyBandedMatrices.SymTridiagonal([[1,2]; Fill(3,∞)], [[1, 2]; zeros(∞)]),
+                LazyBandedMatrices.SymTridiagonal([[1,2]; fill(3,∞)], [[1, 2]; zeros(∞)]))
+            @test abs.(ql(A).L[1:10,1:10]) ≈ abs.(ql(A[1:1000,1:1000]).L[1:10,1:10])
+        end
+
+        A = LazyBandedMatrices.SymTridiagonal([[1,2]; Fill(3,∞)], [[1, 2]; Fill(1,∞)])
+        Q,L = ql(A)
+        @test (Q*L)[1:10,1:10] ≈ A[1:10,1:10]
+
+        @test (L*Q)[1:10,1:10] ≈ LazyBandedMatrices.SymTridiagonal(L*Q)[1:10,1:10]
+    end
+
+    @test_throws ErrorException ql(zeros(∞,∞))
 end

test/test_infqr.jl

@@ -43,23 +43,23 @@ import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout
 
         @testset "col/rowsupport" begin
             A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ℵ₀, 1, 1)
-            F = qr(A);
+            F = qr(A)
             @test MemoryLayout(typeof(F.factors)) isa AdaptiveLayout{BandedColumns{DenseColumnMajor}}
             @test bandwidths(F.factors) == (1,2)
             @test colsupport(F.factors,1) ==  1:2
             @test colsupport(F.factors,5) ==  3:6
             @test rowsupport(F.factors,1) ==  1:3
             @test rowsupport(F.factors,5) ==  4:7
-            Q,R = F;
+            Q,R = F
             @test MemoryLayout(typeof(R)) isa TriangularLayout
             @test colsupport(R,1) == 1:1
             @test colsupport(R,5) == 3:5
             @test rowsupport(R,1) == 1:3
             @test rowsupport(R,5) == 5:7
             @test colsupport(Q,1) == 1:2
-            @test colsupport(Q,5) == 3:6
-            @test rowsupport(Q,1) == 1:3
-            @test rowsupport(Q,5) == 4:7
+            @test colsupport(Q,5) == 1:6
+            @test rowsupport(Q,1) == 1:∞
+            @test rowsupport(Q,5) == 4:∞
         end
 
         @testset "Qmul" begin
@@ -296,4 +296,11 @@ import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout
         b = [1; 2; 3; zeros(∞)]
         @test (qr(A) \ b) ≈ (ul(A) \ b)
     end
+
+    @testset "SymTridiagonal QR" begin
+        A = LazyBandedMatrices.SymTridiagonal([[1,2]; Fill(3,∞)], [[1, 2]; Fill(1,∞)])
+        Q,R = qr(A)
+        @test (Q*R)[1:10,1:10] ≈ A[1:10,1:10]
+        @test (R*Q)[1:10,1:10] ≈ LazyBandedMatrices.SymTridiagonal((R*Q)[1:10,1:10])
+    end
 end