Almost banded adaptive QR (#22)

* Add Semiseparable

* Update for ArrayLayouts 0.2

* Start almost banded adaptive QR

* _qr

* update packages

* Update Project.toml

* Update Project.toml

* debug frozen test

* Revert "debug frozen test"

This reverts commit dc319ff.

* Revert "Revert "debug frozen test""

This reverts commit 183c0da.

Author: dlfivefifty

Project.toml

@@ -13,17 +13,19 @@ LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
+SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 
 [compat]
-ArrayLayouts = "0.1"
-BandedMatrices = "0.14"
-BlockArrays = "0.11"
-BlockBandedMatrices = "0.7.1"
-FillArrays = "0.8.4"
-InfiniteArrays = "0.6.1"
-LazyArrays = "0.15"
-LazyBandedMatrices = "0.2"
-MatrixFactorizations = "0.2"
+ArrayLayouts = "0.2.1"
+BandedMatrices = "0.15.1"
+BlockArrays = "0.12"
+BlockBandedMatrices = "0.8"
+FillArrays = "0.8.6"
+InfiniteArrays = "0.7"
+LazyArrays = "0.16.4"
+LazyBandedMatrices = "0.2.5"
+MatrixFactorizations = "0.3.1"
+SemiseparableMatrices = "0.0.1"
 julia = "1.3"
 
 [extras]

src/InfiniteLinearAlgebra.jl

@@ -1,24 +1,24 @@
 module InfiniteLinearAlgebra
-using BlockArrays, BlockBandedMatrices, BandedMatrices, LazyArrays, LazyBandedMatrices, 
+using BlockArrays, BlockBandedMatrices, BandedMatrices, LazyArrays, LazyBandedMatrices, SemiseparableMatrices,
         FillArrays, InfiniteArrays, MatrixFactorizations, ArrayLayouts, LinearAlgebra
 
 import Base: +, -, *, /, \, ^, OneTo, getindex, promote_op, _unsafe_getindex, print_matrix_row, size, axes,
             AbstractMatrix, AbstractArray, Matrix, Array, Vector, AbstractVector, Slice,
             show, getproperty, copy, map
 import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted
 
-import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, sublayout
+import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, sublayout, _qr
 import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns,
                         _default_banded_broadcast
 import FillArrays: AbstractFill, getindex_value      
 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,
-                    factorize, sub_materialize, LazyLayout, LazyArrayStyle, lazy_getindex,
+                    factorize, sub_materialize, LazyLayout, LazyArrayStyle, layout_getindex,
                     @lazymul, applylayout, ApplyLayout, PaddedLayout, materialize!, zero!, MulAddStyle,
                     LazyArray, LazyMatrix, LazyVector
-import MatrixFactorizations: ql, ql!, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ
+import MatrixFactorizations: ql, ql!, _ql, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ
 
 import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange
 
@@ -29,6 +29,7 @@ import LazyBandedMatrices: MulBandedLayout, BroadcastBandedLayout
 import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, _BlockSkylineMatrix, blockstart, blockstride,
         BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal
 
+import SemiseparableMatrices: AbstractAlmostBandedLayout, _almostbanded_qr!
 
 LazyArrays.@lazymul BandedMatrix{<:Any,<:Any,<:OneToInf}
 *(A::BandedMatrix{<:Any,<:Any,<:OneToInf}, b::CachedVector) = apply(*,A,b)

src/banded/infbanded.jl

@@ -2,7 +2,7 @@
 # Diagonal
 ###
 
-getindex(D::Diagonal, k::InfAxes, j::InfAxes) = lazy_getindex(D, k, j)
+getindex(D::Diagonal, k::InfAxes, j::InfAxes) = layout_getindex(D, k, j)
 
 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}}}}}

src/infql.jl

@@ -208,7 +208,7 @@ end
 # BlockTridiagonal
 ####
 
-function _ql(A::BlockTriPertToeplitz, d, e)
+function _blocktripert_ql(A, d, e)
     N = max(length(A.blocks.du.args[1])+1,length(A.blocks.d.args[1]),length(A.blocks.dl.args[1]))
     c,a,b = A[Block(N+1,N)],A[Block(N,N)],A[Block(N-1,N)]
     P,τ = blocktailiterate(c,a,b,d,e)
