知道矩阵是对称的和半正定的矩阵求逆的更有效方法

Question

我numpy在 python 中反转协方差矩阵。协方差矩阵是对称的和半正定的。

我想知道是否存在针对对称正半定矩阵优化的算法，速度比numpy.linalg.inv()（当然，如果它的实现可以很容易地从 python 访问！）。我没有设法numpy.linalg在网上找到或搜索一些东西。

编辑：

正如@yixuan 所观察到的，半正定矩阵通常不是严格可逆的。我检查了在我的情况下我只是得到了正定矩阵，所以我接受了一个适用于正定性的答案。无论如何，在 LAPACK 低级例程中，我发现DSY*仅针对 simmetric/hermitian 矩阵进行了优化的例程，尽管它们似乎缺失scipy（也许只是安装版本的问题）。

Answer 1

首先，您需要确保协方差矩阵是正定矩阵（pd）而不是半定矩阵，否则矩阵不可逆。

在没有pd假设的情况下，矩阵求逆通常通过LU分解来完成，而对于pd矩阵，可以使用Cholesky分解，这通常会降低计算成本。

在 Scipy 中，该linalg.solve()函数有一个参数sym_pos，假设矩阵为 pd。下面是一个简单的例子：

import numpy as np
from scipy import linalg
import time

def pd_inv(a):
    n = a.shape[0]
    I = np.identity(n)
    return linalg.solve(a, I, sym_pos = True, overwrite_b = True)

n = 50
A = np.random.rand(n, n)
B = A.dot(A.T)

start = time.clock()
for i in xrange(0, 10000):
    res = linalg.inv(B)

end = time.clock()
print end - start
## 1.324752

start = time.clock()
for i in xrange(0, 10000):
    res = pd_inv(B)

end = time.clock()
print end - start
## 1.109778

在我的机器上，pd_inv()对于小矩阵（~100x100）有一些优势。对于大型矩阵，几乎没有任何区别，有时pd_inv()甚至更慢。

Answer 2

据我所知，对称矩阵没有标准的矩阵逆函数。通常，您需要对稀疏性等进行更多约束，才能为您的求解器获得良好的加速。但是，如果您查看scipy.linalg，您会发现有一些针对 Hermitian（对称）矩阵进行了优化的特征值例程。

例如，当我生成一个随机的 200x200 密集矩阵并求解特征值时，我得到：

from scipy.linalg import inv,pinvh,eig,eigh
B = np.rand(200,200)
B = B+B.T
%timeit inv(B)
1000 loops, best of 3: 915 µs per loop

%timeit pinvh(B)
100 loops, best of 3: 6.93 ms per loop

所以相反没有优势，但是：

%timeit eig(B)
10 loops, best of 3: 39.1 ms per loop

%timeit eigh(B)
100 loops, best of 3: 4.9 ms per loop

特征值的 8 倍加速。

如果您的矩阵是稀疏的，您应该查看 scipy.sparse.linalg，它有一些求解器，其中一些（如 bicg 和 cg）需要 Hermitian 矩阵，因此可能会更快。但是，只有当您的矩阵稀疏时才有意义，仅求解特定的解向量b并且实际上可能不会更快，具体取决于矩阵结构。您真的必须对其进行基准测试才能找到答案。

我问了一个关于C++ 求解器的类似问题，最终发现它真的取决于应用程序，你必须为你的问题选择最好的求解器。

Answer 3

该 API 尚不存在，但您可以使用低级 LAPACK?POTRI例程系列。

\n\n

的文档字符串sp.linalg.lapack.dpotri如下

\n\n

Docstring:     \ninv_a,info = dpotri(c,[lower,overwrite_c])\n\nWrapper for ``dpotri``.\n\nParameters\n----------\nc : input rank-2 array('d') with bounds (n,n)\n\nOther Parameters\n----------------\noverwrite_c : input int, optional\n    Default: 0\nlower : input int, optional\n    Default: 0\n\nReturns\n-------\ninv_a : rank-2 array('d') with bounds (n,n) and c storage\ninfo : int\nCall signature: sp.linalg.lapack.dpotri(*args, **kwargs)\n

