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

Question

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

我想知道是否存在针对对称半正定矩阵进行优化的算法，比numpy.linalg.inv() 更快（当然，如果它的实现可以从 python 轻松访问！）。我没有设法在numpy.linalg 中找到任何东西或在网上搜索。

编辑：

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

Answer 1

据我所知，对称矩阵没有标准的矩阵逆函数。通常，您需要对稀疏性等进行更多限制，以便为您的求解器获得良好的加速。但是，如果您查看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++ solvers 的类似问题，最终发现它确实取决于应用程序，您必须为您的问题选择最佳解决方案。

Answer 2

首先你需要确保协方差矩阵是正定的 (pd)而不是半-定的，否则矩阵是不可逆的。

没有 p.d.假设，矩阵求逆通常由 LU 分解完成，而对于 p.d.矩阵，可以使用 Cholesky 分解，这通常会降低计算成本。

在 Scipy 中，linalg.solve() 函数有一个参数 sym_pos，它假定矩阵是 p.d.。下面是一个简单的示例：

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 3

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

sp.linalg.lapack.dpotri的docstring如下

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) 更适合。

Answer 4

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

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%（请参阅 this question and answer）：

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。