Generalized linear algebraic equation solver

Question

$\newcommand{\d}{\vec{d}}$ $\newcommand{\S}{\vec{S}}$

In Mathematica one can easily solve a linear system given by $$A \vec{S} = \vec{d}$$ where $A$ is a matrix, simply by using S = LinearSolve[A,d].

I don't know if this has a name, but I am wondering if there is a nice way to solve a generalized system of linear equations such as this one (this is a specific example but the form of the matrix on the left is unimportant, only that it is a "matrix of matrices"): $$ \begin{bmatrix} 1 & A & 0 & & & 0 \\[8pt] -A & 1 & A & & & \\[1pt] 0 & -A & 1 & \ddots & & \\[1pt] & & \ddots & \ddots & & \\[1pt] & & & & 1 & A \\[6pt] 0 & & & & -A & 1 \end{bmatrix} \begin{bmatrix} \S_1 \\[3pt] \S_2 \\[3pt] \S_3 \\[3pt] \vdots \\[3pt] \S_{N-1} \\[3pt] \S_{N} \end{bmatrix} = \begin{bmatrix} \d_1 \\[3pt] \d_2 \\[3pt] \d_3 \\[3pt] \vdots \\[3pt] \d_{N-1} \\[3pt] \d_{N} \end{bmatrix} $$ where $1$ is the identity matrix and $A$ is an $n\times n$ matrix. The $\S_i$ and $\d_i$ are $n$-component vectors, the $\d_i$ are known and I want to solve for the $\S_i$.

The "brute force" way to solve this is to Flatten the "vectors of vectors" and ArrayFlatten the "matrix of matrices", then use LinearSolve, then use Partition[S,n] to recover the vectors. But I feel like such systems must not be that uncommon and maybe there is some built in method of doing this.

Answer 1

You can build the system matrix as follows

A = RandomReal[{-1, 1}, {2, 2}];
A = A\[Transpose].A;
n = 10;
AA = Plus[
  KroneckerProduct[SparseArray[{Band[{2, 1}] -> -1, Band[{1, 2}] -> 1}, {n, n}], A],
  KroneckerProduct[IdentityMatrix[n, SparseArray], N@IdentityMatrix[2]]
  ]

The right-hand side can be assembled as follows:

d = RandomReal[{-1, 1}, {n, 2}];
S1 = RandomReal[{-1, 1}, {2}];
Sn = RandomReal[{-1, 1}, {2}];
b = Flatten[d];
b[[;; 2]] += A.S1;
b[[-2 ;;]] -= A.Sn;

Solve as usual (this is a usual linear system). I'd suggest a solver specialized for banded matrices:

x = LinearSolve[AA, b, Method -> "Banded"];

In order to convert back to a list of pairs, you can use

Partition[x,2]