# Cholesky Decomposition of a Sparse Symmetric Positive Semidefinite (SPSD) Singular Matrix
**URL:** https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682
**Category:** Numerics
**Tags:** linearalgebra, numerics, sparse, factorization
**Created:** [September 21, 2024, 2:57pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682 "2024-09-21T14:57:40Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 21, 2024, 2:57pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/1 "2024-09-21T14:57:40Z")
> [@Cholesky decomposition of low-rank positive-semidefinite matrix](https://discourse.julialang.org/t/cholesky-decomposition-of-low-rank-positive-semidefinite-matrix/70397/2):
>
> 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).)
Is there such equivalent for sparse matrix?
I don’t see the pivoting option for `SparseMatrixCSC`.
How would you handle a Sparse SPSD matrix?
Currently I just use `check = false`.
Probably small shift should help but it also changes my results too much.
I wonder if there is something clever to do.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 21, 2024, 7:09pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/2 "2024-09-21T19:09:50Z")
> [@RoyiAvital](#):
>
> Probably small shift should help but it also changes my results too much.
What results? What problem are you trying to solve with Cholesky of a singular semidefinite matrix? (If you are solving normal equations arising from a least-squares problem you might want to consider QR instead, especially in an ill-conditioned case.)
PS. You should generally not revive years-old threads with new questions (“necropost”) … instead, start a new thread, and cross-reference/link the old thread as needed.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 22, 2024, 6:01am UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/3 "2024-09-22T06:01:13Z")
The context I encountered it is in [Solve Large Scale Underdetermined Linear Equation with per Element Equality Constraint](https://discourse.julialang.org/t/119655).
The matrix to decompose, \boldsymbol{B}\_{\mathcal{U}}, is SPSD. When I set `check = false` I get the correct result. Yet I wonder if the missing pivot is crucial in this context.
I thought it would be better in the original thread as it is the first that show up when you do a Google search. So it is better to have high Information / Thread ratio (The same problem for the 2 cases: Dense / Sparse).
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 22, 2024, 11:42am UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/4 "2024-09-22T11:42:12Z")
> [@RoyiAvital](#):
>
> The context I encountered it is in [Solve Large Scale Underdetermined Linear Equation with per Element Equality Constraint](https://discourse.julialang.org/t/119655).
Yes, so why not use sparse QR, which should be the default if you just `A \ b`? This should be better behaved anyway since your problem is apparently ill-conditioned?
> [@RoyiAvital](#):
>
> I thought it would be better in the original thread as it is the first that show up when you do a Google search.
The original thread is not about sparse problems, and was already resolved. By introducing a new variation on the original thread you are diluting it and making it harder to follow, as well as confusing the forum by making an old thread seem new.
For the same reason, if you want to open a general discussion of necroposting (which you can easily find many internet discussions of), it would be better to start a new thread.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 22, 2024, 11:58am UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/5 "2024-09-22T11:58:45Z")
> [@stevengj](#):
>
> Yes, so why not use sparse QR, which should be the default if you just `A \ b`? This should be better behaved anyway since your problem is apparently ill-conditioned?
Wouldn’t that be slower than going the Cholesky path?
Indeed the matrix is rank deficient, yet still it is SPSD.
Hence I assumed a matching factorization will be more efficient than a general.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 22, 2024, 4:38pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/6 "2024-09-22T16:38:05Z")
> [@RoyiAvital](#):
>
> Hence I assumed a matching factorization will be more efficient than a general.
First, QR for least-squares _is_ specialized to the problem — you apply QR to A, _not_ to A^T A (which you never compute).
Second, while QR on A is generally slower than forming A^T A and doing Cholesky (the “normal equations” approach), the latter can be _much_ less accurate when A is badly conditioned (i.e. the columns are nearly linearly dependent, i.e. A^T A is nearly singular). If you care about getting the _right_ answer, and not just the _fastest_ answer, I would use QR.
See also [Efficient way of doing linear regression - #33 by stevengj](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/33)
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 22, 2024, 5:22pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/7 "2024-09-22T17:22:16Z")
Indeed, as always all you said is accurate and well explained.
I’m aware about the numerical consequence of the squaring.
I have mentioned them when balancing the 2 alternatives.
In the context of this specific question.
How would you handle the decomposition of rank deficient SPSD sparse matrix?
If I understand correctly, pivoting is not an option?
Does it make sense only use `check = false`?
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 22, 2024, 6:23pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/8 "2024-09-22T18:23:33Z")
> [@RoyiAvital](#):
>
> How would you handle the decomposition of rank deficient SPSD sparse matrix?
> If I understand correctly, pivoting is not an option?
> Does it make sense only use `check = false`?
As I said in the original thread, the problem is that roundoff errors usually make it impossible to distinguish semidefinite matrices from slightly indefinite matrices. I looked through the CHOLMOD documentation and I don’t see any option to specify a tolerance to ignore slightly negative pivots. Pivoting doesn’t help if it’s indefinite. So the problem is that even with `check=false` the Cholesky factorization might simply get stuck partway through.
You could try an L D L^T factorization instead. But even if you have a factorization that succeeds, I wouldn’t use the normal equations in an ill-conditioned case: as soon as you form A^T A you have lost too many digits.
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [September 22, 2024, 11:14pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/9 "2024-09-22T23:14:20Z")
Instead of the QR of A, you could try the QR of `[A ; λI]` for some λ \> 0. Start with `x = 0` and for each k, define `bₖ = [b - A * x ; 0]`. Solve the least-squares problem for `Δx` and update `x = x + Δx`.
intermediate solutions will converge to the minimum-norm solution of the rank-deficient problem. That is the method of Golub and Riley. See, e.g., [1] for more information, meaningful stopping conditions (that depend on the rank), etc.
Of course, you could also just use iterative methods such as LSQR or LSMR as they identify the min-norm solution.
[1] [https://onlinelibrary.wiley.com/doi/abs/10.1002/(SICI)1099-1506(199803/04)5:2\<79::AID-NLA126\>3.0.CO;2-4](https://onlinelibrary.wiley.com/doi/abs/10.1002/(SICI)1099-1506(199803/04)5:2%3C79::AID-NLA126%3E3.0.CO;2-4).
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 22, 2024, 11:35pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/10 "2024-09-22T23:35:54Z")
(There are also direct methods for the minimum-norm solution, e.g. using the LQ factorization.)
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [September 23, 2024, 2:06am UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/11 "2024-09-23T02:06:49Z")
They are the same as for least squares; they compute the QR of the transposed matrix.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 23, 2024, 4:34am UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/12 "2024-09-23T04:34:17Z")
If I solve the system \boldsymbol{A} \boldsymbol{x} = \boldsymbol{b} for \boldsymbol{A} which is rank deficient SPSD matrix.
Assume I use Cholesky decomposition on \boldsymbol{A} and \boldsymbol{b} \notin R \left( \boldsymbol{A} \right).
What solution will I get? Will it be the LS solution?
Let’s assume the matrix is SPSD. Assume it is not coming from the normal equation.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 23, 2024, 12:08pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/13 "2024-09-23T12:08:47Z")
> [@RoyiAvital](#):
>
> What solution will I get? Will it be the LS solution?
You will get garbage in floating-point arithmetic, if the Cholesky factorization even completes. In exact arithmetic, you would get division by zero if b is not in the range (column space) of A.
(I’m assuming “use Cholesky” means that you apply Cholesky elimination steps to both sides and then try to run backsubstitution / triangular solves. This process is seeking an exact solution which does not exist in your example.)
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 23, 2024, 12:23pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/14 "2024-09-23T12:23:18Z")
What I means is the general case where one have an rank deficient SPSD matrix.
In the **dense case** , I’d use the pivot choice as @stevengj described in the linked thread.
In the **sparse case** pivoting is not an option. Hence only doing Cholesky decomposition without any check/
In Julia:
```julia
# `mA` - Some rank deficient SPSD sparse Matrix
oFct = cholesky(mA; check = false);
vX = oFct \ vB;
```
In my code it yielded the correct solution (If I don’t use `check = false` the decomposition fails).
I wonder when will it fail? Assuming the the decomposition did no yield `NaN`.
My first guess was when \boldsymbol{b} \notin R \left( \boldsymbol{A} \right) it might fail.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 23, 2024, 12:35pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/15 "2024-09-23T12:35:52Z")
> [@RoyiAvital](#):
>
> My first guess was when \boldsymbol{b} \notin R \left( \boldsymbol{A} \right) it might fail.
The `\` has to fail in this case — in exact arithmetic, you are seeking a solution that does not exist with singular triangular solves (which will involve dividing by zero). In floating-point arithmetic, who knows what garbage it might produce if the zeros get replaced by some roundoff error.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 23, 2024, 12:47pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/16 "2024-09-23T12:47:12Z")
I see.
Summarizing the case of solving a linear system with rank deficient SPSD matrix:
1. Dense Case: Use the Bunch Kaufman decomposition (Sometimes called LDLt as well, see [Shouldn’t `bunchkaufman()` Be Named `ldlt()`](https://discourse.julialang.org/t/77641)).
2. Sparse Case: Use the LDLt decomposition.
In case \boldsymbol{b} \in R \left( \boldsymbol{A} \right), one may use Cholesky Decomposition with `check = false`. Yet the decomposition is not guaranteed to work still.
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [September 23, 2024, 1:15pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/17 "2024-09-23T13:15:19Z")
I always call Bunch-Kaufman LBLt, because B is block diagonal. That’s in contrast with the “signed Cholesky” / SQD factorization LDLt.
Note that Riley’s method above was initially proposed for your use case: A SPSD. Perform
(A + λI) Δx = b - Ax ; x = x + Δx
repeatedly, and you will get the min-norm solution of Ax = b.
Alternatively, we recently proposed method MINARES [2] that performs particularly well when A is symmetric and singular, and b is not in the range of A. It will be cheaper than LSMR on the least-squares problem.
[2] [[2310.01757] MinAres: An Iterative Solver for Symmetric Linear Systems](https://arxiv.org/abs/2310.01757)
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 23, 2024, 1:37pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/18 "2024-09-23T13:37:01Z")
@dpo , I noticed.
Wouldn’t this iterative approach be slower?
For relative small sparse system with simple structure (In this case, up to 9 elements per row where is row is at the order of ~1e6 elements) I’d guess the direct methods are faster.
I guess I need to set a good pre conditioner.
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [September 23, 2024, 4:46pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/19 "2024-09-23T16:46:54Z")
It’s worth benchmarking. If A is that sparse, matrix-vector products will be very fast.
I just wanted to correct a statement I made above: MINARES will be cheaper than LSMR on Ax=b, not necessarily cheaper than on the LS problem.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 23, 2024, 5:31pm UTC](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/20 "2024-09-23T17:31:44Z")
> [@dpo](#):
>
> It’s worth benchmarking. If A is that sparse, matrix-vector products will be very fast.
[I did benchmark](https://discourse.julialang.org/t/119655/8). They were order of magnitude slower.