Adaptive QR works!

Author: dlfivefifty

src/InfiniteLinearAlgebra.jl

@@ -1,26 +1,45 @@
 module InfiniteLinearAlgebra
 using BlockArrays, BlockBandedMatrices, BandedMatrices, LazyArrays, FillArrays, InfiniteArrays, MatrixFactorizations, LinearAlgebra
 
-import Base: +, -, *, /, \, OneTo, getindex, promote_op, _unsafe_getindex, print_matrix_row, size,
+import Base: +, -, *, /, \, ^, OneTo, getindex, promote_op, _unsafe_getindex, print_matrix_row, size, axes,
             AbstractMatrix, AbstractArray, Matrix, Array, Vector, AbstractVector,
             show, getproperty
+import Base.Broadcast: BroadcastStyle
+
 import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange
 import FillArrays: AbstractFill
 import BandedMatrices: BandedMatrix, _BandedMatrix, bandeddata
-import LinearAlgebra: lmul!,  ldiv!, matprod, qr, QRPackedQ, AbstractTriangular, AbstractQ, adjoint, transpose,
-                        QR
-import LazyArrays: CachedArray, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes
-import MatrixFactorizations: ql, ql!, QLPackedQ, getL, reflector!, reflectorApply!,
-                            QL
+import LinearAlgebra: lmul!,  ldiv!, matprod, qr, AbstractTriangular, AbstractQ, adjoint, transpose
+import LazyArrays: CachedArray, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, ApplyStyle, LazyArrayApplyStyle, LazyArrayStyle,
+                    CachedMatrix, CachedArray, resizedata!, MemoryLayout
+import MatrixFactorizations: ql, ql!, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ
 
 import BlockArrays: BlockSizes, cumulsizes, _find_block, AbstractBlockVecOrMat, sizes_from_blocks
 
-import BandedMatrices: BandedMatrix, bandwidths
+import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!
 
 import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, AbstractBlockSizes, cumulsizes, _BlockSkylineMatrix, BlockSizes, blockstart, blockstride,
         BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal
 
 
+# Fix ∞ BandedMatrix
+ApplyStyle(::typeof(*), ::Type{<:BandedMatrix{<:Any,<:Any,<:OneToInf}}, _::Type{<:AbstractArray}...) =
+    LazyArrayApplyStyle()
+
+# BroadcastStyle(::Type{<:BandedMatrix{<:Any,<:Any,<:OneToInf}}) = LazyArrayStyle{2}()
+
+^(A::BandedMatrix{T,<:Any,<:OneToInf}, p::Integer) where T =
+    if p < 0 
+        inv(A)^(-p)
+    elseif p == 0
+        Eye{T}(∞)
+    elseif p == 1
+        copy(A)
+    else
+        A*A^(p-1)
+    end
+    
+
 if VERSION < v"1.2-"
     import Base: has_offset_axes
     require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))

src/infqr.jl

