Inverse of symmetric tridiagonal block Toeplitz matrix

Question

There is a triagonal block matrix $M$ of form:

$$ M = \begin{bmatrix} A & B^T & 0 & 0 & \cdots & 0 & 0 \\ B & A & B^T & 0 & \cdots & 0 & 0 \\ 0 & B & A & B^T & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & B & A \end{bmatrix} $$ where $A, B$ are real-valued square matrices of the same size.

Also, $A$ is positive definite and symmetric. Later makes $M$ symmetric as well.

My interest is in closed-form solution for elements of $M^{-1}$.

From "Explicit inverses of some tridiagonal matrices" C.M. da Fonseca, J. Petronilho, I am aware of the closed-form solution for tridiagonal toeplitz matrix of form:

$$ T = \begin{bmatrix} a & b & 0 & 0 & \cdots & 0 & 0 \\ c & a & b & 0 & \cdots & 0 & 0 \\ 0 & c & a & b & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & c & a \end{bmatrix} $$ where $a, b, c$ - scalars.

$$ (T^{-1})_{ij} = \begin{cases} (-1)^{i+j}\frac{b^{j-i}}{\left(\sqrt{bc}\right)^{j-i+1}}\frac{U_{i-1}(d)U_{n-j}(d)}{U_{n}(d)} \quad \text{if} \quad i \le j \\ (-1)^{i+j}\frac{c^{i-j}}{\left(\sqrt{bc}\right)^{i-j+1}}\frac{U_{j-1}(d)U_{n-i}(d)}{U_{n}(d)} \quad \text{if} \quad i \gt j \end{cases} $$ where $d = \frac{a}{2\sqrt{bc}}$, $U_{k}(x)$ - Chebyshev polynomials of the second kind.

Don't see a way to extend it to block matrices though. Does anybody know how it can be done? or any alternative way?

Answer 1

After some extended search, found an article "Numerically Stable Algorithms for Inversion of Block Tridiagonal and Banded Matrices" that gives direct algorithm for even more general problem.

$$ M = \begin{bmatrix} A_1 & -B_1 & 0 & 0 & \cdots & 0\\ -B_1^T & A_2 & -B_2 & 0 & \cdots & 0\\ 0 & -B_2^T & A_3 & -B_3 & \cdots & 0\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & A_N \end{bmatrix} $$ where all $A_i$ are symmetric, all $B_i$ are non-singular.

Then, symmetry of $M$, as a whole, is exploited to get decomposition $(M^{-1})_{ij}=U_iV_j^T,\quad j \ge i$.

Note, when inverse of a symmetric matrix exists (and we assume it does) it should be symmetric as well.

$U$ and $V$ are given recursively:

$$ U_1 = I, \quad U_2=B_1^{-1}A_1\\ U_{i+1}=B_i^{-1}(A_iU_i-B^T_{i-1}U_{i-1}),\quad i=2,...,N-1\\ V_N^T=(A_NU_N-B_{N-1}^TU_{N-1}^T)^{-1}, \quad V_{N-1}^T=V_N^TA_NB_{N-1}^{-1}\\ V_i^T=(V_{i+1}^TA_{i+1}-V_{i+2}^TB_{i+1}^T)B_i^{-1}, \quad i=N-2,...,1 $$

Under constraints, that all $A_i$ are equal, all $B_i$ are equal, and $B=B^T$ (which I missed when posting original question!) equations become even simpler:

$$ U_1 = I, \quad U_2 = B^{-1}A \\ U_{i+1} = B^{-1}AU_i - U_{i-1},\quad i=2,...,N-1\\ V_N^T=(AU_N-BU_{N-1}^T)^{-1}, \quad V_{N-1}^T=V_N^TAB^{-1}\\ V_i^T=V_{i+1}^TAB^{-1}-V_{i+2}^T, \quad i=N-2,...,1 $$