Implement reverse Cholesky for symmetric Tridiagonal matrices (#179)

* Implement reverse Cholesky for symmetric Tridiagonal matrices

* Avoid infinite loop if bN is somehow zero (if it can be?)

* Address review

* Fix tests

Author: DanielVandH

src/InfiniteLinearAlgebra.jl

@@ -138,6 +138,7 @@ include("banded/infbanded.jl")
 include("blockbanded/blockbanded.jl")
 include("banded/infqltoeplitz.jl")
 include("banded/infreversecholeskytoeplitz.jl")
+include("banded/infreversecholeskytridiagonal.jl")
 include("infql.jl")
 include("infqr.jl")
 include("inful.jl")

src/banded/infreversecholeskytridiagonal.jl

@@ -0,0 +1,141 @@
+const MAX_TRIDIAG_CHOL_N = 2^21 - 1 # Maximum allowable size of Cholesky factor before terminating to prevent OutOfMemory errors without convergence
+mutable struct LazySymTridiagonalReverseCholeskyFactors{T,M1,M2} <: LazyMatrix{T}
+    const A::M1 # original matrix
+    const L::M2 # LN
+    const ε::T  # adaptive tolerance
+    N::Int      # size of L used for approximating Ln     
+    n::Int      # size of approximated finite section
+end # this should behave like a lower Bidiagonal matrix
+function LazySymTridiagonalReverseCholeskyFactors(A, N, n, L, ε)
+    require_one_based_indexing(A)
+    M1, M2 = typeof(A), typeof(L)
+    T = eltype(L)
+    return LazySymTridiagonalReverseCholeskyFactors{T,M1,M2}(A, L, convert(T, ε), N, n)
+end
+
+function getproperty(C::ReverseCholesky{<:Any,<:LazySymTridiagonalReverseCholeskyFactors}, d::Symbol) # mimic getproperty(::ReverseCholesky{<:Any, <:Bidiagonal}, ::Symbol)
+    Cfactors = getfield(C, :factors)
+    #Cuplo    = 'L' 
+    if d == :U
+        return Cfactors'
+    elseif d == :L || d == :UL
+        return Cfactors
+    else
+        return getfield(C, d)
+    end
+end
+MemoryLayout(::Type{<:LazySymTridiagonalReverseCholeskyFactors}) = BidiagonalLayout{LazyLayout,LazyLayout}()
+
+size(L::LazySymTridiagonalReverseCholeskyFactors) = size(L.A)
+axes(L::LazySymTridiagonalReverseCholeskyFactors) = axes(L.A)
+Base.eltype(L::Type{LazySymTridiagonalReverseCholeskyFactors}) = eltype(L.L)
+
+copy(L::LazySymTridiagonalReverseCholeskyFactors) = LazySymTridiagonalReverseCholeskyFactors(copy(L.A), copy(L.L), L.ε, L.N, L.n)
+copy(U::Adjoint{T,<:LazySymTridiagonalReverseCholeskyFactors}) where {T} = copy(parent(U))'
+
+LazyBandedMatrices.bidiagonaluplo(L::LazySymTridiagonalReverseCholeskyFactors) = 'L'
+
+"""
+    InfiniteBoundsAccessError <: Exception 
+
+Struct for defining an error when accessing a `LazySymTridiagonalReverseCholeskyFactors` object outside of the 
+maximum allowable finite section of size `$MAX_TRIDIAG_CHOL_N × $MAX_TRIDIAG_CHOL_N`.
+"""
+struct InfiniteBoundsAccessError <: Exception
+    i::Int
+    j::Int
+end
+function Base.showerror(io::IO, err::InfiniteBoundsAccessError)
+    print(io, "InfiniteBoundsAccessError: Tried to index reverse Cholesky factory at index (", err.i, ", ", err.j)
+    print(io, "), outside of the maximum allowable finite section of size (", MAX_TRIDIAG_CHOL_N, " × ", MAX_TRIDIAG_CHOL_N, ")")
+end
+
+function getindex(L::LazySymTridiagonalReverseCholeskyFactors, i::Int, j::Int)
+    max(i, j) > MAX_TRIDIAG_CHOL_N && throw(InfiniteBoundsAccessError(i, j))
+    T = eltype(L)
+    if j > i
+        return zero(T)
+    elseif max(i, j) > L.n
+        _expand_factor!(L, max(i, j))
+        return L.L[i, j]
+    else
+        return L.L[i, j]
+    end
+end
+
+function reversecholesky_layout(::SymTridiagonalLayout, ::NTuple{2,OneToInf{Int}}, A, ::NoPivot; kwds...)
+    a, b = A.dv, A.ev
+    T = promote_type(eltype(a), eltype(b), eltype(b[1] / a[1])) # could also use promote_op(/, eltype(a), eltype(b)), but promote_op is fragile apparently 
+    tol = eps(real(T)) # no good way to pass this as a keyword currently, so just hardcode it
+    L = Bidiagonal([zero(T)], T[], :L)
+    chol = LazySymTridiagonalReverseCholeskyFactors(A, 1, 1, L, tol)
+    _expand_factor!(chol, 2^4) # initialise with 2^4 
+    return ReverseCholesky(chol, 'L', 0)
+end
+
+function _expand_factor!(L::LazySymTridiagonalReverseCholeskyFactors, n)
+    L.n ≥ n && return L
+    return __expand_factor!(L::LazySymTridiagonalReverseCholeskyFactors, n)
+end
+
+function compute_ξ(LL::LazySymTridiagonalReverseCholeskyFactors)
+    #=
+    We can show that ||LN' Pn inv(LN') Vb|| = |bN| ||ξ||, where 
+        ξₙ = LN[n, n]νₙ,
+        ξᵢ = LN[i, i]νᵢ + LN[i+1, i]νᵢ₊₁, i = 1, 2, …, n-1, where 
+        νN = 1/LN[N, N] 
+        νᵢ = -(LN[i+1, i] / LN[i, i]) νᵢ₊₁, i = 1, 2, …, N-1.
+    =#
+    L, N, n = LL.L, LL.N, LL.n
+    ν = inv(L[N, N])
+    for i in (N-1):-1:n
+        ν *= -(L[i+1, i] / L[i, i])
+    end
+    ξ = (L[n, n] * ν)^2
+    for i in (n-1):-1:1
+        ξ′ = L[i+1, i] * ν
+        ν *= -(L[i+1, i] / L[i, i])
+        ξ′ += L[i, i] * ν
+        ξ += ξ′^2
+    end
+    bN = LL.A[N, N+1]
+    scale = iszero(bN) ? one(bN) : abs(bN)
+    return scale * sqrt(ξ) # could maybe just return sqrt(ξ), but maybe bN helps for scaling?
+end
+
+function has_converged(LL::LazySymTridiagonalReverseCholeskyFactors)
+    ξ = compute_ξ(LL)
+    return ξ ≤ LL.ε
+end
+
+function _resize_factor!(L::LazySymTridiagonalReverseCholeskyFactors, N=2L.N)
+    L.N = N
+    resize!(L.L.dv, L.N)
+    resize!(L.L.ev, L.N - 1)
+    return L
+end
+
+function _finite_revchol!(L::LazySymTridiagonalReverseCholeskyFactors)
+    # Computes the reverse Cholesky factorisation of L.A[1:L.N, 1:L.N]
+    N = L.N
+    a, b = L.A.dv, L.A.ev
+    ℓa, ℓb = L.L.dv, L.L.ev
+    ℓa[N] = sqrt(a[N])
+    for i in (N-1):-1:1
+        ℓb[i] = b[i] / ℓa[i+1]
+        ℓa[i] = sqrt(a[i] - ℓb[i]^2)
+    end
+    return L
+end
+
+function __expand_factor!(L::LazySymTridiagonalReverseCholeskyFactors, n)
+    L.N > MAX_TRIDIAG_CHOL_N && return L
+    L.n = n
+    L.N < L.n && _resize_factor!(L, 2n)
+    while !has_converged(L) && L.N ≤ MAX_TRIDIAG_CHOL_N
+        _resize_factor!(L)
+        _finite_revchol!(L)
+    end
+    !has_converged(L) && L.N > MAX_TRIDIAG_CHOL_N && @warn "Reverse Cholesky algorithm failed to converge. Returned results may not be accurate." maxlog = 1
+    return L
+end

