知道矩阵对称且为正半定数的更有效的矩阵求逆方法 [英] More efficient way to invert a matrix knowing it is symmetric and positive semi-definite

问题描述

我正在使用numpy反转python中的协方差矩阵.协方差矩阵是对称的并且是正半定的.

I'm inverting covariance matrices with numpy in python. Covariance matrices are symmetric and positive semi-definite.

我想知道是否存在一种针对对称正半定矩阵优化的算法，该算法比numpy.linalg.inv()快(当然，是否可以从python轻松访问它的实现！).我无法在numpy.linalg中找到任何内容或在网上搜索.

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.

正如@yixuan所观察到的，正半定矩阵通常不是严格可逆的.我检查了一下我是否得到了正定矩阵，因此我接受了一个对正定有效的答案.无论如何，在LAPACK低级例程中，我发现DSY*例程已针对Simmetric/hermitian矩阵进行了优化，尽管它们似乎在scipy中丢失了(也许只是安装版本的问题).

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).

推荐答案

该API尚不存在，但您可以为其使用低级LAPACK ?POTRI例程族.

The API doesn't exist yet but you can use the low level LAPACK ?POTRI routine family for it.

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输出.如果为零，则意味着无论正定性如何，都成功地求解了该方程式.因为这是低级呼叫，所以不会执行其他检查.

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

如果比较速度

>>>> %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.

So you get a quite nonnegligible speed benefit. If you are working with complex numbers, then you have to use zpotri instead.

问题是您是否需要逆运算.如果您需要计算B⁻¹ * A，则可能不需要，因为solve(B,A)对此更好.

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.

这篇关于知道矩阵对称且为正半定数的更有效的矩阵求逆方法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！