Suport ∞-UL decomposition for BlockToeplitz (#35)

* Suport ∞-UL decomposition for BlockToeplitz

* add tests

* Support adaptive QR

* Use MatrixFactorizations 0.5

* tridiagonal Toeplitz ul

* compare with adaptiveqr

* Update Project.toml

* Drop Julia 1.3

* avoid special lazy BandedMatrix overload

* Update Project.toml

* ArrayLayouts 0.4

* update for mulreduce

* using works

* Update runtests.jl

* Update runtests.jl

* Update Project.toml

* tests pass again

* Update Project.toml

Author: dlfivefifty

.travis.yml

@@ -5,7 +5,6 @@ os:
   - osx
   - windows
 julia:
-  - 1.3
   - 1.5
   - nightly
 matrix:

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.3.5"
+version = "0.4.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -16,17 +16,18 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 
 [compat]
-ArrayLayouts = "0.3.5"
-BandedMatrices = "0.15.13"
-BlockArrays = "0.12.8"
-BlockBandedMatrices = "0.8.7"
-FillArrays = "0.8.11, 0.9"
-InfiniteArrays = "0.7.3"
-LazyArrays = "0.16.13"
-LazyBandedMatrices = "0.2.11"
-MatrixFactorizations = "0.4.1, 0.5"
-SemiseparableMatrices = "0.1"
-julia = "1.3"
+ArrayLayouts = "0.4.3"
+BandedMatrices = "0.15.17"
+BlockArrays = "0.12.12"
+BlockBandedMatrices = "0.9"
+FillArrays = "0.9.4"
+InfiniteArrays = "0.8"
+LazyArrays = "0.17.3"
+LazyBandedMatrices = "0.3.1"
+MatrixFactorizations = "0.5.2"
+SemiseparableMatrices = "0.1.2"
+julia = "1.5"
+
 
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

src/InfiniteLinearAlgebra.jl

@@ -4,25 +4,26 @@ using BlockArrays, BlockBandedMatrices, BandedMatrices, LazyArrays, LazyBandedMa
 
 import Base: +, -, *, /, \, ^, OneTo, getindex, promote_op, _unsafe_getindex, size, axes,
             AbstractMatrix, AbstractArray, Matrix, Array, Vector, AbstractVector, Slice,
-            show, getproperty, copy, map, require_one_based_indexing, similar
+            show, getproperty, copy, map, require_one_based_indexing, similar, inv
 import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted
 
-import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, sublayout, _qr,
-                        MatLmulVec, MatLmulMat, AbstractQLayout, materialize!
+import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, TridiagonalLayout, 
+                        sublayout, _qr, __qr, MatLmulVec, MatLmulMat, AbstractQLayout, materialize!, diagonaldata, subdiagonaldata, supdiagonaldata,
+                        _bidiag_forwardsub!, mulreduce
 import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns,
-                        _default_banded_broadcast
+                        _default_banded_broadcast, banded_similar
 import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row
 import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange, AbstractInfUnitRange, InfAxes
 import LinearAlgebra: lmul!,  ldiv!, matprod, qr, AbstractTriangular, AbstractQ, adjoint, transpose
 import LazyArrays: applybroadcaststyle, CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, ApplyStyle, LazyArrayApplyStyle, LazyArrayStyle,
-                    resizedata!, MemoryLayout, mulapplystyle, LmulStyle, RmulStyle,
+                    resizedata!, MemoryLayout,
                     factorize, sub_materialize, LazyLayout, LazyArrayStyle, layout_getindex,
-                    @lazymul, applylayout, ApplyLayout, PaddedLayout, zero!, MulAddStyle,
+                    applylayout, ApplyLayout, PaddedLayout, zero!, MulAddStyle,
                     LazyArray, LazyMatrix, LazyVector, paddeddata
-import MatrixFactorizations: ql, ql!, _ql, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ,
+import MatrixFactorizations: ul, ul!, _ul, ql, ql!, _ql, QLPackedQ, getL, getR, getU, reflector!, reflectorApply!, QL, QR, QRPackedQ,
                             QRPackedQLayout, AdjQRPackedQLayout, QLPackedQLayout, AdjQLPackedQLayout, LayoutQ
 
-import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange, blockcolsupport
+import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange, blockcolsupport, BlockLayout, AbstractBlockLayout
 
 import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!, _banded_qr, _BandedMatrix
 
@@ -34,9 +35,6 @@ import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, _BlockSkylineMat
 
 import SemiseparableMatrices: AbstractAlmostBandedLayout, _almostbanded_qr!
 
-LazyArrays.@lazymul BandedMatrix{<:Any,<:Any,<:OneToInf}
-*(A::BandedMatrix{<:Any,<:Any,<:OneToInf}, b::CachedVector) = apply(*,A,b)
-
 
 # BroadcastStyle(::Type{<:BandedMatrix{<:Any,<:Any,<:OneToInf}}) = LazyArrayStyle{2}()
 
@@ -61,6 +59,7 @@ include("blockbanded/blockbanded.jl")
 include("banded/infqltoeplitz.jl")
 include("infql.jl")
 include("infqr.jl")
+include("inful.jl")
 
 #######
 # block broadcasted

src/banded/infbanded.jl

@@ -20,6 +20,8 @@ _prepad(p, a::Zeros{T,1}) where T = Zeros{T}(length(a)+p)
 _prepad(p, a::Ones{T,1}) where T = Ones{T}(length(a)+p)
 _prepad(p, a::AbstractFill{T,1}) where T = Fill{T}(getindex_value(a), length(a)+p)
 
+banded_similar(T, (m,n)::Tuple{Int,Infinity}, (l,u)::Tuple{Int,Int}) = BandedMatrix{T}(undef, (n,m), (u,l))'
+
 function BandedMatrix{T}(kv::Tuple{Vararg{Pair{<:Integer,<:AbstractVector}}},
                          ::NTuple{2,Infinity},
                          (l,u)::NTuple{2,Integer}) where T
@@ -306,7 +308,7 @@ for Typ in (:ConstRows, :PertConstRows)
     end
 end
 
-
+const TridiagonalToeplitzLayout = Union{SymTridiagonalLayout{FillLayout},TridiagonalLayout{FillLayout}}
 const BandedToeplitzLayout = BandedColumns{ConstRows}
 const PertToeplitzLayout = BandedColumns{PertConstRows}
 const PertTriangularToeplitzLayout{UPLO,UNIT} = TriangularLayout{UPLO,UNIT,BandedColumns{PertConstRows}}
@@ -383,9 +385,42 @@ function _bandedfill_mul(M::MulAdd, ::Tuple{InfAxes,InfAxes}, ::Tuple{InfAxes,In
     ret
 end
 
-mulapplystyle(::BandedColumns{FillLayout}, ::PertToeplitzLayout) = LazyArrayApplyStyle()
-mulapplystyle(::PertToeplitzLayout, ::BandedColumns{FillLayout}) = LazyArrayApplyStyle()
-mulapplystyle(::BandedColumns{FillLayout}, ::BandedToeplitzLayout) = LazyArrayApplyStyle()
-mulapplystyle(::BandedToeplitzLayout, ::BandedColumns{FillLayout}) = LazyArrayApplyStyle()
-mulapplystyle(::AbstractQLayout, ::BandedToeplitzLayout) = LazyArrayApplyStyle()
-mulapplystyle(::AbstractQLayout, ::PertToeplitzLayout) = LazyArrayApplyStyle()
+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)
+
+
+function _bidiag_forwardsub!(M::Ldiv{<:Any,<:PaddedLayout})
+    A, b_in = M.A, M.B
+    dv = diagonaldata(A)
+    ev = subdiagonaldata(A)
+    b = paddeddata(b_in)
+    N = length(b)
+    b[1] = bj1 = dv[1]\b[1]
+    @inbounds for j = 2:N
+        bj  = b[j]
+        bj -= ev[j - 1] * bj1
+        dvj = dv[j]
+        if iszero(dvj)
+            throw(SingularEbception(j))
+        end
+        bj   = dvj\bj
+        b[j] = bj1 = bj
+    end
+
+    @inbounds for j = N+1:length(b_in)
+        iszero(bj1) && break
+        bj = -ev[j - 1] * bj1
+        dvj = dv[j]
+        if iszero(dvj)
+            throw(SingularEbception(j))
+        end
+        bj   = dvj\bj
+        b_in[j] = bj1 = bj
+    end
+
+    b_in
+end

src/blockbanded/blockbanded.jl

@@ -1,18 +1,6 @@
 sizes_from_blocks(A::AbstractVector, ::Tuple{OneToInf{Int}}) = (map(length,A),)
 
 
-
-MemoryLayout(::Type{<:PseudoBlockArray{<:Any,N,<:Arr}}) where {N,Arr} = 
-    MemoryLayout(Arr)
-
-struct BlockLayout{LAY} <: MemoryLayout end
-
-MemoryLayout(::Type{<:BlockArray{<:Any,N,<:Arr}}) where {N,Arr} = 
-    blocklayout(MemoryLayout(Arr))
-
-blocklayout(_) = BlockLayout{UnknownLayout}()    
-blocklayout(::LazyLayout) = BlockLayout{LazyLayout}()    
-
 # for LazyLay in (:(BlockLayout{LazyLayout}), :(TriangularLayout{UPLO,UNIT,BlockLayout{LazyLayout}} where {UPLO,UNIT}))
 #     @eval begin
 #         combine_mul_styles(::$LazyLay) = LazyArrayApplyStyle()

src/blockbanded/infblocktridiagonal.jl

@@ -1,26 +1,27 @@
 const BlockTriPertToeplitz{T} = BlockMatrix{T,Tridiagonal{Matrix{T},Vcat{Matrix{T},1,Tuple{Vector{Matrix{T}},Fill{Matrix{T},1,Tuple{OneToInf{Int}}}}}},
                                         NTuple{2,BlockedUnitRange{Vcat{Int,1,Tuple{Vector{Int},InfStepRange{Int,Int}}}}}}
 
+const BlockTridiagonalToeplitzLayout = BlockLayout{TridiagonalLayout{FillLayout},DenseColumnMajor}
 
 function BlockTridiagonal(adjA::Adjoint{T,BlockTriPertToeplitz{T}}) where T
     A = parent(adjA)
-    BlockTridiagonal(Matrix.(adjoint.(A.blocks.du)), 
-                     Matrix.(adjoint.(A.blocks.d)), 
+    BlockTridiagonal(Matrix.(adjoint.(A.blocks.du)),
+                     Matrix.(adjoint.(A.blocks.d)),
                      Matrix.(adjoint.(A.blocks.dl)))
-end                                       
+end
 
 for op in (:-, :+)
     @eval begin
-        function $op(A::BlockTriPertToeplitz{T}, λ::UniformScaling) where T 
+        function $op(A::BlockTriPertToeplitz{T}, λ::UniformScaling) where T
             TV = promote_type(T,eltype(λ))
-            BlockTridiagonal(Vcat(convert.(AbstractVector{Matrix{TV}}, A.blocks.dl.args)...), 
-                             Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.d, Ref(λ)).args)...), 
+            BlockTridiagonal(Vcat(convert.(AbstractVector{Matrix{TV}}, A.blocks.dl.args)...),
+                             Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.d, Ref(λ)).args)...),
                              Vcat(convert.(AbstractVector{Matrix{TV}}, A.blocks.du.args)...))
         end
         function $op(λ::UniformScaling, A::BlockTriPertToeplitz{V}) where V
             TV = promote_type(eltype(λ),V)