test/test_infbanded.jl

@@ -63,13 +63,13 @@ using Base: oneto
         T = LazyBandedMatrices.Tridiagonal(Fill(1,∞), Zeros(∞), Fill(3,∞))
         @test T[2:∞,3:∞] isa SubArray
         @test exp.(T) isa BroadcastMatrix
-        @test exp.(T)[2:∞,3:∞] isa SubArray
+        @test exp.(T)[2:∞,3:∞] isa BroadcastMatrix
         @test exp.(T[2:∞,3:∞]) isa BroadcastMatrix
 
         B = LazyBandedMatrices.Bidiagonal(Fill(1,∞), Zeros(∞), :U)
         @test B[2:∞,3:∞] isa SubArray
         @test exp.(B) isa BroadcastMatrix
-        @test exp.(B)[2:∞,3:∞] isa SubArray
+        @test exp.(B)[2:∞,3:∞] isa BroadcastMatrix
 
         @testset "algebra" begin
             T = Tridiagonal(Fill(1,∞), Fill(2,∞), Fill(3,∞))

test/test_infreversecholesky.jl

@@ -1,16 +1,60 @@
-using InfiniteLinearAlgebra, MatrixFactorizations, ArrayLayouts, LinearAlgebra, Test
+using InfiniteLinearAlgebra, LazyBandedMatrices, FillArrays, MatrixFactorizations, ArrayLayouts, LinearAlgebra, Test, LazyArrays
 
