What happens when I use a conjugate gradient solver with a symmetric positive semi-definite matrix?

Question

I have a symmetric positive semi-definite matrix, i.e., a laplacian and wonder what may happen when I use a CG solver, that is an algorithm for positive definite matrices.

• What happens when the matrix is only positive semi-definite?

• Can I avoid problems by somehow substracting the eigenvector for the eigenvalue 0?

For the laplacian I know the eigenvector(s) for eigenvalue 0, because they are the indicator vectors for the connected components of the domain.

I know there are some numeric tricks like adding a small value to the main diagonal, but I guess this will either make the solution worse when the value is too large or the algorithm will still be instable when the value is too small.

Answer 1

You first have to make sure that your system is solvable: this happens iff the right-hand side $b$ is orthogonal to the kernel of $A$. If $A$ has a dimension-1 kernel spanned by $v$, you need to have $v^* b=0$. If that is not the case, then go back to the modelling stage and ask yourself if what you did makes sense.

Convergence of iterative methods for singular systems is tricky in general, but for positive-semidefinite matrices and CG all works nicely: CG will converge on your singular system without any modifications. More precisely, it will converge to the unique solution orthogonal to the kernel ($v^Tx=0$), which is also the minimum-norm one.

For a proof, see Convergence of the Conjugate Gradient Method on Singular Systems. Intuitively, what happens is that, after a suitable orthogonal change of basis, you are solving the system $$ \begin{bmatrix} A_{11} & 0\\ 0 & 0 \end{bmatrix} \begin{bmatrix} x_1\\ x_2 \end{bmatrix} = \begin{bmatrix} b_1\\0 \end{bmatrix}, $$ with $A_{11}$ positive definite. In this form, it is not complicated to see that all the iterates that CG produces have a zero in their second block: $(x_k)_2=0$, for each $k$. In addition, the iterates coincide with the ones that you'd get by applying CG to the (nonsingular) linear system $A_{11}x_1=b_1$.

Numerically, you may want to re-orthogonalize $x_k \leftarrow x_k - vv^Tx_k$ every few steps to make sure that your iterates stay orthogonal to the kernel.

Answer 2

To at least partially answer the question:

You have to make sure the system is solvable in the first place. As mentined in the accepted answer, "the system $Ax=b$ is solvable iff $b \bot Ker(A)$", but that was never obvious to me.

The obvious condition of solvability is $b \in Im(A)$. Next, it holds, that $Im(A)$ and $Ker(A^\top)$ are each other's orthogonal complements (for any matrix $A$). And since $A$ is symmetric in our case, $Ker(A^\top) = Ker(A)$.