inf QR Q mul

Author: dlfivefifty

README.md

@@ -2,13 +2,21 @@
 
 A Julia repository for linear algebra with infinite banded and block-banded matrices
 
-This currently implements the infinite-dimensional QL decomposition for perturbations of Toeplitz operators. Here is an example:
+
+## Infinite-dimensional QR factorization
+
+This currently supports the infinite-dimensional QR factorization for banded matrices, also known as the adaptive QR decomposition as the entries of the QR decomposition are determined lazily. 
+
+## Infinite-dimensional QL factorization
+
+
+This currently supports the infinite-dimensional QL factorization for perturbations of Toeplitz operators. Here is an example:
 ```julia
 # Bull head matrix
 A = BandedMatrix(-3 => Fill(7/10,∞), -2 => Fill(1,∞), 1 => Fill(2im,∞))
 ql(A - 5*I)
 ```
-The infinite-dimensional QL decomposition is a subtly thing: its defined when the operator has non-positive Fredholm index, and if the Fredholm index is not zero, it may not be unique. For the Bull head matrix `A`, here are plots of `ql(A-λ*I).L[1,1]` alongside the image of the symbol `A`, which depicts the essential spectrum of `A` and where the Fredholm index changes. Note we have two plots as the regions with negative Fredholm index  have multiple QL decompositions. Where the Fredholm index is positive, the QL decomposition doesn't exist and is depected in black.
+The infinite-dimensional QL factorization is a subtly thing: its defined when the operator has non-positive Fredholm index, and if the Fredholm index is not zero, it may not be unique. For the Bull head matrix `A`, here are plots of `ql(A-λ*I).L[1,1]` alongside the image of the symbol `A`, which depicts the essential spectrum of `A` and where the Fredholm index changes. Note we have two plots as the regions with negative Fredholm index  have multiple QL factorizations. Where the Fredholm index is positive, the QL factorization doesn't exist and is depected in black.
 
 <img src=https://github.com/JuliaMatrices/InfiniteLinearAlgebra.jl/raw/master/images/ql1.png width=500 height=400>
 <img src=https://github.com/JuliaMatrices/InfiniteLinearAlgebra.jl/raw/master/images/ql2.png width=500 height=400>

src/InfiniteLinearAlgebra.jl

@@ -10,8 +10,9 @@ import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange
 import FillArrays: AbstractFill
 import BandedMatrices: BandedMatrix, _BandedMatrix, bandeddata, bandwidths
 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, mulapplystyle, LmulStyle, RmulStyle
+import LazyArrays: CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, ApplyStyle, LazyArrayApplyStyle, LazyArrayStyle,
+                    CachedMatrix, CachedArray, resizedata!, MemoryLayout, mulapplystyle, LmulStyle, RmulStyle,
+                    colsupport, rowsupport, triangularlayout
 import MatrixFactorizations: ql, ql!, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ
 
 import BlockArrays: BlockSizes, cumulsizes, _find_block, AbstractBlockVecOrMat, sizes_from_blocks

src/infqr.jl

@@ -34,13 +34,36 @@ struct AdaptiveQRFactors{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: Abstra
     data::AdaptiveQRData{T,DM,M}
 end
 
+struct AdaptiveLayout{M} <: MemoryLayout end
+MemoryLayout(::Type{AdaptiveQRFactors{T,DM,M}}) where {T,DM,M} = AdaptiveLayout{typeof(MemoryLayout(DM))}()
+triangularlayout(::Type{Tri}, ::ML) where {Tri, ML<:AdaptiveLayout} = Tri{ML}()
+
 size(F::AdaptiveQRFactors) = size(F.data.data)
 bandwidths(F::AdaptiveQRFactors) = bandwidths(F.data.data)
+function colsupport(F::AdaptiveQRFactors, j)
+    partialqr!(F.data, j)
+    colsupport(F.data.data, j)
+end
+
+function rowsupport(F::AdaptiveQRFactors, j)
+    partialqr!(F.data, j+bandwidth(F,2))
+    rowsupport(F.data.data, j)
+end
+
 function getindex(F::AdaptiveQRFactors, k::Int, j::Int)
     partialqr!(F.data, j)
     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)
+
+
+Base.replace_in_print_matrix(A::AdaptiveQRFactors, i::Integer, j::Integer, s::AbstractString) =
+    i in colsupport(A,j) ? s : Base.replace_with_centered_mark(s)
+Base.replace_in_print_matrix(A::UpperTriangular{<:Any,<:AdaptiveQRFactors}, i::Integer, j::Integer, s::AbstractString) =
+    i in colsupport(A,j) ? s : Base.replace_with_centered_mark(s)
+
 struct AdaptiveQRTau{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: AbstractVector{T}
     data::AdaptiveQRData{T,DM,M}
 end
@@ -58,4 +81,37 @@ getR(Q::QR, ::Tuple{OneToInf{Int},OneToInf{Int}}) where T = UpperTriangular(Q.fa
 function qr(A::BandedMatrix{<:Any,<:Any,<:OneToInf}) 
     data = AdaptiveQRData(A)
     QR(AdaptiveQRFactors(data), AdaptiveQRTau(data))
+end
+
+
+function lmul!(A::QRPackedQ{<:Any,<:AdaptiveQRFactors}, B::CachedVector{T,Vector{T},<:Zeros{T,1}}) where T
+    require_one_based_indexing(B)
+    mA, nA = size(A.factors)
+    sB = length(B.data)
+    mB = length(B)
+    if mA != mB
+        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
+    end
+    Afactors = A.factors
+
+    resizedata!(B, maximum(colsupport(A,sB)))
+    b = B.data    
+    
+    sB_2 = length(b)
+
+    @inbounds begin
+        for k = sB:-1:1
+            cs = colsupport(Afactors,k)
+            vBj = b[k]
+            for i = (k+1:sB_2) ∩ cs
+                vBj += conj(Afactors[i,k])*b[i]
+            end
+            vBj = A.τ[k]*vBj
+            b[k] -= vBj
+            for i = (k+1:sB_2) ∩ cs
+                b[i] -= Afactors[i,k]*vBj
+            end
+        end
+    end
+    B
 end

test/test_infqr.jl

@@ -1,6 +1,7 @@
 using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, InfiniteArrays, MatrixFactorizations, LazyArrays, FillArrays
-import BandedMatrices: _BandedMatrix, _banded_qr!, colsupport, BandedColumns
-import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData
+import LazyArrays: colsupport, rowsupport, MemoryLayout, DenseColumnMajor, TriangularLayout
+import BandedMatrices: _BandedMatrix, _banded_qr!, BandedColumns
+import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
 
 
 @testset "Adaptive QR" begin
@@ -35,6 +36,39 @@ import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData
         @test size(R) == (∞,∞)
         @test R[1,1] == -sqrt(2)
     end
+
+    @testset "col/rowsupport" begin
+        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)
+        @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;
+        @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
+    end
+
+    @testset "Qmul" begin
+        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)).data ≈ qr(A[1:4,1:3]).Q*[1,2,3]
+        @test Q[1,1] ≈ -1/sqrt(2)
+        @test Q[200_000,200_000] ≈ -1.0
+        @test Q[1:101,1:100] ≈ qr(A[1:101,1:100]).Q[:,1:100]
+
+        # lmul!(Q', Base.copymutable(b))
+    end
 end