I'm inverting covariance matrices with
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
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.