Implement reverse Cholesky for symmetric Tridiagonal matrices #179

DanielVandH · 2024-06-03T12:49:07Z

This PR implements reverse Cholesky for symmetric Tridiagonal matrices. The algorithm is based on Section 3.2 here. In particular, letting

$$\boldsymbol A = \begin{bmatrix} \boldsymbol A_N & b_N\boldsymbol E_{N1} \\ b_N\boldsymbol E_{1N} & \boldsymbol A_{\infty} \end{bmatrix},$$

where $\boldsymbol A_N \in \mathbb R^{N \times N}$ and $\boldsymbol E_{ij}$ is the matrix whose only non-zero entry is $e_{ij} = 1$, we can show that

$$\|\boldsymbol L_N^\top \boldsymbol P_n\boldsymbol L_N^{-\top} \boldsymbol E_{N1}\| = \|\boldsymbol \xi_n\|,$$

where $\boldsymbol A_N = \boldsymbol L_N^\top \boldsymbol L_N$, $\boldsymbol P_n = \mathrm{diag}(\boldsymbol I_n, \boldsymbol 0)$, $n$ is the size of the finite section we need to access, and

$$\xi_n = [\boldsymbol L_N]_{n, n}\nu_n, \quad \xi_i = [\boldsymbol L_N]_{i,i}\nu_i + [\boldsymbol L_N]_{i+1,i}\boldsymbol \nu_{i+1}, \quad i=1,2,\ldots,n-1,$$

where

$$\nu_N = \frac{1}{[\boldsymbol L_N]_{N,N}}, \quad \nu_i = -\left(\frac{[\boldsymbol L_N]_{i+1,i}}{[\boldsymbol L_N]_{i,i}}\right)\nu_{i+1}, \quad i=1,2,\ldots,N-1.$$

Thus, for a given $\varepsilon > 0$ (here I default to $\varepsilon=|b_N|\sqrt{10^{-5}})$, if $|\boldsymbol \xi_n| < \varepsilon$ then we say that the $n \times n$ finite section has been well approximated and the code for expanding stops. When the factor is indexed at a greater $n$ than previously requested, the factor expands as needed. The reverse Cholesky for finite matrices is computed simply according to

$$\boldsymbol A_n = \begin{bmatrix} a_1 & b_1\boldsymbol e_1^\top \\ b_1\boldsymbol e_1 & \boldsymbol A_{n-1} \end{bmatrix} = \begin{bmatrix} \sqrt{a_1-b_1\ell_{n-1,11}^{-2}} & b_1\ell_{n-1,11}^{-1}\boldsymbol e_1^\top \\ \boldsymbol 0 & \boldsymbol L_{n-1}^\top \end{bmatrix}\begin{bmatrix} \sqrt{a_1 - b_1^2\ell_{n-1,11}^{-2}} & \boldsymbol 0^\top \\ b_1\ell_{n-1,11}^{-1}\boldsymbol e_1 & \boldsymbol L_{n-1} \end{bmatrix},$$

where $\boldsymbol A_{n-1} = \boldsymbol L_{n-1}^\top\boldsymbol L_{n-1}$.

codecov · 2024-06-03T12:57:38Z

Codecov Report

Attention: Patch coverage is 95.45455% with 4 lines in your changes missing coverage. Please review.

Project coverage is 83.03%. Comparing base (f2fe965) to head (b4a38a7).
Report is 1 commits behind head on master.

Files	Patch %	Lines
src/banded/infreversecholeskytridiagonal.jl	95.45%	4 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #179      +/-   ##
==========================================
+ Coverage   82.32%   83.03%   +0.71%     
==========================================
  Files          11       12       +1     
  Lines        1420     1509      +89     
==========================================
+ Hits         1169     1253      +84     
- Misses        251      256       +5

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

dlfivefifty · 2024-06-04T01:30:16Z

@TSGut @MikaelSlevinsky FYI

dlfivefifty

Here are some comments.

I think the unifying the different designs can happen later, separately as a different PR, but it should be on our radar.

I'm leaning towards having the factors be basic types (BandedMatrix, Bidiagonal) that wrap a data type that stores the data but probably it needs more thought.

src/banded/infreversecholeskytridiagonal.jl

test/test_infreversecholesky.jl

DanielVandH · 2024-06-04T09:21:57Z

Not sure what's up with the tests failing. I can't reproduce that error locally...

DanielVandH · 2024-06-04T09:42:25Z

Oh, I can now reproduce it (actually, a different error?) after upgrading to LazyArrays 2.0.4, e.g.

