Adaptive Cholesky (#76)

* Start adaptive Cholesky

* ∞ cholesky works

* adaptive cholesky works

* Speed up Cholesky solve

* Update test_infcholesky.jl

* test pass

* update Project

* Update Project.toml

Author: dlfivefifty

.github/workflows/ci.yml

@@ -11,7 +11,7 @@ jobs:
       matrix:
         version:
           - '1.5'
-          - '^1.6.0-0'
+          - '1'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.5.5"
+version = "0.5.6"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -16,13 +16,13 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 
 [compat]
-ArrayLayouts = "0.6"
-BandedMatrices = "0.16"
-BlockArrays = "0.14.5, 0.15"
+ArrayLayouts = "0.7"
+BandedMatrices = "0.16.9"
+BlockArrays = "0.15"
 BlockBandedMatrices = "0.10"
 FillArrays = "0.11"
 InfiniteArrays = "0.10.6"
-LazyArrays = "0.20.2, 0.21"
+LazyArrays = "0.21"
 LazyBandedMatrices = "0.5"
 MatrixFactorizations = "0.8"
 SemiseparableMatrices = "0.2.3"

src/InfiniteLinearAlgebra.jl

@@ -9,12 +9,12 @@ 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
+                        _bidiag_forwardsub!, mulreduce, RangeCumsum, _factorize, transposelayout, ldiv!, lmul!, mul
 import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns,
                         _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: lmul!,  ldiv!, matprod, qr, AbstractTriangular, AbstractQ, adjoint, transpose
+import LinearAlgebra: matprod, qr, AbstractTriangular, AbstractQ, adjoint, transpose, AdjOrTrans
 import LazyArrays: applybroadcaststyle, CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, ApplyStyle, LazyArrayApplyStyle, LazyArrayStyle,
                     resizedata!, MemoryLayout,
                     factorize, sub_materialize, LazyLayout, LazyArrayStyle, layout_getindex,
@@ -25,7 +25,7 @@ import MatrixFactorizations: ul, ul!, _ul, ql, ql!, _ql, QLPackedQ, getL, getR,
 
 import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange, blockcolsupport, BlockLayout, AbstractBlockLayout, BlockSlice
 
-import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!, _banded_qr, _BandedMatrix
+import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!, _banded_qr, _BandedMatrix, banded_chol!
 
 import LazyBandedMatrices: ApplyBandedLayout, BroadcastBandedLayout, _krontrav_axes, _block_interlace_axes, LazyBandedLayout,AbstractLazyBandedLayout
 
@@ -66,6 +66,7 @@ include("banded/infqltoeplitz.jl")
 include("infql.jl")
 include("infqr.jl")
 include("inful.jl")
+include("infcholesky.jl")
 
 
 end # module

src/banded/infqltoeplitz.jl

@@ -88,6 +88,9 @@ struct ProductQ{T,QQ<:Tuple} <: AbstractQ{T}
     Qs::QQ
 end
 
+ArrayLayouts.@layoutmatrix ProductQ
+ArrayLayouts.@_layoutlmul ProductQ
+
 ProductQ(Qs::AbstractMatrix...) = ProductQ{mapreduce(eltype,promote_type,Qs),typeof(Qs)}(Qs)
 
 adjoint(Q::ProductQ) = ProductQ(reverse(map(adjoint,Q.Qs))...)
@@ -109,8 +112,7 @@ function _productq_mul(A::ProductQ{T}, x::AbstractVector{S}) where {T,S}
     lmul!(A, Base.copymutable(convert(AbstractVector{TS},x)))
 end
 
-(*)(A::ProductQ, x::AbstractVector) = _productq_mul(A, x)
-(*)(A::ProductQ, x::LazyVector) = _productq_mul(A, x)
+mul(A::ProductQ, x::AbstractVector) = _productq_mul(A, x)
 
 
 # LQ where Q is a product of orthogonal operations

src/infcholesky.jl

@@ -0,0 +1,101 @@
+
+mutable struct AdaptiveCholeskyFactors{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: LayoutMatrix{T}
+    data::CachedMatrix{T,DM,M}
+    ncols::Int
+end
+
+size(U::AdaptiveCholeskyFactors) = size(U.data.array)
+bandwidths(A::AdaptiveCholeskyFactors) = (0,bandwidth(A.data,2))
+AdaptiveCholeskyFactors(A::Symmetric) = AdaptiveCholeskyFactors(cache(parent(A)),0)
+MemoryLayout(::Type{AdaptiveCholeskyFactors{T,DM,M}}) where {T,DM,M} = AdaptiveLayout{typeof(MemoryLayout(DM))}()
+
+
+function partialcholesky!(F::AdaptiveCholeskyFactors{T,<:BandedMatrix}, n::Int) where T
+    if n > F.ncols 
+        _,u = bandwidths(F.data.array)
+        resizedata!(F.data,n+u,n+u);
+        ncols = F.ncols
+        kr = ncols+1:n
+        factors = view(F.data.data,kr,kr)
+        banded_chol!(factors, UpperTriangular)
+        # multiply remaining columns
+        U1 = UpperTriangular(view(F.data.data,n-u+1:n,n-u+1:n))
+        B = view(F.data.data,n-u+1:n,n+1:n+u)
+        ldiv!(U1',B)
+        muladd!(-one(T),B',B,one(T),view(F.data.data,n+1:n+u,n+1:n+u))
+        F.ncols = n
+    end
+    F
+end
+
+function getindex(F::AdaptiveCholeskyFactors, k::Int, j::Int)
+    partialcholesky!(F, max(k,j))
+    F.data.data[k,j]
+end
+
+
+adaptivecholesky(A) = Cholesky(AdaptiveCholeskyFactors(A), :U, 0)
+
+
+ArrayLayouts._cholesky(::SymmetricLayout{<:AbstractBandedLayout}, ::NTuple{2,OneToInf{Int}}, A) = adaptivecholesky(A)
+
+function colsupport(F::AdaptiveCholeskyFactors, j)
+    partialcholesky!(F, maximum(j)+bandwidth(F,2))
+    colsupport(F.data.data, j)
+end
+
+function rowsupport(F::AdaptiveCholeskyFactors, j)
+    partialcholesky!(F, maximum(j)+bandwidth(F,2))
+    rowsupport(F.data.data, j)
+end
+
+colsupport(F::AdjOrTrans{<:Any,<:AdaptiveCholeskyFactors}, j) = rowsupport(parent(F), j)
+rowsupport(F::AdjOrTrans{<:Any,<:AdaptiveCholeskyFactors}, j) = colsupport(parent(F), j)
+
+function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N',<:AdaptiveLayout},<:PaddedLayout})
+    A,B = M.A,M.B
+    T = eltype(M)
+    COLGROWTH = 1000 # rate to grow columns
+    tol = floatmin(real(T))
+
+    require_one_based_indexing(B)
+    mA, nA = size(A)
+    mB = length(B)
+    if mA != mB
+        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
+    end
+    sB = B.datasize[1]
+    l,_ = bandwidths(A)
+
+    jr = 1:min(COLGROWTH,nA)
+
+    P = parent(parent(A))
+
+    @inbounds begin
+        while first(jr) < nA
+            j = first(jr)
+            cs = colsupport(A, last(jr))
+            cs_max = maximum(cs)
+            kr = j:cs_max
+            resizedata!(B, min(cs_max,mB))
+            if (j > sB && maximum(abs,view(B.data,j:last(colsupport(P,j)))) ≤ tol)
+                break
+            end
+            partialcholesky!(P, min(cs_max,nA))
+            U_N = UpperTriangular(view(P.data.data, jr, jr))
+            ldiv!(U_N', view(B.data, jr))
+            jr1 = last(jr)-l+1:last(jr)
+            jr2 = last(jr)+1:last(jr)+l
+            muladd!(-one(T), view(P.data.data, jr1,jr2)', view(B.data,jr1), one(T), view(B.data,jr2))
+            jr = last(jr)+1:min(last(jr)+COLGROWTH,nA)
+        end
+    end
+    B
+end
+
+function ldiv!(R::UpperTriangular{<:Any,<:AdaptiveCholeskyFactors}, B::CachedVector{<:Any,<:Any,<:Zeros{<:Any,1}})
+    n = B.datasize[1]
+    partialcholesky!(parent(R), n)
+    ldiv!(UpperTriangular(view(parent(R).data.data,oneto(n),oneto(n))), view(B.data,oneto(n)))
+    B
+end

src/infqr.jl

@@ -5,7 +5,7 @@ mutable struct AdaptiveQRData{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}}
     ncols::Int
 end
 
-function AdaptiveQRData(::AbstractBandedLayout, A::AbstractMatrix{T}) where T
+function AdaptiveQRData(::Union{SymmetricLayout{<:AbstractBandedLayout},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)
@@ -94,6 +94,7 @@ 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}()
+transposelayout(A::AdaptiveLayout{ML}) where ML = AdaptiveLayout{typeof(transposelayout(ML()))}()
 
 size(F::AdaptiveQRFactors) = size(F.data.data)
 axes(F::AdaptiveQRFactors) = axes(F.data.data)

test/runtests.jl

@@ -489,3 +489,4 @@ include("test_hessenbergq.jl")
 include("test_infql.jl")
 include("test_infqr.jl")
 include("test_inful.jl")
+include("test_infcholesky.jl")

test/test_infcholesky.jl

@@ -0,0 +1,24 @@
+using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, ArrayLayouts, Test
+
+@testset "infinite-cholesky" begin
+    S = Symmetric(BandedMatrix(0 => 1:∞, 1=> Ones(∞)))
+    b = [1; zeros(∞)]
+    @test cholesky(S) \ b ≈ qr(S) \ b ≈ S \ b
+
+    # go past adaptive 
+    b = [randn(10_000); zeros(∞)]
+    @test cholesky(S) \ b ≈ qr(S) \ b ≈ S \ b
+
+    @testset "singularly perturbed" begin
+        ε = 0.0001
+        S = Symmetric(BandedMatrix(0 => 2 .+ ε*(1:∞), 1=> -Ones(∞)))
+        b = [1; zeros(∞)]
+        @test cholesky(S) \ b ≈  S \ b
+    end
+
+    @testset "long bandwidths" begin
+        S = Symmetric(BandedMatrix(0 => 1:∞, 50=> Ones(∞)))
+        b = [1; zeros(∞)]
+        @test cholesky(S) \ b ≈ qr(S) \ b ≈ S \ b
+    end
+end

test/test_infqr.jl

@@ -148,6 +148,12 @@ import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout
             b = [[1,2,3]; zeros(∞)]
             @test A \ b == [ones(3); zeros(∞)]
         end
+
+        @testset "Symmetric" begin
+            A = Symmetric(BandedMatrix(0 => 1:∞, 1=> Ones(∞)))
+            Ã = BandedMatrix(0 => 1:∞, 1=> Ones(∞), -1=> Ones(∞))
+            @test qr(A).R[1:10,1:10] ≈ qr(Ã).R[1:10,1:10]
+        end
     end
 
     @testset "almost-banded" begin