Solving a set of linear equations with block structure and weak coupling

Question

I have a standard set of linear equations $Ax=b$ where the Hessian matrix $A$ has the special block structure as shown:

$A= \begin{pmatrix} T & U\\ U^T & V \end{pmatrix}$,

$x= \begin{pmatrix} c\\ d \end{pmatrix}$

Where $A$ is square Hessian matrix and each dimension is of size $n = c_n + d_n$. This compound dimension is quite large, typically $10^6-10^{10}$. The gradient vector $b$ has the same block structure as $x$. This set of equations is solved for a full Newton step to an iterative procedure.

One can think of the matrices $T$ and $V$ as two different basis expansions of a given problem with the coupling between these two bases represented by $U$. This set of linear equations is usually solved through CG or an associated method as the matrix $A$ is too large to store.

What is interesting about this is that the "coupling" block $U$ is quite small. Often, the $c$ and $d$ blocks are solved alternatively without consideration of $U$; however, this can lead to somewhat slow convergence.

It would appear that this can be solved in a somewhat decoupled manner where were $Ax$ Hessian-vector guesses can be split as follows:

$Ax_1 = T*c + U*d$

$Ax_2 = U^T*c + V*d$

For example, $Ax_1$ could be updated on every third iteration of $Ax_2$. As the computation of $Ax_1$ is much more costly than $Ax_2$ this would be quite advantageous.

I assume that this kind of system is well described somewhere, but have been unable to google the correct combination of words. Does anyone know what this kind of system of equations is called?

Also, if anyone has a good recommendation for material on approximate Newton methods for ill-conditioned systems that would be quite helpful as well.

---

Full disclosure: I asked the same question here with a 100 point bounty without luck. I decided to try again here as this community has more experience with practical solutions.

Answer 1

Consider the alternating directions method $$ x_{k+1}=T^{-1}(c-Uy_{k}),\qquad y_{k+1}=V^{-1}(d-U^{T}x_{k}). $$ If $\|U\|$ is sufficiently small (in a rigorously defined sense described below) then the sequence $\{x_{k},y_{k}\}$ is guaranteed to converge to a solution.

But conjugate gradients will always converge faster. To see this, define the matrix $P$ as the version of $A$ in which the coupling is eliminated altogether, $$ P=\begin{bmatrix}T & 0\\ 0 & V \end{bmatrix}, $$ and note that the alternating directions method is really just the preconditioned iterations $$ Pz_{k+1}+(A-P)z_{k}=b\iff z_{k+1}=b-P^{-1}(A-P)z_{k}. $$ It is a well-known theorem that preconditioned conjugate gradient (PCG) with $P$ as the preconditioner converges as fast or faster than the corresponding fixed-point iterations. But of course this may not matter if the alternating directions method already convergences very quickly.

Theorem. Define $\alpha=\sigma_{\max}(U)$ and $\beta=\sqrt{\lambda_{\min}(T)\lambda_{\min}(V)}$. If $\alpha<\beta$, then the alternating directions iteration converges within $$ \frac{\beta}{\beta-\alpha}\log\epsilon^{-1}\text{ iterations,} $$ and PCG with $P$ as the preconditioner converges within $$ \sqrt{\frac{\beta}{2(\beta-\alpha)}}\log2\epsilon^{-1}\text{ iterations.} $$

Proof. After preconditioning, the matrix $P^{-1/2}AP^{-1/2}=I+G$ where $G=P^{-1/2}(A-P)P^{-1/2}$ is written explicitly as $$ G=\begin{bmatrix}0 & T^{-1/2}UV^{-1/2}\\ V^{-1/2}U^{T}T^{-1/2} & 0 \end{bmatrix}, $$ and the spectral radius for the iteration matrix is $$\begin{align*} \rho(G) & =\|T^{-1/2}UV^{-1/2}\|\\ & \le\|U\|\sqrt{\|T^{-1}\|\|V^{-1}\|} = \frac{\alpha}{\beta} \end{align*}$$ where we note that if $A$ is positive definite, then $\|A^{-1}\|^{-1}=\lambda_{\min}(A)$. If the spectral radius is strictily less than 1, then the iterated sequence is convergent. Finally, the condition number for $P^{-1/2}AP^{-1/2}$ is equal to $\kappa=(1+\rho(G))/(1-\rho(G))=(\beta+\alpha)/(\beta-\alpha)$. Gradient descent converges in $\frac{1}{2}(\kappa+1)\log\epsilon^{-1}$ iterations, and PCG gradients converges in $\frac{1}{2}\sqrt{\kappa+1}\log2\epsilon^{-1}$ iterations.

---

References

• Shewchuk, Jonathan Richard. "An introduction to the conjugate gradient method without the agonizing pain." (1994).
• Saad, Yousef. Iterative methods for sparse linear systems. Siam, 2003.
• Greenbaum, Anne. Iterative methods for solving linear systems. Vol. 17. Siam, 1997.
• Kelley, C.T. Iterative methods for linear and nonlinear equations. Siam, 1995.