@@ -226,7 +226,7 @@ function _ql(A::BlockTriPertToeplitz, d, e)
             Vcat(F.τ,mortar(Fill(τ,∞)))), P[Block(1,1)], P[Block(1,2)]
 end
 
-ql(A::BlockTriPertToeplitz) = _ql(A, A[Block(2,3)], A[Block(3,3)])[1]
+ql(A::BlockTriPertToeplitz) = _blocktripert_ql(A, A[Block(2,3)], A[Block(3,3)])[1]
 
 ql(A::Adjoint{T,BlockTriPertToeplitz{T}}) where T = ql(BlockTridiagonal(A))
 

src/infqr.jl

@@ -10,24 +10,47 @@ function AdaptiveQRData(::AbstractBandedLayout, A::AbstractMatrix{T}) where T
     data = BandedMatrix{T}(undef,(2l+u+1,0),(l,l+u)) # pad super
     AdaptiveQRData(CachedArray(data,A), Vector{T}(), 0)
 end
+
+function AdaptiveQRData(::AbstractAlmostBandedLayout, A::AbstractMatrix{T}) where T 
+    l,u = almostbandwidths(A)
+    r = almostbandedrank(A)
+    data = AlmostBandedMatrix{T}(undef,(2l+u+1,0),(l,l+u),r) # pad super
+    
+    AdaptiveQRData(CachedArray(data,A,(0,0)), Vector{T}(), 0)
+end
+
 AdaptiveQRData(A::AbstractMatrix{T}) where T = AdaptiveQRData(MemoryLayout(typeof(A)), A)
 
 function partialqr!(F::AdaptiveQRData{<:Any,<:BandedMatrix}, n::Int)
     if n > F.ncols 
-        l,u = bandwidths(F.data.array)
-        resizedata!(F.data,n+l,n+u+l);
+        l,u = bandwidths(F.data.data)
+        resizedata!(F.data,n+l,n+u);
         resize!(F.τ,n);
         ñ = F.ncols
         τ = view(F.τ,ñ+1:n);
         if l ≤ 0 
             zero!(τ)
         else
-            factors = view(F.data.data,ñ+1:n+l,ñ+1:n);
-            _banded_qr!(factors, τ);
-            # multiply remaining columns
-            n̄ = max(ñ+1,n-l-u+1) # first column interacting with extra data
-            Q = QRPackedQ(view(F.data.data,n̄:n+l,n̄:n), view(F.τ,n̄:n))
-            lmul!(Q',view(F.data.data,n̄:n+l,n+1:n+u+l))
+            factors = view(F.data.data,ñ+1:n+l,ñ+1:n+u);
+            _banded_qr!(factors, τ, n-ñ)
+        end
+        F.ncols = n
+    end
+    F
+end
+
+function partialqr!(F::AdaptiveQRData{<:Any,<:AlmostBandedMatrix}, n::Int)
+    if n > F.ncols 
+        l,u = almostbandwidths(F.data.data)
+        resizedata!(F.data,n+l,n+l+u);
+        resize!(F.τ,n);
+        ñ = F.ncols
+        τ = view(F.τ,ñ+1:n)
+        if l ≤ 0 
+            zero!(τ)
+        else
+            factors = view(F.data.data,ñ+1:n+l,ñ+1:n+l+u)
+            _almostbanded_qr!(factors, τ, n-ñ)
         end
         F.ncols = n
     end
@@ -82,11 +105,17 @@ end
 getR(Q::QR, ::Tuple{OneToInf{Int},OneToInf{Int}}) where T = UpperTriangular(Q.factors)
 
 
-function _banded_qr(::NTuple{2,OneToInf{Int}}, A)
+function adaptiveqr(A)
     data = AdaptiveQRData(A)
     QR(AdaptiveQRFactors(data), AdaptiveQRTau(data))
 end
 
+_qr(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A)
+_qr(::AbstractAlmostBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A)
+
+partialqr!(F::QR, n) = partialqr!(F.factors, n)
+partialqr!(F::AdaptiveQRFactors, n) = partialqr!(F.data, n)
+
 #########
 # lmul!
 #########
@@ -142,11 +171,12 @@ function lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AdaptiveQRFactors}}, B::C
 
     @inbounds begin
         while first(jr) < nA
+            j = first(jr)
             cs = colsupport(A.factors, last(jr))
             cs_max = maximum(cs)