Expression: exp.(T)[2:∞, 3:∞] isa SubArray
  MethodError: BroadcastStyle(::Type{SubArray{Float64, 2, LazyBandedMatrices.Tridiagonal{Float64, Fill{Float64, 1, Tuple{OneToInf{Int64}}}, Zeros{Float64, 1, Tuple{OneToInf{Int64}}}, Fill{Float64, 1, Tuple{OneToInf{Int64}}}}, Tuple{InfiniteArrays.InfUnitRange{Int64}, InfiniteArrays.InfUnitRange{Int64}}, false}}) is ambiguous.

  Candidates:
    BroadcastStyle(::Type{<:SubArray{<:Any, 2, <:BandedMatrices.AbstractBandedMatrix, <:Tuple{AbstractUnitRange{Int64}, AbstractUnitRange{Int64}}}})
      @ BandedMatrices C:\Users\User\.julia\packages\BandedMatrices\as2PJ\src\generic\broadcast.jl:40
    BroadcastStyle(::Type{<:SubArray{<:Any, 2, <:Any, <:Tuple{var"#s2", var"#s1"} where {var"#s2"<:(Union{Base.Slice{OneToInf{T}}, InfiniteArrays.AbstractInfUnitRange{T}, InfiniteArrays.InfStepRange{T}} where T<:Integer), var"#s1"<:(Union{Base.Slice{OneToInf{T}}, InfiniteArrays.AbstractInfUnitRange{T}, InfiniteArrays.InfStepRange{T}} where T<:Integer)}}})
      @ InfiniteArrays C:\Users\User\.julia\packages\InfiniteArrays\7v9AI\src\infrange.jl:489

  Possible fix, define
    BroadcastStyle(::Type{<:SubArray{<:Any, 2, <:BandedMatrices.AbstractBandedMatrix, <:Tuple{var"#s2", var"#s1"} where {var"#s2"<:Union{InfiniteArrays.AbstractInfUnitRange{Int64}, Base.Slice{OneToInf{T}} where T<:Integer}, var"#s1"<:Union{InfiniteArrays.AbstractInfUnitRange{Int64}, Base.Slice{OneToInf{T}} where T<:Integer}}}})

  Stacktrace:
    [1] combine_styles(c::SubArray{Float64, 2, LazyBandedMatrices.Tridiagonal{Float64, Fill{Float64, 1, Tuple{OneToInf{Int64}}}, Zeros{Float64, 1, Tuple{OneToInf{Int64}}}, Fill{Float64, 1, Tuple{OneToInf{Int64}}}}, Tuple{InfiniteArrays.InfUnitRange{Int64}, InfiniteArrays.InfUnitRange{Int64}}, false})
      @ Base.Broadcast .\broadcast.jl:460
    [2] broadcasted
      @ .\broadcast.jl:1341 [inlined]
    [3] _broadcastarray2broadcasted
      @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\lazybroadcasting.jl:73 [inlined]
    [4] _broadcastarray2broadcasted
      @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\lazybroadcasting.jl:79 [inlined]
    [5] _broadcasted
      @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\lazybroadcasting.jl:80 [inlined]
    [6] sub_materialize(::LazyArrays.BroadcastLayout{typeof(exp)}, A::SubArray{Float64, 2, BroadcastMatrix{Float64, typeof(exp), Tuple{LazyBandedMatrices.Tridiagonal{Float64, Fill{Float64, 1, Tuple{OneToInf{Int64}}}, Zeros{Float64, 1, Tuple{OneToInf{Int64}}}, Fill{Float64, 1, Tuple{OneToInf{Int64}}}}}}, Tuple{InfiniteArrays.InfUnitRange{Int64}, InfiniteArrays.InfUnitRange{Int64}}, false})
      @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\lazybroadcasting.jl:97
    [7] sub_materialize
      @ C:\Users\User\.julia\packages\ArrayLayouts\feBnP\src\ArrayLayouts.jl:126 [inlined]
    [8] layout_getindex
      @ C:\Users\User\.julia\packages\ArrayLayouts\feBnP\src\ArrayLayouts.jl:132 [inlined]
    [9] getindex(A::BroadcastMatrix{Float64, typeof(exp), Tuple{LazyBandedMatrices.Tridiagonal{Float64, Fill{Float64, 1, Tuple{OneToInf{Int64}}}, Zeros{Float64, 1, Tuple{OneToInf{Int64}}}, Fill{Float64, 1, Tuple{OneToInf{Int64}}}}}}, kr::InfiniteArrays.InfUnitRange{Int64}, jr::InfiniteArrays.InfUnitRange{Int64})
      @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\feBnP\src\ArrayLayouts.jl:147
   [10] macro expansion
      @ C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Test\src\Test.jl:669 [inlined]
   [11] macro expansion
      @ C:\Users\User\.julia\dev\InfiniteLinearAlgebra.jl\test\test_infbanded.jl:66 [inlined]
   [12] macro expansion
      @ C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Test\src\Test.jl:1577 [inlined]
   [13] macro expansion
      @ C:\Users\User\.julia\dev\InfiniteLinearAlgebra.jl\test\test_infbanded.jl:35 [inlined]
   [14] macro expansion
      @ C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Test\src\Test.jl:1577 [inlined]
   [15] top-level scope
      @ C:\Users\User\.julia\dev\InfiniteLinearAlgebra.jl\test\test_infbanded.jl:7
