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?
Best Answer
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 $$