-            kr = first(jr):cs_max
+            kr = j:cs_max
             resizedata!(B, min(cs_max,mB))
-            if (first(jr) > sB && maximum(abs,view(B.data,colsupport(A.factors,first(jr)))) ≤ tol)
+            if (j > sB && maximum(abs,view(B.data,j:last(colsupport(A.factors,j)))) ≤ tol)
                 break
             end
             partialqr!(A.factors.data, min(cs_max,nA))

test/runtests.jl

@@ -60,7 +60,7 @@ end
 end
 
 
-@testset "Algebra" begin 
+# @testset "Algebra" begin 
     @testset "BandedMatrix" begin
         A = BandedMatrix(-3 => Fill(7/10,∞), -2 => 1:∞, 1 => Fill(2im,∞))
         @test A isa BandedMatrix{ComplexF64}
@@ -170,7 +170,7 @@ end
     # BlockTridiagonal(Zeros.(1:∞,2:∞),
     #         (n -> Diagonal(((n+2).+(0:n)))/ (2n + 2)).(0:∞),
     #         Zeros.(2:∞,1:∞))
-end
+# end
 
 include("test_hessenbergq.jl")
 include("test_infql.jl")

test/test_infqr.jl

@@ -1,8 +1,10 @@
-using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, InfiniteArrays, MatrixFactorizations, LazyArrays, FillArrays, SpecialFunctions, Test
-import LazyArrays: colsupport, rowsupport, MemoryLayout, DenseColumnMajor, TriangularLayout, resizedata!
+using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, InfiniteArrays, MatrixFactorizations, LazyArrays,
+        FillArrays, SpecialFunctions, Test, SemiseparableMatrices
+import LazyArrays: colsupport, rowsupport, MemoryLayout, DenseColumnMajor, TriangularLayout, resizedata!, arguments
 import LazyBandedMatrices: BroadcastBandedLayout
 import BandedMatrices: _BandedMatrix, _banded_qr!, BandedColumns
-import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
+import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout, adaptiveqr
+import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout
 
 
 @testset "Adaptive QR" begin
@@ -26,7 +28,7 @@ import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
     end
 
     @testset "AdaptiveQRFactors" begin
-        A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1) 
+        A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1)
         F = qr(A);
         @test F.factors[1,1] ≈ -sqrt(2)
         @test F.factors[100,100] ≈ qrunblocked(A[1:101,1:100]).factors[100,100]
@@ -39,7 +41,7 @@ import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
     end
 
     @testset "col/rowsupport" begin
-        A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1) 
+        A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1)
         F = qr(A);
         @test MemoryLayout(typeof(F.factors)) isa AdaptiveLayout{BandedColumns{DenseColumnMajor}}
         @test bandwidths(F.factors) == (1,2)
@@ -60,7 +62,7 @@ import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
     end
 
     @testset "Qmul" begin
-        A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1) 
+        A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1)
         Q,R = qr(A);
         b = Vcat([1.,2,3],Zeros(∞))
         @test lmul!(Q, Base.copymutable(b)).datasize[1] == 4
@@ -91,7 +93,7 @@ import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
         @test qr(A[1:3000,1:3000]).Q'b[1:3000] ≈ (F.Q'b)[1:3000]
         @time J = A \ Vcat([besselj(1,z)], Zeros(∞))
         @test J[1:2000] ≈ [besselj(k,z) for k=0:1999]
-    
+
         z = 10_000; # the bigger z the longer before we see convergence
         A = BandedMatrix(0 => -2*(0:∞)/z, 1 => Ones(∞), -1 => Ones(∞))
         @time J = A \ Vcat([besselj(1,z)], Zeros(∞))
@@ -132,4 +134,77 @@ import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
         @test F.Q[1:10,1:10] == Eye(10)
         @test F.R[1:10,1:10] == A[1:10,1:10]
     end