∞-Toeplitz: Error During Test at C:\Users\User\.julia\dev\InfiniteLinearAlgebra.jl\test\test_infbanded.jl:71
  Test threw exception
  Expression: exp.(B)[2:∞, 3:∞] isa SubArray
  MethodError: BroadcastStyle(::Type{SubArray{Float64, 2, LazyBandedMatrices.Bidiagonal{Float64, Fill{Float64, 1, Tuple{OneToInf{Int64}}}, Zeros{Float64, 1, Tuple{OneToInf{Int64}}}}, Tuple{InfiniteArrays.InfUnitRange{Int64}, InfiniteArrays.InfUnitRange{Int64}}, false}}) is ambiguous.

  Candidates:
    BroadcastStyle(::Type{<:SubArray{<:Any, 2, <:BandedMatrices.AbstractBandedMatrix, <:Tuple{AbstractUnitRange{Int64}, AbstractUnitRange{Int64}}}})
      @ BandedMatrices C:\Users\User\.julia\packages\BandedMatrices\as2PJ\src\generic\broadcast.jl:40
    BroadcastStyle(::Type{<:SubArray{<:Any, 2, <:Any, <:Tuple{var"#s2", var"#s1"} where {var"#s2"<:(Union{Base.Slice{OneToInf{T}}, InfiniteArrays.AbstractInfUnitRange{T}, InfiniteArrays.InfStepRange{T}} where T<:Integer), var"#s1"<:(Union{Base.Slice{OneToInf{T}}, InfiniteArrays.AbstractInfUnitRange{T}, InfiniteArrays.InfStepRange{T}} where T<:Integer)}}})
      @ InfiniteArrays C:\Users\User\.julia\packages\InfiniteArrays\7v9AI\src\infrange.jl:489

  Possible fix, define
    BroadcastStyle(::Type{<:SubArray{<:Any, 2, <:BandedMatrices.AbstractBandedMatrix, <:Tuple{var"#s2", var"#s1"} where {var"#s2"<:Union{InfiniteArrays.AbstractInfUnitRange{Int64}, Base.Slice{OneToInf{T}} where T<:Integer}, var"#s1"<:Union{InfiniteArrays.AbstractInfUnitRange{Int64}, Base.Slice{OneToInf{T}} where T<:Integer}}}})

  Stacktrace:
    [1] combine_styles(c::SubArray{Float64, 2, LazyBandedMatrices.Bidiagonal{Float64, Fill{Float64, 1, Tuple{OneToInf{Int64}}}, Zeros{Float64, 1, Tuple{OneToInf{Int64}}}}, Tuple{InfiniteArrays.InfUnitRange{Int64}, InfiniteArrays.InfUnitRange{Int64}}, false})
      @ Base.Broadcast .\broadcast.jl:460
    [2] broadcasted
      @ .\broadcast.jl:1341 [inlined]
    [3] _broadcastarray2broadcasted
      @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\lazybroadcasting.jl:73 [inlined]
    [4] _broadcastarray2broadcasted
      @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\lazybroadcasting.jl:79 [inlined]
    [5] _broadcasted
      @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\lazybroadcasting.jl:80 [inlined]
    [6] sub_materialize(::LazyArrays.BroadcastLayout{typeof(exp)}, A::SubArray{Float64, 2, BroadcastMatrix{Float64, typeof(exp), Tuple{LazyBandedMatrices.Bidiagonal{Float64, Fill{Float64, 1, Tuple{OneToInf{Int64}}}, Zeros{Float64, 1, Tuple{OneToInf{Int64}}}}}}, Tuple{InfiniteArrays.InfUnitRange{Int64}, InfiniteArrays.InfUnitRange{Int64}}, false})
      @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\lazybroadcasting.jl:97
    [7] sub_materialize
      @ C:\Users\User\.julia\packages\ArrayLayouts\feBnP\src\ArrayLayouts.jl:126 [inlined]
    [8] layout_getindex
      @ C:\Users\User\.julia\packages\ArrayLayouts\feBnP\src\ArrayLayouts.jl:132 [inlined]
    [9] getindex(A::BroadcastMatrix{Float64, typeof(exp), Tuple{LazyBandedMatrices.Bidiagonal{Float64, Fill{Float64, 1, Tuple{OneToInf{Int64}}}, Zeros{Float64, 1, Tuple{OneToInf{Int64}}}}}}, kr::InfiniteArrays.InfUnitRange{Int64}, jr::InfiniteArrays.InfUnitRange{Int64})
      @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\feBnP\src\ArrayLayouts.jl:147
   [10] macro expansion
      @ C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Test\src\Test.jl:669 [inlined]
   [11] macro expansion
      @ C:\Users\User\.julia\dev\InfiniteLinearAlgebra.jl\test\test_infbanded.jl:71 [inlined]
   [12] macro expansion
      @ C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Test\src\Test.jl:1577 [inlined]
   [13] macro expansion
      @ C:\Users\User\.julia\dev\InfiniteLinearAlgebra.jl\test\test_infbanded.jl:35 [inlined]
   [14] macro expansion
      @ C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Test\src\Test.jl:1577 [inlined]
   [15] top-level scope
      @ C:\Users\User\.julia\dev\InfiniteLinearAlgebra.jl\test\test_infbanded.jl:7

