# Safe application of iterative methods on diagonally dominant matrices

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.

• Could you post a link or a citation to the highly cited work? – Geoff Oxberry Mar 4 '12 at 5:12
• It could be retreived from: citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.164.1421 Since you are not expected to read the whole work, take a look at p.7 (bottom). I suppose the choice of iterative solvers is justified, but I feel a better (or, at least, different) rationale is needed. – usero Mar 4 '12 at 14:22
• I wonder whether these guys are from the same community as the combinatorial preconditioners. – shuhalo Mar 5 '12 at 14:26

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

• Thanks for answering. Something similar was what I've considered: namely, with the weighted Laplacian structured as the $D^{-1}L$ above, it could be shown that its eigenvalues are within $[0, 2)$, hence with the spectral radius within $(0, 2)$ (one eigenvalue is greater than $0$, and at least one is $0$). Therefore, the spectral radius of the iteration matrix $I-D^{-1}L$ is less then $1$, hence with convergent Jacobi. Perhaps the above assumption on the spectral radius of $I-D^{-1}L$ (excluding $0$) is not safe? – usero Mar 10 '12 at 11:51
• I think the spectra of $D^{-1}L$ should be in $[0,2]$, that is closed at $2$. I do not know how you can get $2$ excluded. From my point of view, the (Gershgorin circle theorem) [en.wikipedia.org/wiki/Gershgorin_circle_theorem] can only give the estimate including $2$. If it is the case, the estimate of the spectral radius of $I-D^{-1}L$ is $\leq 1$ with the equality achievable with the vectors in the kernel of $L$. I think the convergence you want is that in the orthogonal complement space $V_1$ as noted in the above 'answer'. – Hui Zhang Mar 10 '12 at 12:07
• You could take a look at Lemma 1.7 (v) of math.ucsd.edu/~fan/research/cb/ch1.pdf The matrix $D^{−1}L$ could be regarded as a weighted Laplacian on a complete graph, hence with excluded $2$. I guess it is a sufficient argument for the convergence proof?...........Does your approach require other pre/post-processing of the iterates beyond centering $c$. I'm asking because you introduced $V_0$ And regarding the spectra of $I-D^{-1}L-V_0V_0'$: given that the spectral radius ($sr$) of $I-D^{-1}L$ is $(0, 1]$, the addition of $-\frac{1}{n}$, would yield $sr<1$. Isn't this a good enough argument? – usero Mar 12 '12 at 8:26
• Hi, thanks for pointing to a good book. But I found I can not take a quick look. About your last argument, it apppears almost the same as the "answer" above. Just be careful, you are not adding $\frac{1}{n}$ but $\frac{1}{n}\mathbf{1}_{n\times n}$, so it is not a simple addition to the $sr$ of $I-D^{-1}L$. Generally, the $sr$ of sum of two matrices are not simple sum of the $sr$'s of the individual matrices. – Hui Zhang Mar 12 '12 at 10:05
• Good that you pointed that out. Does your approach require other pre/post-processing of the iterates beyond centering c. I'm asking because you introduced $V_0$, and I thought that you're talking about projecting out the null-space. If so, is projection the null-space out really necessary for convergence? – usero Mar 12 '12 at 10:27

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.

• You could perform an orthogonal change of basis so that there is a zero on the diagonal (find any orthogonal matrix $Q$ in which the first column is the constant vector). Under this transformation $A_1 = Q^T A Q$, the matrix $A_1$ is still symmetric positive semi-definite, but the first diagonal entry is 0 so direct application of Jacobi would fail. Since $A_1$ is dense, you wouldn't do this in practice, but this shows that the basis matters. If $Z$ is an orthogonal basis for the null space, projected GMRES is just solving $(I-Z)P^{-1} A x = (I-Z)P^{-1} b$. – Jed Brown Mar 5 '12 at 15:45
• Hmm, it seems I replied to a comment that was deleted. I'll leave the comment here in case it's useful. – Jed Brown Mar 5 '12 at 15:46
• Thanks for the answer, it's on much higher specialized level then I expected. Therefore, I'll need some guides on: 1) how to project out the null space at each iteration? 2) In my understanding, you stated that the Jacobi application to the system as stated primarily might not converge to the exact solution (i.e. the iterands are not getting better solution estimates). It is therefore suggested to choose different preconditioners? If so, does that practically imply a dynamic check on behaviour with $diag(A)$, and change if problem occurs (with the above case of the linear system)? – usero Mar 5 '12 at 15:50
• My 1) from above should be regarded as: given the Jacobi iteration with the system primarily posted, is it needed to project out the nullspace, and, if so, how could one incorporate it within the update $X_{k+1}=D^{-1}(b-(A-D)X_k)$? Postprocessing the iterate $X_{k+1}$, and considering the postprocessed version for $X_{k}$? – usero Mar 5 '12 at 15:59
• In a reasonable basis, Jacobi should be stable. It can also use 1 on the diagonal if the diagonal matrix element is 0, the projection still removes the null space. Are you planning to use a Krylov method like CG or GMRES? If not, why not? If you are, you just need an orthogonal basis for the null space. You only have the constant mode in your null space, so an orthogonal projector into the null space is $N=ZZ^T$ where $Z$ is the column vector. The orthogonal projector that removes the null space is thus $I-N$. (My first comment had a mistake, if $Z$ is the basis, $N=I-ZZ^T$ is the projector.) – Jed Brown Mar 5 '12 at 16:04