# "more efficient way to invert a matrix knowing it is symmetric and positive semi-definite" Code Answer

3

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.

By nicolaas on May 8 2022