Heinz - 15 days ago 7

Python Question

I'm inverting covariance matrices with

`numpy`

I wondered if there exists an algorithm optimised for symmetric positive semi-definite matrices, faster than

`numpy.linalg.inv()`

`numpy.linalg`

As observed by @yixuan, positive semi-definite matrices are not in general strictly invertible. I checked that in my case I just got positive definite matrices, so I accepted an answer that works for positive definiteness. Anyway, in the LAPACK low-level routines I found

`DSY*`

`scipy`

Answer

The API doesn't exist yet but you can use the low level LAPACK `?POTRI`

routine family for it.

The docstring of `sp.linalg.lapack.dpotri`

is as follows

```
Docstring:
inv_a,info = dpotri(c,[lower,overwrite_c])
Wrapper for ``dpotri``.
Parameters
----------
c : input rank-2 array('d') with bounds (n,n)
Other Parameters
----------------
overwrite_c : input int, optional
Default: 0
lower : input int, optional
Default: 0
Returns
-------
inv_a : rank-2 array('d') with bounds (n,n) and c storage
info : int
Call signature: sp.linalg.lapack.dpotri(*args, **kwargs)
```

The most important is the `info`

output. If it is zero that means it solved the equation succesfully **regardless of positive definiteness**. Because this is low-level call no other checks are performed.

```
>>>> M = np.random.rand(10,10)
>>>> M = M + M.T
>>>> # Make it pos def
>>>> M += (1.5*np.abs(np.min(np.linalg.eigvals(M))) + 1) * np.eye(10)
>>>> zz , _ = sp.linalg.lapack.dpotrf(M, False, False)
>>>> inv_M , info = sp.linalg.lapack.dpotri(zz)
>>>> # lapack only returns the upper or lower triangular part
>>>> inv_M = np.triu(inv_M) + np.triu(inv_M, k=1).T
```

Also if you compare the speed

```
>>>> %timeit sp.linalg.lapack.dpotrf(M)
The slowest run took 17.86 times longer than the fastest. This could mean that an intermediate result is being cached.
1000000 loops, best of 3: 1.15 µs per loop
>>>> %timeit sp.linalg.lapack.dpotri(M)
The slowest run took 24.09 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 2.08 µs per loop
>>>> III = np.eye(10)
>>>> %timeit sp.linalg.solve(M,III, sym_pos=1,overwrite_b=1)
10000 loops, best of 3: 40.6 µs per loop
```

So you get a quite nonnegligible speed benefit. If you are working with complex numbers, then you have to use `zpotri`

instead.

The question is whether you need the inverse or not. You probably don't if you need to compute `B⁻¹ * A`

because `solve(B,A)`

is better for that.