-            BlockTridiagonal(Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.dl.args))...), 
-                             Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, Ref(λ), A.blocks.d).args)...), 
+            BlockTridiagonal(Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.dl.args))...),
+                             Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, Ref(λ), A.blocks.d).args)...),
                              Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.du.args))...))
         end
         $op(adjA::Adjoint{T,BlockTriPertToeplitz{T}}, λ::UniformScaling) where T = $op(BlockTridiagonal(adjA), λ)
@@ -31,8 +32,8 @@ end
 *(a::AbstractVector, b::AbstractFill{<:Any,2,Tuple{OneTo{Int},OneToInf{Int}}}) = ApplyArray(*,a,b)
 
 sizes_from_blocks(A::Diagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.diag, 1), size.(A.diag,2)
-
 sizes_from_blocks(A::Tridiagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.d, 1), size.(A.d,2)
+sizes_from_blocks(A::Bidiagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.dv, 1), size.(A.dv,2)
 
 axes_print_matrix_row(_, io, X, A, i, ::AbstractVector{<:Infinity}, sep) = nothing
 axes_print_matrix_row(::NTuple{2,BlockedUnitRange}, io, X, A, i, ::AbstractVector{<:Infinity}, sep) = nothing
@@ -46,7 +47,7 @@ function BlockSkylineSizes(A::BlockTriPertToeplitz, (l,u)::NTuple{2,Int})
     for J=1:N
         block_starts[max(1,J-u),J] = J == 1 ? 1 :
                             block_starts[max(1,J-1-u),J-1]+block_sizes[J-1]*block_strides[J-1]
