# Solving linear system of the form $ABx=b$

I would like to solve a linear system of the form $ABx=b$ in parallel, where $A$ and $B$ are large, sparse matrices. Currently, I am forming the system matrix $AB$ explicitly, but I would like to avoid the expensive sparse matrix-matrix product. I can think of two alternatives, but I have questions about each of them.

First, I could implement the mapping $x\mapsto ABx$ as a sequence of two sparse matrix-vector products and then apply a Krylov method. While this is easy to set up, I am worried that out-of-the-box preconditioners like (S)SOR or ILU preconditioners may no longer work.

Second, I could solve the system $ABx=b$ by first solving $Ay=b$ followed by $Bx=y$, since both matrices are square. However, I am not sure how to choose the tolerances in this case. It seems that any tolerance $\varepsilon$ for the original system can be achieved by choosing tolerances $\varepsilon_1$ for the system $Ay=b$ and $\varepsilon_2$ for the system $Bx=y$ such that

\begin{align} \|b-ABx\| &=\|b-Ay+Ay-ABx\|\\ &\le\|b-Ax\|+\|A\|\|y-Bx\|\\ &< \varepsilon_1+\|A\|\varepsilon_2\\ &< \varepsilon, \end{align}

but this does not seem practical, because $\|A\|$ is unknown, and $\varepsilon_2$ would likely be very small.

Are you aware of any standard methods to solve this kind of linear system?

• Have you considered using sparse direct factorization methods on $A$ or $B$ or both rather than doing this all iteratively? How big and sparse are your matrices? How badly conditioned are they? – Brian Borchers Jul 24 '18 at 1:57
• The number of unknowns is in the millions or tens of millions for large cases. The number of nonzeros per row should be less than 30 for each of the matrices $A$ and $B$, but I'll have to check that. I do not know the condition numbers, but an SOR-preconditioned GMRES solver usually converges within a few dozen iterations, at least for cases with a few hundred thousand unknowns. Also, the systems are solved within a Newton iteration, so approximations with coarse tolerances are often sufficient. – cthl Jul 24 '18 at 2:29
• If you're willing to use $\| A \|_{1}$ (which is easy to compute with no more work than a single matrix-vector multiplication) to get a bound on the 1-norm of the residual (or just use $\| A \|_{1}$ as an overestimate of $\| A \|_{2}$), then you can implement the second method. – Brian Borchers Jul 24 '18 at 2:50
• That's an interesting idea. I will experiment with the overestimation of $\|A\|_2$ by $\|A\|_1$. – cthl Jul 24 '18 at 3:02
• If you don't get a good enough solution after one round of this, you can always repeat the process- with a good initial solution it should finish the job quickly. – Brian Borchers Jul 24 '18 at 3:44

Defining the auxiliary variable $y=Bx$ yields the following algebraically equivalent expanded system, $$\underbrace{\begin{bmatrix} 0 & A \\ B & -I \end{bmatrix}}_{K} \underbrace{\begin{bmatrix} x \\ y \end{bmatrix}}_{u} = \underbrace{\begin{bmatrix} b \\ 0 \end{bmatrix}}_{f},$$ which you could solve with GMRES or another nonsymmetric Krylov method.
To build a preconditioner, notice that $$\begin{bmatrix}A^{-1} & 0\\ 0 & I\end{bmatrix}\begin{bmatrix} 0 & A \\ B & -I \end{bmatrix} \begin{bmatrix}B^{-1} & 0 \\ 0 & I\end{bmatrix}= \begin{bmatrix}0 & I \\ I & -I\end{bmatrix},$$ which is extremely well-conditioned (condition number: $\phi^2\approx 2.618$). So if you have good preconditioners $\widetilde{A}$ and $\widetilde{B}$ for $A$ and $B$, you can create good left and right preconditioners for the expanded system as follows: $$P_L:=\begin{bmatrix}\widetilde{A} & 0 \\ 0 & I\end{bmatrix}, \quad\quad P_R:=\begin{bmatrix}\widetilde{B} & 0 \\ 0 & I\end{bmatrix}.$$ Then you solve the preconditioned system: $$Ku=f \quad \Leftrightarrow \quad (P_L^{-1} K P_R^{-1})P_R u = P_L^{-1} f.$$
1. Solve $P_L w = f$ for $w$. (apply inverse of left preconditioner)
2. Solve $(P_L^{-1} K P_R^{-1})v = w$ for $v$. (GMRES)
3. Solve $P_R u = v$ for $u$. (apply inverse of right preconditioner)
4. Extract first component of $u$ to get $x$.
• Thank you for the detailed answer. This makes a lot of sense. I will need to see how I can construct the block matrix efficiently from the sparse matrices $A$ and $B$, but it should be possible. As for the preconditioning, I am not sure how difficult it would be to implement the above strategy. Since the original system $ABx=b$ is not too ill-conditioned, do you think a simple SOR or ILU preconditioner would still work reasonably well when applied to the full block system? – cthl Jul 24 '18 at 14:48