numpy - 加速矩阵求逆

Question

我需要按顺序运行此代码数千次:

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成为 prior_mu , V_posterior成为 prior_V , a_posterior成为 prior_a , 和 b_posterior成为 prior_b为下一次通话。 y 和 x 在每次调用中都不同。

这个函数非常慢——运行大约需要 8 秒。这是因为规模。我有 ~5000 个参数，所以 prior_mu是 (1, 5000), prior_V是 (5000,5000) 并且是对称的 positive-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

Answer 1

就稳定性而言，写 solve(A, unit_matrix) 几乎总是更好。而不是 inv(A) .不过，这对性能没有帮助。

这里线性代数的性能几乎肯定是由底层 LAPACK 库修复的。库存 ATLAS 可能最慢，OpenBLAS 或 MKL 更好，有时更好。

但是，我很确定这里的主要改进确实是算法。首先，对于 PSD 矩阵，Cholecky (cholesky/cho_solve) 应该更好。其次，您似乎正在进行一级更新( x.T @x )，通常可以在 N**2 中实现通过 Shermann-Morrison 公式的一些变体进行运算，而不是 N**3用于直接反转。