I have implemented a blockwise matrix inversion. When comparing the blockwise inversion to an inverse of the entire matrix I am seeing deviations from the correct inverse. What should the expectation for blockwise inversion be regarding stability relative to a standard inverse? This causes an issue because multiplying the blockwise results by the original matrix no longer results in the identity.
Matrix being inverted:
[0.02788, 0.00679, 0.00425, 0.00515
0.00679, 0.14084, 0.00497, 0.00289
0.00425, 0.00497, 0.03055, 0.005
0.00515, 0.00289, 0.005, 0.05109]
Inverse Result:
[37.5317721, -1.5862237, -4.4296639, -3.2600532,
-1.5862237, 7.2126689, -0.9269588, -0.1573844,
-4.4296639, -0.9269588, 33.962595, -2.8248443,
-3.2600532, -0.1573844, -2.8248443, 20.1872839]
Blockwise Inverse Result:
[37.4861497, -1.5808969, -4.4373105, -3.1229364,
-1.5808969, 7.2152228, -0.9702984, -0.1593941,
-4.4373105, -0.9702984, 33.9787159, -2.8388331,
-3.1229364, -0.1593941, -2.8388331, 20.1639422]
My guess is that this block formula won't work very well in general, just because of how many different matrix multiplies and inversions are going on.
As a sample bad case, consider the matrix, $$ X = \begin{bmatrix} 1 & 0 & 1 & 1 \\ 0 & \epsilon & 1 & 1 \\ 1 & 1 & \epsilon & 1 \\ 1 & 1 & 1 & -\epsilon \end{bmatrix} $$
This matrix is invertible as $\epsilon \to 0$, and in this limit the condition number approaches something near $5.3$. So for small epsilon, the original matrix is fairly well conditioned and stable algorithms should be able to give a good solution.
However, The condition number of the condition number of the $(1,1)$ block is $1/\epsilon$ which goes to infinity as $\epsilon \to 0$. This means the inverse of the $(1,1$) block may not even be representable in finite precision arithmetic, and so the block formula will not work well.