# Matrix inversion accuracy

# Theory

Inversion of large matrices is difficult. Newman (1974) proposes a method for testing inversion accuracy.

Let \(A\) be the (square, non-singular) matrix to be inverted. \(X\) is the approximation of \(A^{-1}\) obtained by calculation. \(N(Y)={sqrt{sum_{i,j}{y_{i,j}^2}}\) is the Euclidean norm of a matrix \(Y=(y_{i,j})\).

\(R = I - AX\) is the difference (ideally zero) between the identity matrix and the product of \(A\) and its calculated inverse.

The norm of \(R\) is not a good indicator of inversion accuracy (Newman, 1974, section 3). On the other hand, the norm of \(A^{-1} - X\) is, but \(A^{-1}\) is unknown. Newman shows that \[\frac{N(XR)}{1 + N(XR)} \leq N(A^{-1} - X) \leq \frac{N(XR)}{1 - N(XR)},\] which allows to bound the value of \(N(A^{-1} - X)\) between the calculated values of \(X\) and \(R\).

The aim is for this value to be small, ideally zero if \(A^{-1} = X\), relative to the norm of \(X\).

# Example

The `size`

parameter sets the size (number of rows and columns) of the matrices used, here a 1000 x 1000 matrix.

`size <- 1000`

A matrix is created:

`A <- matrix(rnorm(size^2), nrow = size)`

It is inverted numerically (\(X\) is calculated), then \(R\) and \(XR\) are calculated:

```
X <- solve(A)
AX <- A %*% X
R <- diag(size) - AX
XR <- X %*% R
```

The Euclidean norm function is defined:

```
norm.matrix <- function(X) {
sqrt(sum(X^2))
}
```

Matrix norms are calculated:

```
# Matrix to invert
norm.matrix(A)
```

`## [1] 1000.567`

```
# Numerical inversion
norm.matrix(X)
```

`## [1] 166.33`

```
# Inversion deviation
norm_R <- norm.matrix(R)
norm_XR <- norm.matrix(XR)
```

Finally, the bounds of the norm of \(A^{-1} - X\) is provided:

`norm_XR / (1 + norm_R)`

`## [1] 9.210669e-11`

`norm_XR / (1 - norm_R)`

`## [1] 9.210669e-11`

The norm of \(A^{-1} - X\) is very close to 0: matrix inversion in R works well.

# Reference

Newman, M. (1974). How to Determine the Accuracy of the Output of a Matrix Inversion Program. *JOURNAL OF RESEARCH of the National Bureau of Standards - B. Mathematical Sciences*, 788(2), 65–68.