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.


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

  • 1
    $\begingroup$ 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. $\endgroup$ May 5, 2020 at 6:40
  • $\begingroup$ Ah I understand now. Thanks, I’ll give it a shot $\endgroup$
    – Justin
    May 5, 2020 at 14:50

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

  • $\begingroup$ 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? $\endgroup$
    – Justin
    May 4, 2020 at 17:06
  • $\begingroup$ No, you don't need to compute $A^{-T}$. You just solve the $2\times 2$ system in my displayed equation. $\endgroup$ May 4, 2020 at 17:25
  • $\begingroup$ Can you check my edited post? I am not sure if my logic w.r.t. solving the $2\times 2$ system works. $\endgroup$
    – Justin
    May 5, 2020 at 3:08
  • $\begingroup$ @Justin Done, I have commented under it. $\endgroup$ May 5, 2020 at 6:41

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.