Solving triangular matrix equations on a GPU

Suppose I have these two $$N\times N$$ lower triangular banded matrices:

$$A = \begin{bmatrix} a_0 & & \\ a_1 & a_0 & \\ a_2 & a_1 & a_0 \\ a_3 & a_2 & a_1 & a_0 \\ & \ddots & \ddots & \ddots & \ddots \\ & & a_3 & a_2 & a_1 & a_0 \end{bmatrix}$$ $$B = \begin{bmatrix} 1 & & \\ b_1 & 1 & \\ b_2 & b_1 & 1 \\ b_3 & b_2 & b_1 & 1 \\ & \ddots & \ddots & \ddots & \ddots \\ & & b_3 & b_2 & b_1 & 1 \end{bmatrix}$$

Let $$x$$ and $$y$$ be vectors of size $$N$$. Given vector $$x$$, i need to obtain $$y$$ through the following matrix operations:

$$y = \left(B^T\right)^{-1}A^TB^{-1}Ax$$

For the inverse operations, I am using forward and backward substitutions for $$B^{-1}$$ and $$\left(B^T\right)^{-1}$$, respectively. I have several independent $$x$$ vectors that I need to process, and I am doing all of this on a GPU with CUDA. Each CUDA thread is responsible for a single $$x$$ to $$y$$ operation.

However, I want to parallelize this by applying several CUDA threads to operate on a single $$x$$ vector. This lets me use shared memory as "temporary work vectors" and minimize global memory accesses. I also won't ever have a case where vectors of size $$N$$ are too large to fit in shared memory.

The problem is that it appears tricky to parallelize the forward/substitutions, I was wondering if

1) Is there a mathematical trick to simplify the above matrix equations? As in, only having to do one inverse operation instead of two.

Or

2) Is there a better way to invert my two B matrices? Like a more parallel friendly direct solver, or iterative solver, that lets multiple CUDA threads operate on the same vector?

Side note, the coefficient $$b_1$$ will vary somewhere between 1 and 2, whereas $$b_2$$ and $$b_3$$ both vary between 0 and 1. Thus B is not diagonally dominant.

EDIT: following up to Federico's suggestion, if my system can be rewritten as:

$$\begin{bmatrix} 0 & B \\ B^T & A \end{bmatrix} \begin{bmatrix} y \\ z \end{bmatrix} =\begin{bmatrix} b \\ 0 \end{bmatrix}$$

where $$b = Ax$$, and if I invert the $$2\times 2$$ block system, I get the following:

$$\begin{bmatrix} y \\ z \end{bmatrix} = \left.\begin{bmatrix} 0 & B \\ B^T & A^T \end{bmatrix}\right. ^{-1} \begin{bmatrix} b \\ 0 \end{bmatrix}$$ Could the inverted $$2\times 2$$ be rewritten as: $$\left.\begin{bmatrix} 0 & B \\ B^T & A^T \end{bmatrix}\right. ^{-1} = \left(\frac{1}{0\cdot A^T - B^TB}\right) \begin{bmatrix} A^T & -B^T \\ -B & 0 \end{bmatrix}\\ =-\left(B^TB\right)^{-1}\begin{bmatrix} A^T & -B^T \\ -B & 0 \end{bmatrix}$$ Thus the solution $$y$$ would be as simple as: $$y = -\left(B^TB\right)^{-1}A^Tb = -\left(B^TB\right)^{-1}A^TAx$$ Does this make logical sense? It doesn't seem clear to me how the above is equivalent to $$y = \left(B^T\right)^{-1}A^TB^{-1}Ax$$

Edit2: actually the above doesn’t make sense at all. I think what I did only applies if $$A$$ and $$B$$ were scalars. I’m still wondering if it’s possible to do only one inverse/solve operation (even if it means losing the triangular sparsity pattern), thus opening up the possibility of using preconditioned iterative solvers which might be more amenable to GPUs

• As you suspect, the formula for the inverse of a $2\times 2$ system only holds for scalars, not blocks; so your computations are invalid, unfortunately. But you are mis-interpreting my suggestion; don't try to compute that inverse symbolically. Take that $2N\times 2N$ matrix and pass it directly to your solver. That's your "only one inverse/solve operation"; the only catch is that the matrix you are going to pass to your solver is a bit larger. This won't be a huge issue, since the number of nonzeros has a more important role than the size when dealing with sparse matrices. – Federico Poloni May 5 at 6:40
• Ah I understand now. Thanks, I’ll give it a shot – Justin May 5 at 14:50

1 Answer

1) Is there a mathematical trick to simplify the above matrix equations? As in, only having to do one inverse operation instead of two.

Yes: Schur complement formulation. Your system is equivalent to the larger one

$$\begin{bmatrix} 0 & B\\ B^T & A^T \end{bmatrix} \begin{bmatrix} y\\ z \end{bmatrix} = \begin{bmatrix} -b\\0 \end{bmatrix}$$ with $$b=Ax$$ (which you can simply compute) and $$z=-B^{-1}b$$.

The idea is that one constructs a block-$$2\times 2$$ matrix whose Schur complement is the system matrix to be inverted, $$-BA^{-T}B^T$$ (up to a sign).

This matrix is no longer triangular, but I would guess that by reversing the order of unknowns in some of the blocks you can reduce it to a triangular or anti-triangular matrix, which you can then solve by direct substitution. If you are lucky maybe you can end up with a banded block-$$2\times 2$$ matrix. So this answer is not the end of the story; try some more manipulations to get a simpler form.

• Thanks for your answer, i can definitely see how this can be expressed as a Schur Complement. However, isn’t A^-T a dense matrix that needs to be calculated for explicitly? So doesn’t that mean I still need two inverse operations, one for A and then one for the Schur Complement? – Justin May 4 at 17:06
• No, you don't need to compute $A^{-T}$. You just solve the $2\times 2$ system in my displayed equation. – Federico Poloni May 4 at 17:25
• Can you check my edited post? I am not sure if my logic w.r.t. solving the $2\times 2$ system works. – Justin May 5 at 3:08
• @Justin Done, I have commented under it. – Federico Poloni May 5 at 6:41