Safe application of iterative methods on diagonally dominant matrices

Question

Suppose the following linear system is given $$Lx=c,\tag1$$ where $L$ is the weighted Laplacian known to be positive $semi-$definite with a one dimensional null space spanned by $1_n=(1,\dots,1)\in\mathbb{R}^n$, and the translation variance of $x\in\mathbb{R}^{n}$, i.e., $x+a1_n$ does not change the function value (whose derivative is $(1)$). The only positive entries of $L$ are on its diagonal, which is a summation of the absolute values of the negative off-diagonal entries.

I found in one highly cited academic work in its field that, although $L$ is $not~strictly$ diagonally dominant, methods such as Conjugate Gradient, Gauss-Seidl, Jacobi, could still be safely used to solve $(1)$. The rationale is that, because of translation invariance, one is safe to fix one point (eg. remove the first row and column of $L$ and the first entry from $c$ ), thus converting $L$ to a $strictly$ diagonally dominant matrix. Anyway, the original system is solved in the full form of $(1)$, with $L\in\mathbb{R}^{n\times n}$.

Is this assumption correct, and, if so, what are the alternative rationale? I'm trying to understand how the convergence of the methods still hold.

If Jacobi method is convergent with $(1)$, what could one state about the spectral radius $\rho$ of the iteration matrix $D^{-1}(D-L)$, where $D$ is the diagonal matrix with entries of $L$ on its diagonal? Is $\rho(D^{-1}(D-L)\leq1$, thus different from the general convergence guarantees for $\rho(D^{-1}(D-L))<1$? I'm asking this since the eigenvalues of the Laplacian matrix $D^{-1}L$ with ones on the diagonal should be in range $[0, 2]$.

From the original work:

......................................

At each iteration, we compute a new layout (x(t +1), y(t + 1)) by solving the following linear system: $$ L · x(t + 1) = L(x(t),y(t)) · x(t) \\ L · y(t + 1) = L(x(t),y(t)) · y(t) \tag 8$$ Without loss of generality we can fix the location of one of the sensors (utilizing the translation degree of freedom of the localized stress) and obtain a strictly diagonally dominant matrix. Therefore, we can safely use Jacobi iteration for solving (8)

.......................................

In the above, the notion of "iteration" is related to the underlying minimization procedure, and is not to be confused with Jacobi iteration. So, the system is solved by Jacobi (iteratively), and then the solution is bought to the right-hand side of (8), but now for another iteration of the underlying minimization. I hope this clarifies the matter.

Note that I found Which iterative linear solvers converge for positive semidefinite matrices? , but am looking for a more elaborate answer.

Answer 1

The Jacobi iteration can be proved convergent.

The first thing you should make sure is that $c^T \mathbf{1}_n = 0$, which is the condition for existence of solution (I assume $L=L^T$, otherwise you need $c\in (\mathrm{Ker} L^T)^\perp$) because you said $V_0:=\mathrm{Ker} L = \mathrm{span}\{\mathbf{1}_n\}$. We will use the convention that $V_0$ is also the matrix with columns being orthonormal basis of it. In your case, $V_0:=\mathbf{1}_n/\sqrt{n}$.

Then, for the errors of the Jacobi iteration on the original system, you have $$ e_1 = (I-D^{-1}L)e_0 = (I-D^{-1}L) (P e_0 + V_0a)=(I-D^{-1}L) P e_0 + V_0a,$$ where $P:=I-V_0V_0'$ is the orthogonal projection onto $V_1:=V_0^\perp$. From the above iteration, we know that $$ P e_1 = P (I-D^{-1}L) P e_0, $$
from which we have the iteration matrix $S$ in $V_1$, $$ S: = P (I-D^{-1}L) P. $$ Not that $S$ has the same spectra (except zeros) with the following matrix $$ \tilde{S}:= (I-D^{-1}L) P P=(I-D^{-1}L) P=(I-D^{-1}L)(I-V_0V_0')\\ =I-D^{-1}L-V_0V_0'.$$ We want the spectral radius of $S$ less than one to prove the convergence.

The following quote is old and kept only for reference. See after for the new proof.

In your case, $V_0V_0'=\frac{1}{n}\mathbf{1}_{n\times n}.$ And you can verify that $D^{-1}L+V_0V_0'$ is strictly diagonal-dominant by using the assumption that the entries of $L$ are positive on the diagonal and negative otherwise. To show the eigenvalues of $D^{-1}L+V_0V_0'$ are real, we note that the matrix is self-adjoint under the inner product $<x,y>:=y^TDx.$

If $V_0$ is not in your specific form, I have not found an answer to the convergence question. Could someone clarifies this?

Note that $V_0$ is the eigen-vector corresponding to the eigenvalue $1$ of $I-D^{-1}L$. Based on the observation, we call Theorem 2.1 from Eigenvalues of rank-one updated matrices with some applications by Jiu Ding and Ai-Hui Zhou.

Theorem 2.1 Let $u$ and $v$ be two $n$-dimensional column vectors such that $u$ is an eigenvector of $A$ associated with eigenvalue $\lambda_1$. Then, the eigenvalues of $A+uv^T$ are $\{\lambda_1+u^Tv,\lambda_2,\ldots,\lambda_n\}$ counting algebraic multiplicity.

Then, we know that the spectra of $\tilde{S}$ is the same as $I-D^{-1}L$ except that the eigenvalue $1$ in the latter is shifted by $-1$ into the eigenvalue zero in the former. Since $\rho(I-D^{-1}L)\subset (-1,1]$, we have $\rho(\tilde{S})\subset (-1,1)$.

Answer 2

Krylov methods never explicitly use the dimensionality of the space they iterate in, therefore you can run them on singular systems so long as you keep the iterates in the non-null subspace. This is normally done by projecting out the null space at each iteration. There are two things that can go wrong, the first is much more common than the second.

1. The preconditioner is unstable when applied to the singular operator. Direct solvers and incomplete factorization may have this property. As a practical matter, we just choose different preconditioners, but there are more principled ways to design preconditioners for singular systems, e.g. Zhang (2010).
2. At some iteration, $x$ is in the non-null subspace, but $A x$ lives entirely in the null space. This is only possible with nonsymmetric matrices. Unmodified GMRES breaks down in this scenario, but see Reichel and Ye (2005) for breakdown free variants.

To solve singular systems using PETSc, see KSPSetNullSpace(). Most methods and preconditioners can solve singular systems. In practice, the small null space for PDEs with Neumann boundary conditions are almost never a problem as long as you inform the Krylov solver of the null space and choose a reasonable preconditioner.

From the comments, it sounds like you are specifically interested in Jacobi. (Why? Jacobi is useful as a multigrid smoother, but there are much better methods to use as solvers.) Jacobi applied to $A x = b$ does not converge when the vector $b$ has a component in the null space of $A$, however, the part of the solution orthogonal to the null space does converge, so if you project the null space out of each iterate, it converges. Alternatively, if you choose a consistent $b$ and initial guess, the iterates (in exact arithmetic) do not accumulate components in the null space.