Heinz Heinz - 15 days ago 7
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).

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.