The error propagation in calculating the inverse using a matrix decomposition

Question

I have been trying to calculate the matrix inverse of some large matrix with entries ranging by orders of magnitude. I tried to use the matrix decomposition to simplify the computation, where a matrix \begin{equation} (AB) \end{equation} could be decomposed into two matrices, $A$, and $B$. The matrix $A$ is some diagonal matrix of large magnitude i.e. $diag(1,10^5, 10^{10}, 10^{15},10^{20},...)$, which reduced the large magnitudes in matrix $B$. Further, because $A$ is diagonal, its matrix inverse could be computed exactly.

However, this method failed. Because, though $ A^{-1} $ can be computed exactly and $ B^{-1} $ appeared to be simpler, the numerical error in $ B^{-1}A^{-1} $ may be larger than the numerical error in $ (AB)^{-1} $ In this case, since \begin{equation} A^{-1}A=I \end{equation} was exact, the numerical error in
\begin{equation} \max(B^{-1} B) \end{equation} got propagated from $10^{50}$ to $10^{50000}$ in $\max(B^{-1} A^{-1} A B)$, while $\max((AB)^{-1}AB)$ had an error of $10^{50}$.

However, most of the matrix inverse algorithms I looked at involved some types of decomposition, such as the QR decomposition, Singular value decomposition. Thus, the numerical error might propagate especially when the matrix entry ranged in large orders. But why their matrix algorithms were said to be stabler than Gauss–Jordan elimination?

What are other algorithms that can be used to compute for a large matrix with large entries, that do not have such issues?

Answer 1

Irrespective of how you compute an approximate inverse $K\approx M^{-1}$, there is a limit to the (normwise) accuracy up to which $KM \approx I$ can hold: just because of the fact that $K$ and $M$ are approximated by their truncation to the machine precision $u$, the best inequality you can obtain is $$ \|KM - I\| \leq O(u) \|K\|\|M\|. $$ So this explains why your computations fail: your matrices are horribly ill-conditioned. In particular, it is very well possible that you can't get a better error than what you already have obtained for $AB$.

In particular, when can we expect that the accuracy of inverting $B$ and $AB$ are the same? When $A$ is well-conditioned, that is, when the condition number $\kappa(A) = \|A\| \|A^{-1}\|$ is close to 1 (note that it cannot be smaller than 1, since $1=\|I\|= \|AA^{-1}\|\leq \|A\|\|A^{-1}\|$).

This explains why the classical factorizations 'work'. For a QR decomposition, $A=Q$ is orthogonal and has $\kappa(A)=1$. For a PLU decomposition, $L=A$ is lower triangular with $L_{ii}=1$ and $\max |L_{ij}| \leq 1$. This alone is not sufficient by itself to ensure that $\kappa(L)$ is small, but as a matter of fact it is moderate for most matrices: the stability of Gaussian elimination is a tricky problem, but it is extremely rare to encounter counterexamples. Your matrix $A$, on the other hand, has very large condition number, so decomposing a matrix as the product $AB$ can have a detrimental effect on stability.