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.
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.
The context I encountered it is in Solve Large Scale Underdetermined Linear Equation with per Element Equality Constraint.
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).
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?
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.
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.
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
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?
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.
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.
(There are also direct methods for the minimum-norm solution, e.g. using the LQ factorization.)
They are the same as for least squares; they compute the QR of the transposed matrix.
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.
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.)
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:
# `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.
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.
I see.
Summarizing the case of solving a linear system with rank deficient SPSD matrix:
bunchkaufman() Be Named ldlt()).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.
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
@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.
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.
I did benchmark. They were order of magnitude slower.