# The error propagation in calculating the inverse using a matrix decomposition

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?

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.
• Thank you, this is very helpful. However, 1. (notation) could you add an explanation for $O(u)$ the machine precision, please? It looked like the more precession it have($1000$ compare to $10$) the higher the upper bound for $||KM-I||$, or was it ($1/1000$ compare to $1/100$). 2. For the matrix condition number to be close to $1$, the matrix elements has to be approximately in the same order? Does that mean the SVD could not help in this case, since the original entry ranged in large orders, and that one of the matrix in SVD is probably going to contain a very large and a very small number? Apr 1 at 10:02
• $O(u)$ means a quantity that goes to 0 as a multiple of $u$, which is usually small ($2.2\times 10^{-16}$ for Float64). Usually this hides polynomial factors in the dimension of the matrix; in this case I think the bound can be proved to hold with $\frac{n}{1-nu}u$ in place of $O(u)$. There are many more details in Higham's book Accuracy and Stability of Numerical Algorithms, which has a chapter on inversion. Apr 1 at 13:53
• And no, a matrix with entries of different sizes is not always ill-conditioned. For instance, any small perturbation of (a multiple of) $I$ is always going to have condition number very close to 1, even if it contains entries with very different magnitudes. Apr 1 at 14:07