Implement finite section based adaptive QL (#125)

* basic implementation of finite section based QL

* import AbstractCachedMatrix

* bug fixes + basic tests work

* more tests and lowertriangular layout

* simplify loops somewhat

* stop checking full finite sections + refactor

* refactor before working on matching QL structure

* change underlying structs

* work towards expanding tau and factors simultaneously

* Update infql.jl

* explicitly import some lazyarrays functions + type declarations

* big style overhaul

* Update src/infql.jl

Co-authored-by: Sheehan Olver <solver@mac.com>

* tidying, test fix, todo comment

* move getindex from here to LazyArrays

* Update Project.toml

* change Q implementation to be lazy

* import the default getQ

* bug fix and tidy lazy implementation

* fix adjoint

* reinstatate some changes

* add tests, remove a function

* add a test for complex entries

* add test for ldiv behavior of Q

* LazyVector -> LayoutVector

* nzzeros -> last(colsupport

* restore non-tridiagonal functionality

* Update test_infql.jl

* Add matrix tests, remove unnecessary * overrides

* lazy multiplication by Q

* Update Project.toml

---------

Co-authored-by: Sheehan Olver <solver@mac.com>

Author: TSGut

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.6.17"
+version = "0.6.18"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -25,7 +25,7 @@ BlockBandedMatrices = "0.11.5, 0.12"
 DSP = "0.7"
 FillArrays = "0.13, 1"
 InfiniteArrays = "0.12"
-LazyArrays = "0.22, 1"
+LazyArrays = "1.1.1"
 LazyBandedMatrices = "0.8.7"
 MatrixFactorizations = "0.9.6, 1"
 SemiseparableMatrices = "0.3"

src/InfiniteLinearAlgebra.jl

@@ -12,17 +12,17 @@ import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted
 import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, TridiagonalLayout,
                         sublayout, _qr, __qr, MatLmulVec, MatLmulMat, AbstractQLayout, materialize!, diagonaldata, subdiagonaldata, supdiagonaldata,
                         _bidiag_forwardsub!, mulreduce, RangeCumsum, _factorize, transposelayout, ldiv!, lmul!, mul, CNoPivot
-import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns,
+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 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,
                     factorize, sub_materialize, LazyLayout, LazyArrayStyle, layout_getindex,
-                    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,
+                    applylayout, ApplyLayout, PaddedLayout, CachedLayout, AbstractCachedVector, AbstractCachedMatrix, cacheddata, zero!, MulAddStyle, ApplyArray,
+                    LazyArray, LazyMatrix, LazyVector, paddeddata, arguments, resizedata!
+import MatrixFactorizations: ul, ul!, _ul, ql, ql!, _ql, QLPackedQ, getL, getR, getQ, getU, reflector!, reflectorApply!, QL, QR, QRPackedQ,
                             QRPackedQLayout, AdjQRPackedQLayout, QLPackedQLayout, AdjQLPackedQLayout, LayoutQ
 
 import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange, blockcolsupport, BlockLayout, AbstractBlockLayout, BlockSlice

src/banded/hessenbergq.jl

@@ -112,7 +112,7 @@ function materialize!(L::MatLmulVec{<:HessenbergQLayout{'L'}})
     Q, x = L.A,L.B
     T = eltype(Q)
     t = Array{T}(undef, 2)
-    nz = nzzeros(x,1)
+    nz = last(colsupport(x))
     for n = 1:length(Q.q)
         v = view(x, n:n+1)
         mul!(t, Q.q[n], v)
@@ -126,7 +126,7 @@ function materialize!(L::MatLmulVec{<:HessenbergQLayout{'U'}})
     Q, x = L.A,L.B
     T = eltype(Q)
     t = Array{T}(undef, 2)
-    for n = min(length(Q.q),nzzeros(x,1)):-1:1
+    for n = min(length(Q.q),last(colsupport(x))):-1:1
         v = view(x, n:n+1)
         mul!(t, Q.q[n], v)
         v .= t

src/infql.jl

@@ -92,23 +92,9 @@ function ql_hessenberg!(B::InfBandedMatrix{TT}; kwds...) where TT
     QLHessenberg(_BandedMatrix(H, ℵ₀, l, 1), Vcat( LowerHessenbergQ(F.Q).q, F∞.q))
 end
 
-getindex(Q::QLPackedQ{T,<:InfBandedMatrix{T}}, i::Int, j::Int) where T =
-    (Q'*[Zeros{T}(i-1); one(T); Zeros{T}(∞)])[j]'
-getindex(Q::QLPackedQ{<:Any,<:InfBandedMatrix}, I::AbstractVector{Int}, J::AbstractVector{Int}) =
-    [Q[i,j] for i in I, j in J]
-
 getL(Q::QL, ::NTuple{2,InfiniteCardinal{0}}) = LowerTriangular(Q.factors)
 getL(Q::QLHessenberg, ::NTuple{2,InfiniteCardinal{0}}) = LowerTriangular(Q.factors)
 
-# number of structural non-zeros in axis k
-nzzeros(A::AbstractArray, k) = size(A,k)
-nzzeros(::Zeros, k) = 0
-nzzeros(B::Vcat, k) = sum(size.(B.args[1:end-1],k))
-nzzeros(B::CachedArray, k) = max(B.datasize[k], nzzeros(B.array,k))
-function nzzeros(B::AbstractMatrix, k)
-    l,u = bandwidths(B)
-    k == 1 ? size(B,2) + l : size(B,1) + u
-end
 
 function materialize!(M::Lmul{<:QLPackedQLayout{<:BandedColumns},<:PaddedLayout})
     A,B = M.A,M.B
@@ -123,7 +109,7 @@ function materialize!(M::Lmul{<:QLPackedQLayout{<:BandedColumns},<:PaddedLayout}
     D = Afactors.data
     for k = 1:∞
         ν = k
-        allzero = k > nzzeros(B,1) ? true : false
+        allzero = k > last(colsupport(B)) ? true : false
         for j = 1:nB
             vBj = B[k,j]
             for i = max(1,ν-u):k-1
@@ -157,7 +143,7 @@ function materialize!(M::Lmul{<:AdjQLPackedQLayout{<:BandedColumns},<:PaddedLayo
     l,u = bandwidths(Afactors)
     D = Afactors.data
     @inbounds begin
-        for k = nzzeros(B,1)+u:-1:1
+        for k = last(colsupport(B))+u:-1:1
             ν = k
             for j = 1:nB
                 vBj = B[k,j]
@@ -175,15 +161,6 @@ function materialize!(M::Lmul{<:AdjQLPackedQLayout{<:BandedColumns},<:PaddedLayo
     B
 end
 
-function _lmul_cache(A::Union{AbstractMatrix{T},AbstractQ{T}}, x::AbstractVector{S}) where {T,S}
-    TS = promote_op(matprod, T, S)
-    lmul!(A, cache(convert(AbstractVector{TS},x)))
-end
-
-(*)(A::QLPackedQ{T,<:InfBandedMatrix}, x::AbstractVector) where {T} = _lmul_cache(A, x)
-(*)(A::AdjointQtype{T,<:QLPackedQ{T,<:InfBandedMatrix}}, x::AbstractVector) where {T} = _lmul_cache(A, x)
-(*)(A::QLPackedQ{T,<:InfBandedMatrix}, x::LazyVector) where {T} = _lmul_cache(A, x)
-(*)(A::AdjointQtype{T,<:QLPackedQ{T,<:InfBandedMatrix}}, x::LazyVector) where {T} = _lmul_cache(A, x)
 
 
 function blocktailiterate(c,a,b, d=c, e=a)
@@ -246,7 +223,7 @@ function lmul!(adjA::AdjointQtype{<:Any,<:QLPackedQ{<:Any,<:InfBlockBandedMatrix
     l = 2l+1
     u = 2u+1
     @inbounds begin
-        for k = nzzeros(B,1)+u:-1:1
+        for k = last(colsupport(B))+u:-1:1
             ν = k
             for j = 1:nB
                 vBj = B[k,j]
@@ -290,7 +267,7 @@ function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N',BandedColumns{Per
         throw(DimensionMismatch("second dimension of left hand side A, $n, and length of right hand side b, $(length(b)), must be equal"))
     end
     data = triangulardata(A)
-    nz = nzzeros(b,1)
+    nz = last(colsupport(b))
     @inbounds for j in 1:n
         iszero(data[j,j]) && throw(SingularException(j))
         bj = b[j] = data[j,j] \ b[j]
@@ -333,8 +310,6 @@ function _ql(::SymTridiagonalLayout, ::NTuple{2,OneToInf{Int}}, A, args...; kwds
     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...)
@@ -378,3 +353,169 @@ function LazyBandedMatrices._SymTridiagonal(::Tuple{TriangularLayout{'L', 'N', P
     ev = [A[k,k+1] for k=1:m-1]
     SymTridiagonal([dv; Fill(dv[end],∞)], [ev; Fill(ev[end],∞)])
 end
+
+###
+# Experimental adaptive finite section QL
+###
+mutable struct AdaptiveQLTau{T} <: AbstractCachedVector{T}
+    data::Vector{T}
+    M::AbstractMatrix{T}
+    datasize::Integer
+    tol::Real
+    AdaptiveQLTau{T}(D, M, N::Integer, tol) where T = new{T}(D, M, N, tol)
+end
+mutable struct AdaptiveQLFactors{T} <: AbstractCachedMatrix{T}
+    data::BandedMatrix{T}
+    M::AbstractMatrix{T}
+    datasize::Tuple{Int, Int}
+    tol::Real
+    AdaptiveQLFactors{T}(D, M, N::Tuple{Int, Int}, tol) where T = new{T}(D, M, N, tol)
+end
+
+size(::AdaptiveQLFactors) = (ℵ₀, ℵ₀)
+size(::AdaptiveQLTau) = (ℵ₀, )
+
+# adaptive QL accepts optional tolerance
+function ql(A::InfBandedMatrix{T}, tol = eps(float(real(T)))) where T
+    factors, τ = initialadaptiveQLblock(A, tol)
+    QL(AdaptiveQLFactors{T}(factors, A, size(factors), tol),AdaptiveQLTau{T}(τ, A, length(τ), tol))
+end
+
+# Computes the initial data for the finite section based QL decomposition
+function initialadaptiveQLblock(A::AbstractMatrix{T}, tol) where T
+    maxN = 10000   # Prevent runaway loop
+    j = 50         # We initialize with a 50 × 50 block that is adaptively expanded
+    Lerr = one(real(T))
+    N = j
+    checkinds = max(1,j-bandwidth(A,1)-bandwidth(A,2))
+    @inbounds Ls = ql(A[checkinds:N,checkinds:N]).L[2:j-checkinds+1,2:j-checkinds+1]
+    @inbounds while Lerr > tol
+        # compute QL for small finite section and large finite section
+        Ll = ql(A[checkinds:2N,checkinds:2N]).L[2:j-checkinds+1,2:j-checkinds+1]
+        # compare bottom right sections and stop if desired level of convergence achieved
+        Lerr = norm(Ll-Ls,2)
+        if N == maxN
+            error("Reached max. iterations in adaptive QL without convergence to desired tolerance.")
+        end
+        Ls = Ll
+        N = 2*N
+    end
+    F = ql(A[1:(N÷2),1:(N÷2)])
+    return (F.factors[1:50,1:50], F.τ[1:50])
+end
+
+# Resize and filling functions for cached implementation
+function resizedata!(K::AdaptiveQLFactors, nm...)
+    nm = maximum(nm)
+    νμ = K.datasize[1]
+    if nm > νμ
+        olddata = copy(K.data)
+        K.data = similar(K.data, nm, nm)
+        K.data[axes(olddata)...] = olddata
+        inds = νμ:nm
+        cache_filldata!(K, inds)
+        K.datasize = size(K.data)
+    end
+    K
+end
+
+function resizedata!(K::AdaptiveQLTau, nm...)
+    nm = maximum(nm)
+    νμ = K.datasize
+    if nm > νμ
+        resize!(K.data,nm)
+        cache_filldata!(K, νμ:nm)
+        K.datasize = size(K.data,1)
+    end
+    K
+end
+
+function cache_filldata!(A::AdaptiveQLFactors{T}, inds::UnitRange{Int}) where T
+    j = maximum(inds)
+    maxN = 1000*j # Prevent runaway loop
+    Lerr = one(real(T))
+    N = j
+    checkinds = max(1,j-bandwidth(A.M,1)-bandwidth(A.M,2))
+    @inbounds Ls = ql(A.M[checkinds:N,checkinds:N]).L[2:j-checkinds+1,2:j-checkinds+1]
+    @inbounds while Lerr > A.tol
+        # compute QL for small finite section and large finite section
+        Ll = ql(A.M[checkinds:2N,checkinds:2N]).L[2:j-checkinds+1,2:j-checkinds+1]
+        # compare bottom right sections and stop if desired level of convergence achieved
+        Lerr = norm(Ll-Ls,2)
+        if N == maxN
+            error("Reached max. iterations in adaptive QL without convergence to desired tolerance.")
+        end
+        Ls = Ll
+        N = 2*N
+    end
+    A.data = ql(A.M[1:(N÷2),1:(N÷2)]).factors[1:j,1:j]
+end
+
+function cache_filldata!(A::AdaptiveQLTau{T}, inds::UnitRange{Int}) where T
+    j = maximum(inds)
+    maxN = 1000*j
+    Lerr = one(real(T))
+    N = j
+    checkinds = max(1,j-bandwidth(A.M,1)-bandwidth(A.M,2))
+    @inbounds Ls = ql(A.M[checkinds:N,checkinds:N]).L[2:j-checkinds+1,2:j-checkinds+1]
+    @inbounds while Lerr > A.tol
+        # compute QL for small finite section and large finite section
+        Ll = ql(A.M[checkinds:2N,checkinds:2N]).L[2:j-checkinds+1,2:j-checkinds+1]
+        # compare bottom right sections and stop if desired level of convergence achieved
+        Lerr = norm(Ll-Ls,2)
+        if N == maxN
+            error("Reached max. iterations in adaptive QL without convergence to desired tolerance.")
+        end
+        Ls = Ll
+        N = 2*N
+    end
+    A.data = ql(A.M[1:(N÷2),1:(N÷2)]).τ[1:j]
+end
+
+# TODO: adaptively build L*b using caching and forward-substitution
+*(L::LowerTriangular{T, AdaptiveQLFactors{T}}, b::LayoutVector) where T = ApplyArray(*, L, b)
+
+MemoryLayout(::AdaptiveQLFactors) = LazyBandedLayout()
+bandwidths(F::AdaptiveQLFactors) = bandwidths(F.data)
+
+# Q = \\prod_{i=1}^{\\min(m,n)} (I - \\tau_i v_i v_i^T)
+getindex(Q::QLPackedQ{T,<:AdaptiveQLFactors{T}}, i::Int, j::Int) where T =
+(Q'*[Zeros{T}(i-1); one(T); Zeros{T}(∞)])[j]'
+getindex(Q::QLPackedQ{<:Any,<:AdaptiveQLFactors}, I::AbstractVector{Int}, J::AbstractVector{Int}) =
+    [Q[i,j] for i in I, j in J]
+getindex(Q::QLPackedQ{<:Any,<:AdaptiveQLFactors}, I::Int, J::UnitRange{Int}) =
+    [Q[i,j] for i in I, j in J]
+getindex(Q::QLPackedQ{<:Any,<:AdaptiveQLFactors}, I::UnitRange{Int}, J::Int) =
+    [Q[i,j] for i in I, j in J]
+
+materialize!(M::Lmul{<:QLPackedQLayout{<:LazyLayout},<:PaddedLayout}) = ApplyArray(*,M.A,M.B)
+
+function materialize!(M::Lmul{<:AdjQLPackedQLayout{<:LazyLayout},<:PaddedLayout})
+    adjA,B = M.A,M.B
+    A = parent(adjA)
+    mA, nA = size(A.factors)
+    mB, nB = size(B,1), size(B,2)
+    if mA != mB
+        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
+    end
+    Afactors = A.factors
+    l,u = bandwidths(Afactors)
+    l = 2l+1
+    u = 2u+1
+    @inbounds begin
+        for k = last(colsupport(B))+u:-1:1
+            for j = 1:nB
+                vBj = B[k,j]
+                for i = max(1,k-u):k-1
+                    vBj += conj(Afactors[i,k])*B[i,j]
+                end
+                vBj = conj(A.τ[k])*vBj
+                B[k,j] -= vBj
+                for i = max(1,k-u):k-1
+                    B[i,j] -= Afactors[i,k]*vBj
+                end
+            end
+        end
+    end
+    B
+end