+
+    @testset "almost-banded" begin
+        @testset "one-band" begin
+            A = Vcat(Ones(1,∞), BandedMatrix(0 => -Ones(∞), 1 => 1:∞))
+            @test MemoryLayout(typeof(A)) isa VcatAlmostBandedLayout
+            V = view(A,1:10,1:10)
+            @test MemoryLayout(typeof(V)) isa VcatAlmostBandedLayout
+            @test A[1:10,1:10] isa AlmostBandedMatrix
+            @test AlmostBandedMatrix(V) == Matrix(V) == A[1:10,1:10]
+
+            C = cache(A);
+            @test C[1000,1000] ≡ 999.0
+            F = adaptiveqr(A);
+            partialqr!(F.factors.data,2);
+            @test F.factors.data.data[1:3,1:5] ≈ qr(A[1:3,1:5]).factors
+            partialqr!(F.factors.data,3);
+            @test F.factors.data.data[1:4,1:6] ≈ qr(A[1:4,1:6]).factors
+
+            F = adaptiveqr(A);
+            partialqr!(F.factors.data,10);
+            @test F.factors[1:11,1:10] ≈ qr(A[1:11,1:10]).factors
+            @test F.τ[1:10] ≈ qr(A[1:11,1:10]).τ
+            partialqr!(F.factors.data,20);
+            @test F.factors[1:21,1:20] ≈ qr(A[1:21,1:20]).factors
+
+            @test adaptiveqr(A).R[1:10,1:10] ≈ qr(A[1:11,1:10]).R
+
+            @test qr(A) isa MatrixFactorizations.QR{Float64,<:InfiniteLinearAlgebra.AdaptiveQRFactors}
+            @test factorize(A) isa MatrixFactorizations.QR{Float64,<:InfiniteLinearAlgebra.AdaptiveQRFactors}
+
+            @test (adaptiveqr(A) \ [ℯ; zeros(∞)])[1:1000] ≈ (qr(A) \ [ℯ; zeros(∞)])[1:1000] ≈ (A \ [ℯ; zeros(∞)])[1:1000] ≈ [1/factorial(1.0k) for k=0:999]
+        end
+
+        @testset "two-bands" begin
+            B = BandedMatrix(0 => -Ones(∞), 2 => (1:∞).* (2:∞))
+            A = Vcat(Vcat(Ones(1,∞), ((-1).^(0:∞))'), B)
+            @test MemoryLayout(typeof(A)) isa VcatAlmostBandedLayout
+
+            @test qr(A) isa MatrixFactorizations.QR{Float64,<:InfiniteLinearAlgebra.AdaptiveQRFactors}
+            @test factorize(A) isa MatrixFactorizations.QR{Float64,<:InfiniteLinearAlgebra.AdaptiveQRFactors}
+            u = qr(A) \ [1; zeros(∞)]
+            x = 0.1
+            @test (exp(1 - x)*(-1 + exp(2 + 2x)))/(-1 + exp(4)) ≈ dot(u[1:1000], x.^(0:999))
+            u = qr(A) \ Vcat([ℯ,1/ℯ], zeros(∞))
+            @test u[1:1000] ≈ [1/factorial(1.0k) for k=0:999]
+            u = qr(A) \ Vcat(ℯ,1/ℯ, zeros(∞))
+            @test u[1:1000] ≈ [1/factorial(1.0k) for k=0:999]
+
+            A = Vcat(Ones(1,∞), ((-1.0).^(0:∞))', B)
+            @test MemoryLayout(typeof(A)) isa VcatAlmostBandedLayout
+            u = A \ Vcat(ℯ,1/ℯ, zeros(∞))
+            @test u[1:1000] ≈ [1/factorial(1.0k) for k=0:999]
+
+            A = Vcat(Ones(1,∞), ((-1).^(0:∞))', B)
+            u = A \ Vcat(ℯ,1/ℯ, zeros(∞))
+            @test u[1:1000] ≈ [1/factorial(1.0k) for k=0:999]            
+        end
+
+        @testset "more bands" begin
+            L = Vcat(Ones(1,∞), ((-1).^(0:∞))', 
+                     BandedMatrix(-1 => Ones(∞), 1 => Ones(∞), 2 => 4:2:∞, 3 => Ones(∞), 5 => Ones(∞)))
+            F = qr(L).factors.data;
+            resizedata!(F.data,13,19)
+            @test F.data.data[2,8] == -1
+            F = qr(L);
+            partialqr!(F,10);
+            @test F.factors[1:10,1:10] ≈ qr(L[1:13,1:10]).factors[1:10,1:10]
+            @test qr(L).factors[1:10,1:10] ≈ qr(L[1:13,1:10]).factors[1:10,1:10]
+
+            u = L \ [1; 2; zeros(∞)]
+            @test L[1:1000,1:1000]*u[1:1000] ≈ [1; 2; zeros(998)]
+        end
+    end
 end