MathOverflow will be down for maintenance for approximately 3 hours, starting Monday evening (06/24/2013) at approximately 9:00 PM Eastern time (UTC-4).

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 2009 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 2009 at 22:30

If the matrix is singular then then the images of some vectors may be zero and that could mean some of the steps may end up with division by zero.

-

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$.

-

The conjugate gradient method becomes unstable when the matrix A is singular, i.e., if you compute the output of the k-th iteration for a small change in the initial data, large deviations will occur. A similar large deviation occurs due the computer roundoff error. Even for non singular matrices A that are close to to be singular, for example when one of the eigenvalues of A is smaller by orders of magnitudes than the rest, instability problems can occur: The convergence rate becomes slower, and even convergence can be lost due to roundoff errors.

-

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