用 Cholesky 分解在 numpy 中反转矩阵的效率答案

【问题标题】：efficiency of inverting a matrix in numpy with Cholesky decomposition用 Cholesky 分解在 numpy 中反转矩阵的效率
【发布时间】：2017-11-04 20:00:42
【问题描述】：

我有一个对称的正定矩阵（例如协方差矩阵），我想计算它的逆矩阵。在数学上，我知道使用 Cholesky 分解来反转矩阵会更有效，尤其是在矩阵很大的情况下。但我不确定“numpy.lianlg.inv()”是如何工作的。假设我有以下代码：

import numpy as np

X = np.arange(10000).reshape(100,100)
X = X + X.T - np.diag(X.diagonal()) #  symmetry 
X = np.dot(X,X.T) # positive-definite

# simple inversion:
inverse1 = np.linalg.inv(X) 

# Cholesky decomposition inversion:
c = np.linalg.inv(np.linalg.cholesky(X))
inverse2 = np.dot(c.T,c)

哪个更有效（inverse1 或 inverse2）？如果第二个效率更高，为什么 numpy.linalg.inv() 不使用它呢？

【问题讨论】：

关于你的最后一个问题 - numpy 不知道你的矩阵是对称的，所以不能使用后一种方法。检查矩阵是否对称会很慢。
请注意，inv 也没有利用cholesky 是三角形的事实—它不使用 lapack 的 DTRTRI
你的代码没有为我运行，声称 X 不是正定的（可能是由于溢出）
一般情况下，倒置矩阵是bad idea。 inv 价格昂贵且数值不稳定。通常，您想将逆与向量相乘，即，您想求解方程组。在所有这些情况下，最好只使用 linalg.solve 之类的方法求解系统（告诉 solve 矩阵是对称的且正定的将使 solve 使用 Cholesky）。如果您想多次使用逆运算，请计算并存储分解以供以后使用。
如果您确实想反转 cholesky 因子，请使用 scipy.linalg.lapack.dtrtri

标签： python performance numpy matrix matrix-inverse

【解决方案1】：

使用以下设置：

import numpy as np

N = 100
X = np.linspace(0, 1, N*N).reshape(N, N)
X = 0.5*(X + X.T) + np.eye(N) * N

我通过IPython 的%timeit 得到以下时间：

In [28]: %timeit np.linalg.inv(X)
255 µs ± 30.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [29]: %timeit c = np.linalg.inv(np.linalg.cholesky(X)); np.dot(c.T,c)
414 µs ± 15.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

【讨论】：