numpy - 加速矩阵求逆
问题描述
我需要按顺序运行此代码数千次:
def update_posterior(y, x, prior_mu, prior_V, prior_a, prior_b):
# Building blocks used to keep following calculation cleaner
prior_cov_inverse = np.linalg.inv(prior_V)
x_transpose = x.transpose()
n = len(y)
residuals = y - np.dot(x, prior_mu.transpose())
# Calculation of posterior parameters
V_posterior = np.linalg.inv((prior_cov_inverse + x_transpose * x))
mu_posterior = (V_posterior * (prior_cov_inverse * prior_mu.transpose() + x_transpose * y)).transpose()
a_posterior = prior_a + n/2
b_posterior = np.asscalar(prior_b + (residuals.transpose() * np.linalg.inv((np.identity(n) + x * prior_V * x_transpose)) * residuals)/2)
return mu_posterior, V_posterior, a_posterior, b_posterior
它的工作方式是将函数的输出反馈给它,并在下一次调用中使用mu_posterior
become prior_mu
、V_posterior
become prior_V
、a_posterior
becomeprior_a
和b_posterior
become prior_b
。y 和 x 在每次调用中都不同。
这个函数非常慢——运行大约需要 8 秒。这是因为规模。我有约 5000 个参数,所以prior_mu
(1, 5000)、prior_V
is (5000,5000) 和 symmetricpositive-definite
和prior_a
, 和prior_b
是标量。y 是标量,x 是 (1, 5000)。
以下是按行划分的时间:
3.75s: prior_cov_inverse = np.linalg.inv(prior_V)
3.86s: V_posterior = np.linalg.inv((prior_cov_inverse + x_transpose * x))
0.13s: b_posterior = np.asscalar(prior_b + (residuals.transpose() * np.linalg.inv((np.identity(n) + x * prior_V * x_transpose)) * residuals)/2)
有什么想法可以加快速度吗?我尝试使用 Cholesky 分解,但它更慢?!我的猜测是有一种更有效的方法可以在 Python 中实现 Cholesky 分解。
prior_cov_inverse2 = np.linalg.inv(np.linalg.cholesky(prior_V))
prior_cov_inverse2 = np.dot(prior_cov_inverse2.transpose(), prior_cov_inverse2)
编辑:这里有一些示例数据来说明这个问题......
import numpy as np
prior_mu = np.asmatrix(np.full((1, 5040), 5))
prior_V = np.diagflat(np.asmatrix(np.full((1, 5040), 30))) #usually not diagonal, but always symmetric positive definitive
a = 2
b = 2
y = np.asmatrix([10])
x = np.asmatrix(np.concatenate(([1], np.zeros(5039))))
print(update_posterior(y, x, prior_mu, prior_V, a, b))
编辑二:
通过删除矩阵求逆以支持求解以及使用 Sherman Morrison 公式,我已经能够将其从 ~8s/run 降低到 ~1.4s/run。这是我当前的代码。如果有人对如何进一步加快速度有任何想法,请分享!:)
def update_posterior(y, x, prior_mu, prior_V, prior_a, prior_b, I):
# Building blocks used to keep following calculation cleaner
x_transpose = x.transpose()
n = len(y)
residuals = y - np.dot(x, prior_mu.transpose())
# Calculation of posterior parameters
# Below is equivalent to np.linalg.inv(prior_V_inverse + np.dot(x_transpose, x)) but significantly faster
V_posterior = prior_V - np.true_divide(np.linalg.multi_dot((prior_V, x_transpose, x, prior_V)), 1 + np.matmul(np.matmul(x, prior_V), x_transpose))
# Below is equivalent to mu_posterior = np.dot(V_posterior, (np.matmul(prior_V_inverse, prior_mu.transpose()) + np.matmul(x_transpose, y))).transpose() but significantly faster
mu_posterior = np.dot(V_posterior, np.linalg.solve(prior_V, prior_mu.transpose()) + np.matmul(x_transpose, y)).transpose()
a_posterior = prior_a + n/2
b_posterior = np.asscalar(prior_b + (np.matmul(np.matmul(residuals.transpose(), np.linalg.inv((np.identity(n) + np.matmul(np.matmul(x, prior_V), x_transpose)))), residuals))/2)
return mu_posterior, V_posterior, a_posterior, b_posterior
解决方案
就稳定性而言,写solve(A, unit_matrix)
而不是inv(A)
. 不过,这对性能没有帮助。
这里线性代数的性能几乎肯定是由底层 LAPACK 库修复的。库存 ATLAS 可能最慢,OpenBLAS 或 MKL 更好,有时更好。
但是,我很确定这里的主要改进确实是算法。首先,对于 PSD 矩阵,Cholecky (cholesky/cho_solve) 应该更好。其次x.T @x
,您似乎正在执行一级更新(N**2
N**3
推荐阅读
- node.js - 无法在 ejs 模板中显示图像。我正在使用快速文件上传器、mongodb、express、node、ejs
- java - org.hibernate.boot.registry.classloading.spi.ClassLoadingException:无法加载类 [CredDataTmplEntry] ...原因:java.lang.VerifyError
- flutter - 如何在 Flutter 中的所有应用程序中分配 BottomNavigationBar?
- jenkins - 如何从自动触发詹金斯构建的 Gitlab 更改中排除某些文件
- mysql - 从 wordpress 数据库结构中获得最高投票结果
- elasticsearch - Elasticsearch 7 SQL CLI : ./x-pack-env: No such file or directory 错误
- php - 使用 PHP Laravel 2019 在 postgres 中加入来自不同数据库的 2 个或多个表
- javascript - 如何在 Anguar 中使用通过 state 传递的数据
- swift - 无需 UIDocumentPickerViewController 直接从外部驱动器访问文件
- c - C语言中的12字节整数