\n\n

最重要的是info输出。如果它为零，则意味着无论正定性如何，它都成功求解了方程。因为这是低级调用，所以不执行其他检查。

\n\n

>>>> M = np.random.rand(10,10)\n>>>> M = M + M.T\n>>>> # Make it pos def\n>>>> M += (1.5*np.abs(np.min(np.linalg.eigvals(M))) + 1) * np.eye(10)\n>>>> zz , _ = sp.linalg.lapack.dpotrf(M, False, False)\n>>>> inv_M , info = sp.linalg.lapack.dpotri(zz)\n>>>> # lapack only returns the upper or lower triangular part \n>>>> inv_M = np.triu(inv_M) + np.triu(inv_M, k=1).T\n

\n\n

另外如果你比较速度

\n\n

>>>> %timeit sp.linalg.lapack.dpotrf(M)\nThe slowest run took 17.86 times longer than the fastest. This could mean that an intermediate result is being cached.\n1000000 loops, best of 3: 1.15 \xc2\xb5s per loop\n\n>>>> %timeit sp.linalg.lapack.dpotri(M)\nThe slowest run took 24.09 times longer than the fastest. This could mean that an intermediate result is being cached.\n100000 loops, best of 3: 2.08 \xc2\xb5s per loop\n\n>>>> III = np.eye(10)\n\n>>>> %timeit sp.linalg.solve(M,III, sym_pos=1,overwrite_b=1)\n10000 loops, best of 3: 40.6 \xc2\xb5s per loop\n

\n\n

因此，您将获得相当不可忽视的速度优势。如果您正在处理复数，那么您必须使用zpotri。

\n\n

问题是你是否需要逆。如果您需要计算，您可能不会这样做，B\xe2\x81\xbb\xc2\xb9 * A因为solve(B,A)这样做更好。

\n

Answer 4

我尝试了@percusse 的答案，但是当我对其执行计时时，我发现它比np.linalg.inv（使用 100K 随机正定 4x4np.float64矩阵的样本）慢了大约 33% 。这是我的实现：

from scipy.linalg import lapack

def upper_triangular_to_symmetric(ut):
    ut += np.triu(ut, k=1).T

def fast_positive_definite_inverse(m):
    cholesky, info = lapack.dpotrf(m)
    if info != 0:
        raise ValueError('dpotrf failed on input {}'.format(m))
    inv, info = lapack.dpotri(cholesky)
    if info != 0:
        raise ValueError('dpotri failed on input {}'.format(cholesky))
    upper_triangular_to_symmetric(inv)
    return inv

我尝试分析它，令我惊讶的是，它花费了大约 82% 的时间在调用upper_triangular_to_symmetric（这不是“困难”部分）！我认为发生这种情况是因为它正在执行浮点加法以组合矩阵，而不是简单的副本。

我尝试了一个upper_triangular_to_symmetric大约快 87% 的实现（请参阅此问答）：

from scipy.linalg import lapack

inds_cache = {}

def upper_triangular_to_symmetric(ut):
    n = ut.shape[0]
    try:
        inds = inds_cache[n]
    except KeyError:
        inds = np.tri(n, k=-1, dtype=np.bool)
        inds_cache[n] = inds
    ut[inds] = ut.T[inds]


def fast_positive_definite_inverse(m):
    cholesky, info = lapack.dpotrf(m)
    if info != 0:
        raise ValueError('dpotrf failed on input {}'.format(m))
    inv, info = lapack.dpotri(cholesky)
    if info != 0:
        raise ValueError('dpotri failed on input {}'.format(cholesky))
    upper_triangular_to_symmetric(inv)
    return inv

此版本比np.linalg.inv调用upper_triangular_to_symmetric.