6
$\begingroup$

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?

Improve this question
$\endgroup$
6
  • 3
    $\begingroup$ 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? $\endgroup$
    – Brian Borchers
    Commented Jul 24, 2018 at 1:57
  • $\begingroup$ 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. $\endgroup$
    – cthl
    Commented Jul 24, 2018 at 2:29
  • 1
    $\begingroup$ 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. $\endgroup$
    – Brian Borchers
    Commented Jul 24, 2018 at 2:50
  • $\begingroup$ That's an interesting idea. I will experiment with the overestimation of $\|A\|_2$ by $\|A\|_1$. $\endgroup$
    – cthl
    Commented Jul 24, 2018 at 3:02
  • $\begingroup$ 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. $\endgroup$
    – Brian Borchers
    Commented Jul 24, 2018 at 3:44

1 Answer 1

Reset to default
12
$\begingroup$

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

So, you can proceed in the following steps:

  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$.
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – cthl
    Commented Jul 24, 2018 at 14:48
  • 2
    $\begingroup$ @cthl ILU or SSOR on the full system will probably perform worse because the full system is strongly indefinite. It's worth a try, but I would expect doing different ILUs on A and B, then building from there (as described) would be better. $\endgroup$
    – Nick Alger
    Commented Jul 24, 2018 at 20:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.