logo标签列表

已知矩阵对称且半正定，求矩阵的更高效的求逆方法

pythonnumpyscipylinear-algebra
7

我在使用Python中的numpy翻转协方差矩阵。协方差矩阵是对称且半正定的。

我想知道是否存在一种针对对称半正定矩阵进行优化的算法，比numpy.linalg.inv()更快（当然，如果可以从Python轻松访问其实现！）。我在numpy.linalg或搜索网络时没有找到相关信息。

编辑：

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

“?SY”系列适用于对称矩阵，包括复值矩阵。对于Hermitian矩阵，您需要使用“?HE”系列。关于不存在的问题，只有API不可用。大多数LAPACK函数都已经存在，未来还将增加更多函数。 - percusse
4个回答
8

我试用了@percusse的答案，但是在计时执行时，我发现它比使用100K个随机正定的4x4 np.float64矩阵的np.linalg.inv慢约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

我尝试对其进行分析，惊讶地发现它在调用upper_triangular_to_symmetric时占用了约82%的时间（但这并不是“难”部分）。我认为这是因为它在执行浮点数加法以合并矩阵时，而不是简单的复制。
我尝试了一个实现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快了大约68%，并且仅花费了约42%的时间调用upper_triangular_to_symmetric。

6
首先，您需要确保协方差矩阵是正定的，而不是半定的，否则该矩阵不可逆。如果没有正定假设，矩阵求逆通常通过LU分解来完成，而对于正定矩阵，则可以使用Cholesky分解，这通常可以降低计算成本。在Scipy中，linalg.solve()函数有一个sym_pos参数，可以假设矩阵是正定的。以下是一个快速示例：
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()更慢。
5
据我所知，对于对称矩阵，没有标准的矩阵求逆函数。一般来说，您需要更多的稀疏性等约束条件才能为求解器获得更好的加速。但是，如果您查看 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）需要厄米矩阵，因此可能更快。但是，只有在矩阵稀疏、仅求解特定解向量b并且可能根据矩阵结构实际上不会更快时，才是明智的选择。您真的必须进行基准测试才能找出答案。

我问过类似C++求解器的问题，最终发现这真的取决于应用程序，并且您必须为您的问题选择最佳求解器。

5

目前还没有这个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)。

关于最后的观察，我正在计算加权平均值的协方差矩阵：V = (V_1^{-1} + V_2^{-1} + ...) ^ {-1}。一般来说，V_i不是同时对角化的，所以据我所知，我必须执行所有的倒置操作。 - Giacomo Petrillo
1
我尝试了这个建议，发现它实际上比 np.linalg.inv 更慢，因为你将上三角矩阵转换为对称矩阵的方式（浮点数运算很慢 :)）。请参见我的答案以获得改进。 - Kerrick Staley
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接