-                                
+
         for K=max(1,J-u)+1:J+l
             block_starts[K,J] = block_starts[K-1,J]+size(A[Block(K-1,J)],1)
         end
@@ -72,15 +73,15 @@ function BlockSkylineSizes(A::BlockTriPertToeplitz, (l,u)::NTuple{2,Int})
 end
 
 function BlockBandedMatrix(A::BlockTriPertToeplitz{T}, (l,u)::NTuple{2,Int}) where T
-    data = T[]                      
+    data = T[]
     append!(data,vec(A[Block.(1:1+l),Block(1)]))
     N = max(length(A.blocks.du.args[1])+1,length(A.blocks.d.args[1]),length(A.blocks.dl.args[1]))
     for J=2:N
         append!(data, vec(A[Block.(max(1,J-u):J+l),Block(J)]))
     end
     tl = mortar(Fill(vec(A[Block.(max(1,N+1-u):N+1+l),Block(N+1)]),∞))
-    
+
     B = _BlockSkylineMatrix(Vcat(data,tl), BlockSkylineSizes(A, (l,u)))
-end    
+end
 
 BlockBandedMatrix(A::BlockTriPertToeplitz) = BlockBandedMatrix(A, blockbandwidths(A))

src/infql.jl

@@ -280,6 +280,7 @@ function (*)(A::Adjoint{T,<:QLPackedQ{T,<:InfBlockBandedMatrix}}, x::AbstractVec
 end
 
 ldiv!(F::QLProduct, b::AbstractVector) = ldiv!(F.L, lmul!(F.Q',b))
+ldiv!(F::QLProduct, b::LayoutVector) = ldiv!(F.L, lmul!(F.Q',b))
 
 function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N',BandedColumns{PertConstRows}},<:PaddedLayout})
     A,b = M.A,M.B

src/infqr.jl

@@ -19,7 +19,7 @@ function AdaptiveQRData(::AbstractAlmostBandedLayout, A::AbstractMatrix{T}) wher
     AdaptiveQRData(CachedArray(data,A,(0,0)), Vector{T}(), 0)
 end
 
-function AdaptiveQRData(::AbstractBlockBandedLayout, A::AbstractMatrix{T}) where T
+function AdaptiveQRData(::AbstractBlockLayout, A::AbstractMatrix{T}) where T
     l,u = blockbandwidths(A)
     m,n = axes(A)
     data = BlockBandedMatrix{T}(undef,(m[Block.(1:2l+u+1)],n[Block.(1:0)]),(l,l+u)) # pad super
@@ -149,6 +149,7 @@ end
 
 _qr(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A)
 _qr(::AbstractAlmostBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A)
+__qr(_, ::NTuple{2,Infinity}, A) = adaptiveqr(A)
 _qr(::AbstractBlockBandedLayout, ::NTuple{2,Infinity}, A) = adaptiveqr(A)
 
 partialqr!(F::QR, n) = partialqr!(F.factors, n)
@@ -306,7 +307,9 @@ end
 ldiv!(dest::AbstractVector, F::QR{<:Any,<:AdaptiveQRFactors}, b::AbstractVector) =
     ldiv!(F, copyto!(dest, b))
 ldiv!(F::QR{<:Any,<:AdaptiveQRFactors}, b::AbstractVector) = ldiv!(F.R, lmul!(F.Q',b))
+ldiv!(F::QR{<:Any,<:AdaptiveQRFactors}, b::LayoutVector) = ldiv!(F.R, lmul!(F.Q',b))
 \(F::QR{<:Any,<:AdaptiveQRFactors}, B::AbstractVector) = ldiv!(F.R, F.Q'B)
+\(F::QR{<:Any,<:AdaptiveQRFactors}, B::LayoutVector) = ldiv!(F.R, F.Q'B)
 
 
 factorize(A::BandedMatrix{<:Any,<:Any,<:OneToInf}) = qr(A)

src/inful.jl

@@ -0,0 +1,75 @@
+"""
+_ultailL1(C, A, B)
+
+gives L[Block(1,1)] of a block-tridiagonal Toeplitz operator. Based on Delvaux and Dette 2012.
+"""
+
+# need to solve [1 -B*inv(L); 0 1] * [C A B; 0 C L] == [C L 0; 0 C L]
+# that is 
+# A - B*inv(L)*C == L
+# In the scalar case this is quadratic
+# L^2 - A*L + B*C == 0
+# so that 
+# = L = (A ± sqrt(A^2 - 4B*C))/2
+# we choose sign(A) to maximise the magnitude as we know inv(T)[1,1] -> 0, hence
+# inv(L) -> 0
+_ultailL1(c::Number, a::Number, b::Number) = (a + sign(a)*sqrt(a^2-4b*c))/2
+
+# In the matrix case we write inv(L)*C = R (double check) to make it quadratic
+# A - B*R == inv(C)inv(R)
+# C*A*R - C*B*R^2 == I
+# and do eigen decomposition of R to reduce this to a companion matrix problem.
+
+function _ultailL1(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
+    d = size(A,1)
+    λs = filter!(λ -> abs(λ) ≤ 1, eigvals([zeros(2,2) -B; -C -A], [B zeros(2,2); zeros(2,2) -B]))
+    @assert length(λs) == d
+    V = Matrix{eltype(λs)}(undef, d, d)
+    for (j,λ) in enumerate(λs)
+        V[:,j] = svd(A-λ*B - C/λ).V[:,end] # nullspace(A-λ*B - C/λ)
+    end
+    C*(V*Diagonal(inv.(λs))/V)
+end
+
+function _ul(::TridiagonalToeplitzLayout, J::AbstractMatrix, ::Val{false}; check::Bool = true)
+    C = getindex_value(subdiagonaldata(J))
+    A = getindex_value(diagonaldata(J))
+    B = getindex_value(supdiagonaldata(J))
+    L = _ultailL1(C, A, B)
+    U = B/L
+    UL(Tridiagonal(Fill(convert(typeof(L),C),∞), Fill(L,∞), Fill(U,∞)), OneToInf(), 0)
+end
+
+function _ul(::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)
+    C = getindex_value(subdiagonaldata(blocks(J)))
+    A = getindex_value(diagonaldata(blocks(J)))
+    B = getindex_value(supdiagonaldata(blocks(J)))
+    L = _ultailL1(C, A, B)
+    U = B/L
+    UL(mortar(Tridiagonal(Fill(convert(typeof(L),C),∞), Fill(L,∞), Fill(U,∞))), OneToInf(), 0)
+end
+
+
+_inf_getU(::TridiagonalToeplitzLayout, F::UL) = Bidiagonal(Fill(one(eltype(F)),∞),F.factors.du, :U)
+_inf_getL(::TridiagonalToeplitzLayout, F::UL) = Bidiagonal(F.factors.d,F.factors.dl, :L)
+
+
+function _inf_getU(::BlockTridiagonalToeplitzLayout, F::UL)
+    d = size(F.factors.blocks.du[1],1)
+    II = convert(eltype(F.factors.blocks.du), I(d))
+    mortar(Bidiagonal(Fill(II,∞),F.factors.blocks.du, :U))
+end
+
+_inf_getL(::BlockTridiagonalToeplitzLayout, F::UL) = mortar(Bidiagonal(F.factors.blocks.d,F.factors.blocks.dl, :L))
+
+
+getU(F::UL, ::NTuple{2,Infinity}) = _inf_getU(MemoryLayout(F.factors), F)
+getL(F::UL, ::NTuple{2,Infinity}) = _inf_getL(MemoryLayout(F.factors), F)