目前还没有这个API,但您可以使用低层LAPACK?POTRI例程系列来代替。
sp.linalg.lapack.dpotri
的文档字符串如下:
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)
最重要的是info
输出。如果它为零,那意味着它成功地解决了方程,无论正定性如何。由于这是低级调用,不进行其他检查。
>>>> 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
如果您比较速度
>>>> %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
因此,您可以获得相当可观的速度优势。如果您正在使用复数,则必须使用 zpotri
。
问题在于是否需要求逆矩阵。如果您需要计算 B⁻¹ * A
,那么最好使用 solve(B, A)
。