dlfivefifty · 2024-06-04T10:29:28Z

I guess you figured out that it was fixed on master?

a change in lazy arrays triggered an ambiguity error unfortunately

DanielVandH · 2024-06-04T10:45:21Z

Yeah, didn't notice that it was updated. I brought the changes in but it didn't merge as I expected so I just applied it again in b4a38a7.

DanielVandH added 2 commits June 3, 2024 13:40

Implement reverse Cholesky for symmetric Tridiagonal matrices

f98b8d2

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

414dd5a

dlfivefifty requested changes Jun 4, 2024

View reviewed changes

Address review

5a54456

DanielVandH and others added 2 commits June 4, 2024 10:44

Merge branch 'JuliaLinearAlgebra:master' into symtri-revchol

5abec47

Fix tests

b4a38a7

DanielVandH mentioned this pull request Jun 4, 2024

Infinite UL factorisation #181

Open

dlfivefifty approved these changes Jun 5, 2024

View reviewed changes

dlfivefifty merged commit 038c020 into JuliaLinearAlgebra:master Jun 5, 2024
6 checks passed

TSGut mentioned this pull request Jun 6, 2024

Prevent runaway QL #183

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Implement reverse Cholesky for symmetric Tridiagonal matrices #179

Implement reverse Cholesky for symmetric Tridiagonal matrices #179

Uh oh!

DanielVandH commented Jun 3, 2024 •

edited

Loading

Uh oh!

codecov bot commented Jun 3, 2024 •

edited

Loading

Uh oh!

dlfivefifty commented Jun 4, 2024

Uh oh!

dlfivefifty left a comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

DanielVandH commented Jun 4, 2024

Uh oh!

DanielVandH commented Jun 4, 2024 •

edited

Loading

Uh oh!

dlfivefifty commented Jun 4, 2024

Uh oh!

DanielVandH commented Jun 4, 2024

Uh oh!

Uh oh!

Uh oh!

Implement reverse Cholesky for symmetric Tridiagonal matrices #179

Implement reverse Cholesky for symmetric Tridiagonal matrices #179

Uh oh!

Conversation

DanielVandH commented Jun 3, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov bot commented Jun 3, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

dlfivefifty commented Jun 4, 2024

Uh oh!

dlfivefifty left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

DanielVandH commented Jun 4, 2024

Uh oh!

DanielVandH commented Jun 4, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dlfivefifty commented Jun 4, 2024

Uh oh!

DanielVandH commented Jun 4, 2024

Uh oh!

Uh oh!

Uh oh!

DanielVandH commented Jun 3, 2024 •

edited

Loading

codecov bot commented Jun 3, 2024 •

edited

Loading

DanielVandH commented Jun 4, 2024 •

edited

Loading