# Cholesky decomposition of low-rank positive-semidefinite matrix
**URL:** https://discourse.julialang.org/t/cholesky-decomposition-of-low-rank-positive-semidefinite-matrix/70397
**Category:** General Usage
**Tags:** linearalgebra, numerics
**Created:** [October 26, 2021, 10:47am UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-low-rank-positive-semidefinite-matrix/70397 "2021-10-26T10:47:20Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@Shuvomoy\_Das\_Gupta](https://discourse.julialang.org/u/Shuvomoy_Das_Gupta)
#### Post date: [October 26, 2021, 10:47am UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-low-rank-positive-semidefinite-matrix/70397/1 "2021-10-26T10:47:20Z")
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.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [October 26, 2021, 12:20pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-low-rank-positive-semidefinite-matrix/70397/2 "2021-10-26T12:20:51Z")
You maybe want the [“pivoted Cholesky” factorization](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.cholesky), via `cholesky(A, Val(true))`, probably passing `check=false`. (Under the hood, this calls LAPACK’s [dpstrf](http://www.netlib.org/lapack/explore-html/da/dba/group__double_o_t_h_e_rcomputational_ga31cdc13a7f4ad687f4aefebff870e1cc.html).)
**The central complication** is that a positive-semidefinite matrix can easily become _in_definite (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.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [October 26, 2021, 12:34pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-low-rank-positive-semidefinite-matrix/70397/3 "2021-10-26T12:34:51Z")
For example, here is a rank-3 5x5 positive semidefinite matrix:
```julia
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
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
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.
---
### Author:  [@Shuvomoy\_Das\_Gupta](https://discourse.julialang.org/u/Shuvomoy_Das_Gupta)
#### Post date: [October 26, 2021, 12:54pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-low-rank-positive-semidefinite-matrix/70397/4 "2021-10-26T12:54:43Z")
Thanks so much, @stevengj for the very clear answer! Also, thanks for the `dpstrf` pointer, I will go through the LAPACK documentation.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 21, 2024, 7:10pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-low-rank-positive-semidefinite-matrix/70397/6 "2024-09-21T19:10:40Z")
2 posts were split to a new topic: [Sparse Cholesky decomposition of singular semidefinite matrix?](https://discourse.julialang.org/t/sparse-cholesky-decomposition-of-singular-semidefinite-matrix/119682)