Cholesky decomposition of low-rank positive-semidefinite matrix

[Shuvomoy_Das_Gupta]

Dear All,

I have a rank r positive-semidefinite matrix M, which I am trying to decompose via lower-triangularlar Cholesky factorization M=LL^\top, but I am particularly interested in computing an L that has r positive diagonal entries and n-r columns containing all zero. Such a Cholesky matrix exists due to Theorem 10.9 of [Higham02] (I have stated the result below for completeness).

Is there any way I can do this in Julia using any of the available Cholesky decomposition functions? [Higham 02, page 202] presents an outer product Cholesky algorithm, comprising of r stages, where at each state a rank-1 matrix is subtracted from the matrix iterate. Is there any function in Julia that has a similar implementation? Any tips or suggestions will be much appreciated.

[Higham 02, Theorem 10.9] states that If a positive semidefinite matrix M\in\mathbf{S}^{n} has rank r, then there exists a lower-triangular Cholesky matrix L satisfying M=LL^{\top} with r positive diagonal entries and n-r columns containing all zero.

[Higham02] Higham, Nicholas J. Accuracy and stability of numerical algorithms . Society for industrial and applied mathematics, 2002.

[stevengj]

You maybe want the “pivoted Cholesky” factorization, via cholesky(A, Val(true)), probably passing check=false. (Under the hood, this calls LAPACK’s dpstrf.)

The central complication is that a positive-semidefinite matrix can easily become indefinite (slightly negative eigenvalues or pivots) due to roundoff errors. You can combat this to some extent with pivoted Cholesky by passing a tol argument or a check=false argument to control what the algorithm does when a negative pivot is encountered. See the dpstrf LAPACK documentation for more detail — the cholesky(A, Val(true)) documentation currently isn’t so clear.

[stevengj]

For example, here is a rank-3 5x5 positive semidefinite matrix:

julia> B = rand(3,5); A = Hermitian(B'B);

cholesky(A) throws PosDefException, and cholesky(A, Val(true)) throws RankDeficientException. However, passing check=false forces the factorization to proceed even if it is rank-deficient:

julia> F = cholesky(A, Val(true); check=false)
CholeskyPivoted{Float64, Matrix{Float64}}
U factor with rank 4:
5×5 UpperTriangular{Float64, Matrix{Float64}}:
 1.26174  0.830057  1.04195    0.862304     0.981338
  ⋅       0.708177  0.242843  -0.34719     -0.0424352
  ⋅        ⋅        0.573781   0.122706     0.375544
  ⋅        ⋅         ⋅         1.05367e-8  -1.05367e-8
  ⋅        ⋅         ⋅          ⋅          -2.22045e-16
permutation:
5-element Vector{Int64}:
 3
 1
 4
 5
 2

Notice the slightly negative pivot due to roundoff errors, as well as the small 4th pivot (on the order of the square root of the precision — recall that Cholesky takes square roots of pivots). This is still a factorization of A (up to roundoff), taking into account the permutation F.p:

julia> F.L * F.L' - A[F.p, F.p]
5×5 Matrix{Float64}:
  0.0          0.0  -2.22045e-16  0.0  0.0
  0.0          0.0   0.0          0.0  0.0
 -2.22045e-16  0.0   0.0          0.0  0.0
  0.0          0.0   0.0          0.0  0.0
  0.0          0.0   0.0          0.0  2.22045e-16

However, in general the pivoted-Cholesky algorithm only factorizes a submatrix of A[F.p, F.p], up to the point where the algorithm terminates.

[Shuvomoy_Das_Gupta]

Thanks so much, @stevengj for the very clear answer! Also, thanks for the dpstrf pointer, I will go through the LAPACK documentation.