I'm having problems understanding why the conjugate gradient method breaks down for singular matrices. I've read a good introduction to intuitively understanding the CG method through visualizing the quadratic form, but I'm not really understanding mathematically why this seems to be. Does anyone have a good mathematical explanation or any sources that explain why conjugate gradient iteration doesn't work for singular matrices?

• In what sense would you want it to work? Typically, a system won't have any solution at all if the coefficient matrix is singular. – Darsh Ranjan Dec 5 '09 at 22:04
• My question was geared more towards any mathematical explanation of why singular matrices don't work. My understanding, from the previously referenced article (pg 5), is that singular matrices have multiple solutions rather than a single solution. – john Dec 5 '09 at 22:30

It's important to distinguish between CG as an exact method, as originally conceived by Hestenes and Stiefel and as an iterative Krylov-type method. As an iterative method, CG works fine for singular matrices provided you don't do too many iterations. For ill-conditioned systems loss of A-orthogonality of the search vectors causes the method to slow down (and produce spurious eigenvalues if you're calculating them, rather than a solution). Also, as Hestenes later pointed out (1975), CG will converge to the pseudo-inverse solution if you start with an initial approximation of zero.

Like any Krylov subspace method CG will compute the components of the solution starting with the largest eigenvalues. In some problems the largest few eigenvalues are sufficient to get a reasonably accurate solution. In this case CG is great. If you really need the small eigenvalue contribution, then you need to re-orthogonalize the search vectors, or partially re-orthogonalize them. This is not hard. In any case actual singularity is no real problem for the iterative method. Stability issues were solved many years ago.

You can even use CG to compute an SVD of a sparse rectangular matrix. Look at SVDPACK by Michael Berry.

• I've uploaded some PDF lecture notes mainly on CG. Included are some numerical examples of working with ill-conditioned systems. E.g., how 5 iterations of CG will give you the first 4 eigenvalues of a 10-th order Hilbert matrix to a part in a thousand. Another is computing the first 40 singular values of a sparse-ill conditioned matrix. mesoscopic.mines.edu/mediawiki/images/f/ff/… These are free and you can do whatever you want with them. – JohnS Aug 4 '13 at 21:54

My question was geared more towards any mathematical explanation of why singular matrices don't work. My understanding, from the previously referenced article (pg 5), is that singular matrices have multiple solutions rather than a single solution.

Unfortunately, that's not so. If $A$ is a singular matrix, then for most $b$, the equation $Ax=b$ has no solution. It has a solution if and only if $b$ is in the column space of $A$ (or orthogonal to the kernel of $A$, if $A$ is symmetric), in which case there are infinitely many solutions. This is basic linear algebra.

Let's say $A$ is symmetric positive semidefinite, so we can at least entertain the thought of conjugate gradients. Conjugate gradients can be interpreted as maximizing the quadratic function $$f(x) = b^Tx-\frac12 x^TAx.$$ If $A$ is singular, then this function typically has no maximum value, as I will demonstrate. Since $A$ is singular, there are vectors $v$ such that $Av = 0$. If there is any such vector $v$ that is not orthogonal to $b$, consider substituting $x=tv$ into $f(x)$: $$f(tv) = t(b^Tv) - \frac12t^2(v^TAv) = t(b^Tv) = ct,$$ where $c$ is the nonzero real number $b^Tv$. Obviously, this is unbounded as a function of $t$, so we see that there is no maximum value if $b$ is not orthogonal to the kernel of $A$.

A search for "conjugate gradient singular matrix" took me to this question. While the answer is obviously given by the responses, the question can be refined: Can CG still give a working algorithm if the matrix is singular, but behaves as a symmetric positive definite form on a (large) subspace?

A standard example is given by the finite element discretization of the Neumann problem on a simply connected domain. The constant functions are both the kernel and the cokernel of the Laplacian. On functions with vanishing mean, the Laplacian is still a positive definite symmetric operator, and we would like to leverage this structure.

This is non-trivial and best our numerical method is derived from a fully analytic setting, because this might provide us the convergence analysis as well. --- This appraoch is for example elaborated in

On the Finite Element Solution of the Pure Neumann Problem
Pavel Bochev and R. B. Lehoucq
SIAM Review
Vol. 47, No. 1 (Mar., 2005), pp. 50-66