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

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.

• If a matrix has an eigenvalue $0$, then $\text{det}(A)=0$ and we are dealing with a singular matrix. Solving a system of linear equations with a singular matrix poses problems of its own, independent of whether you use CG or ILU to solve it. – Thijs Steel May 28 '20 at 6:54
• @ThijsSteel Oh, you're right. So for most applications my question should probably be what boundary conditions make a system solvable. But I was still interested in what will actually happen when trying to use a CG solver for such a system. Another follow up question would probably be what kind of solutions the "more general" solvers provide in that case. – allo May 28 '20 at 11:02

## 1 Answer

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 https://arxiv.org/abs/1809.00793 . 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.