@@ -0,0 +1,60 @@
+
+mutable struct AdaptiveQRData{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}}
+    data::CachedMatrix{T,DM,M}
+    τ::Vector{T}
+    ncols::Int
+end
+
+function AdaptiveQRData(::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)
+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)
+        resize!(F.τ,n)
+        ñ = F.ncols
+        factors = view(F.data.data,ñ+1:n+l,ñ+1:n)
+        τ = view(F.τ,ñ+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))
+        F.ncols = n
+    end
+    F
+end
+
+struct AdaptiveQRFactors{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: AbstractMatrix{T}
+    data::AdaptiveQRData{T,DM,M}
+end
+
+size(F::AdaptiveQRFactors) = size(F.data.data)
+function getindex(F::AdaptiveQRFactors, k::Int, j::Int)
+    partialqr!(F.data, j)
+    F.data.data[k,j]
+end
+
+struct AdaptiveQRTau{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: AbstractVector{T}
+    data::AdaptiveQRData{T,DM,M}
+end
+
+size(F::AdaptiveQRTau) = (size(F.data.data,1),)
+
+function getindex(F::AdaptiveQRTau, j::Int)
+    partialqr!(F.data, j)
+    F.data.τ[j]
+end
+
+getR(Q::QR, ::Tuple{OneToInf{Int},OneToInf{Int}}) where T = UpperTriangular(Q.factors)
+
+
+function qr(A::BandedMatrix{<:Any,<:Any,<:OneToInf}) 
+    data = AdaptiveQRData(A)
+    QR(AdaptiveQRFactors(data), AdaptiveQRTau(data))
+end

test/runtests.jl

@@ -40,7 +40,13 @@ import BandedMatrices: bandeddata, _BandedMatrix
     At = BandedMatrix(transpose(A))
     @test Ac[1:10,1:10] ≈ (A')[1:10,1:10] ≈ A[1:10,1:10]'
     @test At[1:10,1:10] ≈ transpose(A)[1:10,1:10] ≈ transpose(A[1:10,1:10])
+
+    A = _BandedMatrix(Fill(1,4,∞),∞,1,2)
+    @test A*A isa ApplyArray
+    @test (A^2)[1:10,1:10] == (A*A)[1:10,1:10] == (A[1:100,1:100]^2)[1:10,1:10]
+    @test (A^3)[1:10,1:10] == (A*A*A)[1:10,1:10] == (A[1:100,1:100]^3)[1:10,1:10]
 end
+
 include("test_hessenbergq.jl")
 include("test_infql.jl")
 

test/test_infqr.jl

@@ -1,44 +1,38 @@
 using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, InfiniteArrays, MatrixFactorizations, LazyArrays, FillArrays
-import BandedMatrices: _BandedMatrix, colsupport
-import MatrixFactorizations: QR
-import LazyArrays: CachedMatrix
-
-struct AdaptiveQRFactors{T,DM<:CachedMatrix{T}} <: AbstractMatrix{T}
-    factors::DM
-    τ::Vector{T}
-end
-
-struct AdaptiveQRTau{T,DM<:CachedMatrix{T}} <: AbstractMatrix{T}
-    factors::DM
-    τ::Vector{T}
+import BandedMatrices: _BandedMatrix, _banded_qr!, colsupport, BandedColumns
+import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData
+
+
+@testset "Adaptive QR" begin
+    @testset "test partialqr!" begin
+        A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1)
+        l,u = bandwidths(A)
+        F = AdaptiveQRData(A);
+        @test bandwidths(F.data) == (1,2)
+        n = 3
+        partialqr!(F,n);
+        F̃ = qrunblocked(Matrix(A[1:n+l,1:n]))
+        @test F̃.factors ≈ F.data[1:n+l,1:n]
+        @test F̃.τ  ≈ F.τ
+        @test triu!(F.data[1:n+1,1:n+2]) ≈ F̃.Q'*A[1:n+1,1:n+2] # test extra columns have been modified
+        n  = 6
+        partialqr!(F,n);
+        F̃ = qrunblocked(Matrix(A[1:n+l,1:n]))
+        @test F̃.factors ≈ F.data[1:n+l,1:n]
+        @test F̃.τ  ≈ F.τ
+        @test triu!(F.data[1:n+1,1:n+2]) ≈ F̃.Q'*A[1:n+1,1:n+2] # test extra columns have been modified
+    end
+
+    @testset "AdaptiveQRFactors" begin
+        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]
+        @test F.τ[1] ≈ 1+sqrt(2)/2
+        @test F.τ[100] ≈ qrunblocked(A[1:101,1:100]).τ[100]
+        Q,R = F;
+        
+    end
 end
 
-
-A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1)
-C = cache(A)
-
-
-V = view(C.data,1:10,1:11)
-V isa BandedMatrices.BandedSubBandedMatrix
-qr!(V)
-qr!(C.data)
-
-F = QR(
-
-
-vierwC.data
-
-@test C[100_000,100_000] === 100_000.0
-
-@time C[100_000,100_000]
-M = 
-@time F = qr!(C.data)
-
-@which C[1:100,1:100]
-
-@which bandwidths(C)
-
-@which cache(A)
-C.data
-
-2
+