Cholesky factorization of a Tridiagonal matrix should not be dense

[ranocha]

The Cholesky factorization $A = G G^T$ of a banded matrix $A$ preserves the lower bandwidth, i.e., $G$ has the same lower bandwidth as $A$. In particular, the Cholesky factorization of a Tridiagonal matrix should be Bidiagonal. However, (at least) Julia 1.10 uses a dense matrix:

julia> using LinearAlgebra

julia> A = Tridiagonal(ones(4), 2 * ones(5), ones(4))
5×5 Tridiagonal{Float64, Vector{Float64}}:
 2.0  1.0   ⋅    ⋅    ⋅ 
 1.0  2.0  1.0   ⋅    ⋅ 
  ⋅   1.0  2.0  1.0   ⋅ 
  ⋅    ⋅   1.0  2.0  1.0
  ⋅    ⋅    ⋅   1.0  2.0

julia> cholesky(A)
Cholesky{Float64, Matrix{Float64}}
U factor:
5×5 UpperTriangular{Float64, Matrix{Float64}}:
 1.41421  0.707107  0.0       0.0       0.0
  ⋅       1.22474   0.816497  0.0       0.0
  ⋅        ⋅        1.1547    0.866025  0.0
  ⋅        ⋅         ⋅        1.11803   0.894427
  ⋅        ⋅         ⋅         ⋅        1.09545

julia> LinearAlgebra.cholcopy(A)
5×5 Matrix{Float64}:
 2.0  1.0  0.0  0.0  0.0
 1.0  2.0  1.0  0.0  0.0
 0.0  1.0  2.0  1.0  0.0
 0.0  0.0  1.0  2.0  1.0
 0.0  0.0  0.0  1.0  2.0

The issue seems to be mostly related to cholcopy, since the in-place factorization looks better (but I did not benchmark it in detail):

julia> cholesky!(A)
Cholesky{Float64, Tridiagonal{Float64, Vector{Float64}}}
U factor:
5×5 UpperTriangular{Float64, Tridiagonal{Float64, Vector{Float64}}}:
 1.41421  0.707107  0.0       0.0       0.0
  ⋅       1.22474   0.816497  0.0       0.0
  ⋅        ⋅        1.1547    0.866025  0.0
  ⋅        ⋅         ⋅        1.11803   0.894427
  ⋅        ⋅         ⋅         ⋅        1.09545

I consider this to be a performance bug since this makes the cholesky decomposition of a Tridiagonal matrix as slow as that of a dense matrix.

[dkarrasch]

I don't think we have a bandwidth-preserving Cholesky-decomposition, it's just that in the in-place method we write zeros into the outer diagonals, which is permitted, but the computations are still performed, so perhaps no run-time difference.

[dkarrasch]

We could easily save memory allocations, though, by making the cholcopy(::Tridiagonal) return a Bidiagonal, as you say.

[araujoms]

The implementation is a bit of a mess: if you call cholesky(::Tridiagonal) a dense Hermitian matrix is created, and LAPACK's potrf is called on it. If instead you call cholesky!(::Tridiagonal) you get instead a generic fallback.

The ideal solution would be to call LAPACK's pttrf instead. Should be easy to do, there's already a wrapper for it.

[dkarrasch]

In the real case, we need the matrix to be symmetric, and for that we have the SymTridiagonal type already, which does create a Hermitian(Bidiagonal). I think we can just add the Tridiagonal case here. I'm working on a PR.

[araujoms]

But cholesky(::SymTridiagonal) goes to the generic fallback, it doesn't use any property to accelerate the factorization.

[araujoms]

To be more concrete, I'm saying we should do this:

function cholesky!(A::Tridiagonal)
    d = real(A.d)
    e = A.du
    LAPACK.pttrf!(d, e)
    map!(sqrt, d)
    @views e .*= d[1:end-1]
    U = Bidiagonal(d, e, :U)
    return Cholesky(U, 'U', convert(BlasInt, 0))
end

(with proper dispatching and error checking, of course)