-
-
-@testset "infreversecholesky" begin
+@testset "infreversecholeskytoeplitz" begin
     @testset "Tri Toeplitz" begin
         A = SymTridiagonal(Fill(3, ∞), Fill(1, ∞))
         U, = reversecholesky(A)
-        @test (U*U')[1:10,1:10] ≈ A[1:10,1:10]
+        @test (U*U')[1:10, 1:10] ≈ A[1:10, 1:10]
     end
 
     @testset "Pert Tri Toeplitz" begin
-        A = SymTridiagonal([[4,5, 6]; Fill(3, ∞)], [[2,3]; Fill(1, ∞)])
-        @test reversecholesky(A).U[1:100,1:100] ≈ reversecholesky(A[1:1000,1:1000]).U[1:100,1:100]
+        A = SymTridiagonal([[4, 5, 6]; Fill(3, ∞)], [[2, 3]; Fill(1, ∞)])
+        @test reversecholesky(A).U[1:100, 1:100] ≈ reversecholesky(A[1:1000, 1:1000]).U[1:100, 1:100]
+    end
+end
+
+@testset "infreversecholeskytridiagonal" begin
+    local LL, L
+    @testset "Test on Toeplitz example first" begin
+        A = SymTridiagonal(Fill(3, ∞), Fill(1, ∞))
+        L = reversecholesky(A)
+        LL = InfiniteLinearAlgebra.reversecholesky_layout(SymTridiagonalLayout{LazyArrays.LazyLayout,LazyArrays.LazyLayout}(), axes(A), A, NoPivot())
+        @test L.L[1:1000, 1:1000] ≈ LL.L[1:1000, 1:1000]
+        @test (LL.L'*LL.L)[1:1000, 1:1000] == (LL.U*LL.L)[1:1000, 1:1000] ≈ A[1:1000, 1:1000]
+    end
+
+    @testset "Basic tests" begin
+        L = LL
+        @test MemoryLayout(L.L) isa BidiagonalLayout
+        @test L.U === L.L'
+        @test L.uplo == 'L'
+        @test L.info == 0
+        @test size(L) == (∞, ∞)
+        @test axes(L) == (1:∞, 1:∞)
+        @test eltype(L) == Float64
+        Lc = copy(L)
+        @test !(Lc === L)
+        @test !(Lc.U === L.U)
+        @test Lc.L[1:1000, 1:1000] == L.L[1:1000, 1:1000]
+        UUc = copy(L.L')
+        @test !(UUc === L.U)
+        @test UUc[1:1000, 1:1000] == L.U[1:1000, 1:1000]
+    end
+
+    @testset "Errors" begin
+        err = InfiniteLinearAlgebra.InfiniteBoundsAccessError(4, 6)
+        @test_throws InfiniteLinearAlgebra.InfiniteBoundsAccessError throw(err)
+        @test_throws InfiniteLinearAlgebra.InfiniteBoundsAccessError L.L[1, InfiniteLinearAlgebra.MAX_TRIDIAG_CHOL_N+1]
+        @test_throws InfiniteLinearAlgebra.InfiniteBoundsAccessError L.L[InfiniteLinearAlgebra.MAX_TRIDIAG_CHOL_N+1, 1]
+        @test_throws InfiniteLinearAlgebra.InfiniteBoundsAccessError L.L[InfiniteLinearAlgebra.MAX_TRIDIAG_CHOL_N+1, InfiniteLinearAlgebra.MAX_TRIDIAG_CHOL_N+1]
+    end
+
+    @testset "Another example" begin
+        A = LazyBandedMatrices.SymTridiagonal(Ones(∞), 1 ./ (2:∞))
+        L = reversecholesky(A)
+        @test (L.U*L.L)[1:1000, 1:1000] ≈ A[1:1000, 1:1000]
+        A = 5I + LazyBandedMatrices.SymTridiagonal(1 ./ (2:∞), Ones(∞))
+        L = reversecholesky(A)
+        @test (L.U*L.L)[1:1000, 1:1000] ≈ A[1:1000, 1:1000]
     end
 end