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.
1 Like
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.
3 Likes
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.
4 Likes
Thanks so much, @stevengj for the very clear answer! Also, thanks for the dpstrf
pointer, I will go through the LAPACK documentation.