Heinz - 10 months ago 60
Python Question

# More efficient way to invert a matrix knowing it is symmetric and positive semi-definite

I'm inverting covariance matrices with

`numpy`
in python. Covariance matrices are symmetric and positive semi-definite.

I wondered if there exists an algorithm optimised for symmetric positive semi-definite matrices, faster than
`numpy.linalg.inv()`
(and of course if an implementation of it is readily accessible from python!). I did not manage to find something in
`numpy.linalg`
or searching the web.

EDIT:

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*`
routines that are optimised for just simmetric/hermitian matrices, although it seems they are missing in
`scipy`
(maybe it is just a matter of installed versions).

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.