• Solving Ax=b
• When A is non-singular
• When A is singular
• Can we use LDLT to check whether a P.S.D. matrix is (numerically) singular or not?
• However Eigen library seems to say no

I was recently caught up with a seemingly classic matrix decomposition problem as follows.

Solving Ax=b

To solve a system of equations $ Ax=b $, which is usually overdetermined but sometimes can be underdetermined, one could usually equivalently solve the corresponding normal equation $ Bx=c $, where $B=A^\intercal A$ and $c=A^\intercal b$. Now we are sure that matrix $B$ is positive semidefinite, i.e., P.S.D.. But we are not sure if $B$ is singular or not, since its singularity totally depends on the singularity of $A$.

When A is non-singular

we know that the system can be efficiently solved by perform Cholyskey Decomposition as $B=LL^\intercal$. In the widely-used Eigen library, this can be done as simple as:

Eigen::MatrixXd x = B.llt().solve(c);

or a more stable version as the LDLT decomposition¹:

Eigen::MatrixXd x = B.ldlt().solve(c);

When A is singular

$B$ becomes a singular P.S.D. Now the question is, can we perform Cholyskey Decomposition on a singular P.S.D.? The answer is yes! As explained in this book (equation 2.78 to 2.85), if we always set the indeterminate variable to 0, we can still get a valid decomposition result. However as explained in section 2.7.4 of that book:

The Cholesky algorithm is unstable for singular positive semidefinite matrices h.
It is also unstable for positive definite matrices h that have one or more eigenvalues close to 0.

So the previous answer should be augmented as “yes, but the decomposition is NOT numerically stable”. But wait a second, who cares about whether the decomposition is stable or not? As long as it can be reliably determined that matrix $B$ is singular, we can simply return that this system is underdetermined and thus fail to find a unique solution. However we really want to avoid using SVD or QR decomposition to check the rank of $B$ (for speed consideration, if we need to repeatedly solve a large number of such kind of equations).

Since we are going to solve the system using LDLT if it is non-singular, now the question becomes:

Can we use LDLT to check whether a P.S.D. matrix is (numerically) singular or not?

If we can, then we can directly perform LDLT on $B$ and then examine the singularity (and return failure if it is indeed singular). My answer is yes and my thoughts are explained as follows. We know there always exists a LDLT decomposition for a given P.S.D. matrix $B$ as $B=LDL^\intercal$. Considering the definition of $L$ and $D$, we know that $det(B)=det(D)=\prod d_i$, where $D=diag(d_1, d_2, …, d_n)$, we can then conclude that:

1. if $B$ is non-singular, then none of the diagnal elements $d_i$ is zero,
   since otherwise the determinant of $B$ will become zero.
2. and similarly, if $B$ is singular, then there is at least one diagnal element $d_i \approx 0$ numerically,
   since otherwise the determinant of $B$ will become non-zero.

Thus in the Eigen library all we need to do is to check:

bool is_B_singular = B.ldlt().vectorD().cwiseAbs().minCoeff() < 1e-8; // if min(abs(d_i))>eps

Or a more numerically stable check:

Eigen::VectorXd D = B.ldlt().vectorD().cwiseAbs();
bool is_B_singular = D.maxCoeff() > (D.minCoeff() * 1e8); // similar to check matrix condition number in some sense

However Eigen library seems to say no

In its documentation of LDLT, it says that:

Remember that Cholesky decompositions are not rank-revealing.
Also, do not use a Cholesky decomposition to determine whether a system of equations has a solution.

OK, I agree that Cholesky decompositions are not rank-revealing (I cannot find a textbook stating this, but I think it relates to numerical stability of the decomposition. So if a matrix is singular, we may not be able to simply count the number of non-zero elements as the rank of the matrix). But why can’t I use it to determine whether the previous $Bx=c$ has a unique solution or not by simply looking at the singularity or $B$? Or maybe the documentation only means that I cannot use the test of $rank([B,c]) == rank(B)$ as mentioned in this Matlab post to test if the system still has solution even if $B$ is singular? (i.e., “$B$ is singular” is not equivalent to “$Bx=c$ has no solution”) If the later is what the documentation means, then Eigen’s domentation is really not saying no to our approach. This is the part that I’m not sure yet and need to ask Eigen’s develop team for confirmation, althought according to my own experiments with Eigen it seems to work pretty well using the above code.

OK, here is what Eigen’s expert says:

Such a test is indeed not reliable, it will only warn you that the problem might not be reliably addressed by a LDLT decomposition and that using a rank revealing QR of SVD might be safer.

So I guess it still make sense to use the above code to check…

Updated:

According to the ceres solver’s documentation:

Cholesky factorization is not a rank-revealing factorization, i.e., it cannot reliably detect when the matrix being factorized is not of full rank.

I think it is NOT OK to use the above code to check…