Solving triangular matrix equations on a GPU

Question

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

Answer 1

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.