linear algebra - Inverse of a symmetric block tridiagonal matrix

I am aware of existent discussion on the inverse of a block tridiagonal matrix on this website (for example, How to invert a block tridiagonal matrix?) and I have been googling articles about this topic, but I feel I may be interested in a slightly different setting and I cannot tell whether the references I looked so far discuss that, so I'm posting here.

Similar to the link above, I am interested in the last block along the diagonal, the block in $A^{-1}$ corresponding to $D_n$ in $A$. However, the size of the blocks may vary. I do not assume each $D_i$ must be of same size and I assume each $D_i$ is $n_i \times n_i$.

$$A = \begin{bmatrix}
D_1 & A_2^{\top} & & \\
A_2 & D_2 & A_3^{\top} & & & \\

& \ddots & \ddots & \ddots \\
& & A_{n-1} & D_{n-1} & A_n^{\top} \\
& & & A_n & D_n \\
\end{bmatrix}$$

One reference I looked at is https://epubs.siam.org/doi/pdf/10.1137/0613045 and Theorem 3.4 in it gives a general formula when $A$ is proper (i.e. the matrices $A_i$ are nonsingular). However, I am unsure if my setting fits the paper, as it says "block is of order n" (pg.8), and I wonder if the "order" here means $\Theta(n)$. If it actually means equal-size diagonal block, then I wonder if anyone could point some other reference to me for different-size block setting. Thank you!

Answer

For convenience, let $$T_k =
\begin{bmatrix}
D_1 & A_2^{\top} & & \\

A_2 & D_2 & A_3^{\top} & & & \\
& \ddots & \ddots & \ddots \\
& & A_{k-1} & D_{k-1} & A_k^{\top} \\
& & & A_k & D_k \\
\end{bmatrix}$$ for $k = 1,2,\ldots,m$, where I've let $m$ be the total number of diagonal blocks in the original matrix. This is to avoid confusion since the diagonal blocks are of size $n_1 \times n_1, \ldots, n_m \times n_m$. Our goal is to compute $T_m^{-1}$ as efficiently as possible.

Trivially, $T_1 = D_1$, so $T_1^{-1} = D_1^{-1}$, which can be computed in $O(n_1^3)$ operations.

Now, suppose we've already computed $T_{k-1}^{-1}$ and we wish to compute $T_k^{-1}$. We can partition $$T_k = \begin{bmatrix}T_{k-1} & Z_k^T \\ Z_k & D_k \end{bmatrix}$$ where $Z_k = \begin{bmatrix}0 & A_k\end{bmatrix}$. To invert $T_k$, we can apply the block matrix inverse formula to get $$T_k^{-1} = \begin{bmatrix}T_{k-1}^{-1} + T_{k-1}^{-1}Z_k^TS_kZ_kT_{k-1}^{-1} & -T_{k-1}^{-1}Z_k^TS_k \\ -S_kZ_kT_{k-1}^{-1}& S_k \end{bmatrix} \quad \text{where} \quad S_k = (D_k-Z_kT_{k-1}^{-1}Z_k^T)^{-1}.$$

With $T_{k-1}^{-1}$ already computed, we require the following steps:

1. Multiply $Z_k$ by $T_{k-1}^{-1}$ by $Z_k^T$ to get $Z_kT_{k-1}^{-1}Z_k^T$ -- $O(n_{k-1}^2n_k + n_{k-1}n_k^2)$ operations


2. Subtract $Z_kT_{k-1}^{-1}Z_k^T$ from $D_k$ to get $D_k - Z_kT_{k-1}^{-1}Z_k^T$ -- $O(n_k^2)$ operations


3. Invert $D_k - Z_kT_{k-1}^{-1}Z_k^T$ to get $S_k$ -- $O(n_k^3)$


4. Multiply $S_k$ by $Z_k$ to get $S_kZ_k$ -- $O(n_{k-1}n_k^2)$ operations


5. Multiply $Z_k^T$ by $S_k$ to get $Z_k^TS_k$ -- $O(n_{k-1}n_k^2)$ operations


6. Multiply $-S_kZ_k$ by $T_{k-1}^{-1}$ to get $-S_kZ_kT_{k-1}^{-1}$ -- $O(n_k^2(n_1+\cdots+n_{k-1}))$ operations


7. Multiply $T_{k-1}^{-1}$ by $-Z_k^TS_k$ to get $-T_{k-1}^{-1}Z_k^TS_k$ -- $O(n_k^2(n_1+\cdots+n_{k-1}))$ operations



8. Multiply $Z_k^T$ by $S_kZ_kT_{k-1}^{-1}$ to get $Z_k^TS_kZ_kT_{k-1}^{-1}$ -- $O(n_k^2(n_1+\cdots+n_{k-1}))$ operations


9. Multiply $T_{k-1}^{-1}$ by $Z_k^TS_kZ_kT_{k-1}^{-1}$ to get $T_{k-1}^{-1}Z_k^TS_kZ_kT_{k-1}^{-1}$ -- $O(n_k^2(n_1+\cdots+n_{k-1}))$ operations


10. Add $T_{k-1}^{-1}$ and $T_{k-1}^{-1}Z_k^TS_kZ_kT_{k-1}^{-1}$ to get $T_{k-1}^{-1}+T_{k-1}^{-1}Z_k^TS_kZ_kT_{k-1}^{-1}$ -- $O((n_1+\cdots+n_{k-1})^2)$ operations

Note that many of the above steps take advantage of the fact that $Z_k = \begin{bmatrix}0 & A_k\end{bmatrix}$ and $S_kZ_k = \begin{bmatrix}0 & S_kA_k\end{bmatrix}$ are $n_k \times (n_1+\cdots+n_{k-1})$ matrices which have all zeros except for a block of size $n_k \times n_{k-1}$.

If all the blocks are the same size $n_1 = \cdots = n_m = n$, then the total cost of computing $T_k^{-1}$ from $T_{k-1}^{-1}$, $A_k$, and $D_k$ is $O((k-1)n^3+(k-1)^2n^2)$. Thus, the total cost of computing $T_m^{-1}$ recursively is $O(m^2n^3+m^3n^2)$ as opposed to $O(m^3n^3)$ by just direct inversion. If the blocks are not all the same size, it's a bit harder to analyze how much faster the above method is compared to direct inversion. However, I suspect the above method is still faster in many cases.