"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

Answers related to “more efficient way to invert a matrix knowing it is symmetric and positive semi-definite”

Only authorized users can answer the Search term. Please sign in first, or register a free account.