Safely produce the inverse of a square matrix, such that the output will exactly be the inverse when the input is full rank. When the input is rank deficient, the output will be the Moore-Penrose pseudoinverse. This is intended to be used where a full-rank matrix is expected, but where rank deficient matrices may be encountered due to roundoff or other small issues. It therefore provides the closest thing to a matrix inverse in these cases, whereas a normal inverse would result in a runtime error.

`inv_a = sminv(a);`

We can do this by taking a singular value decomposition of `a`

, inverting
its singular values, transposing it, and reconstructing, as in:

```
[u, s, v] = svd(b);
inv(b) == v * diag(diag(s).^-1) * u.'
```

Of cource, to be safe, we don't want to invert the singular values when they're below some tolerance, so we check for this.

If `a`

is full rank, then `a * sminv(a) == eye(n)`

, where `n`

is the
dimension of `a`

.

If `a`

is rank deficient, then:

` a * sminv(a) == blkdiag(eye(rank(a)), zeros(n - rank(a)))`

Two optional tolerances can be used to determine if the singular values are big enough for inversion: an absolute tolerance and a relative tolerance. For an absolute tolerance of, e.g., 1e-15, no singular value less than 1e-15 will be inverted. For a relative tolerance of 1e-10, no singular value less than the largest singular value times this relative tolerance will be inverted. If both are provided, the maximum tolerance of the two will be used.

Clearly, this matrix doesn't really produce `a / b`

-- that's undefined
where `b`

has null space. However, the intended use of this function is
for safety in critical applications, where `b`

*should* be full rank but
for whatever reason isn't. It keeps the program from crashing with an
error about division by 0.

a | The matrix to invert |
---|---|

abs_tol | Absolute tolerance; singular values smaller than this number won't be inverted. Can be left empty. |

rel_tol | Relative tolerance; singular values smaller than this number times the largest singular value won't be inverted. By default, 1e-15. |

inv_a | The inverse of |
---|

For the normal case, let's create a full rank matrix and show that matrix inversion works.

```
a = randcov(3);
a * sminv(a)
sminv(a) * a
```

```
ans =
1.0000 0.0000 0.0000
-0.0000 1.0000 -0.0000
0.0000 0.0000 1.0000
ans =
1.0000 0.0000 0.0000
-0.0000 1.0000 -0.0000
0.0000 0.0000 1.0000
```

Now make a matrix with null space and show the results.

```
b = randn(3, 2);
a = b * randcov(2) * b.';
```

Show the dimension of `a`

and its rank.

```
size(a, 1)
rank(a)
```

```
ans =
3
ans =
2
```

Show what happens when we “invert” it.

```
a * sminv(a)
sminv(a) * a
```

```
ans =
0.8235 -0.2174 0.3132
-0.2174 0.7323 0.3857
0.3132 0.3857 0.4442
ans =
0.8235 -0.2174 0.3132
-0.2174 0.7323 0.3857
0.3132 0.3857 0.4442
```

Show that this matches MATLAT's `pinv`

:

```
sminv(a)
pinv(a)
```

```
ans =
0.3772 -0.0152 0.2021
-0.0152 0.1319 0.0830
0.2021 0.0830 0.1715
ans =
0.3772 -0.0152 0.2021
-0.0152 0.1319 0.0830
0.2021 0.0830 0.1715
```

Try a normal inverse -- doesn't work.

`inv(a) * a`

```
[Warning: Matrix is singular to working precision.]
[> In ml2jade (line 254)
In enjaden (line 356)
In deploy_web (line 62)]
ans =
NaN NaN Inf
NaN NaN Inf
NaN NaN Inf
```

`*kf`

v1.0.3 January 18th, 2017

©2017 An Uncommon Lab