当维度太高时，numpy multivariate_normal 错误

Question

我正在做一个家庭作业，我注意到当均值和协方差的维度非常高时，multivariate_normal 将永远占用所有 CPU，而不会产生任何结果。

示例代码sn-p，

cov_true  = np.eye(p)
mean_true = np.zeros(p)
beta_true = multivariate_normal(mean_true, cov_true, size=1).T

当p=5000 时，它将永远运行。 环境，python3.4和python3.5，numpy 1.11.0

这真的是一个错误还是我错过了什么？

Answer 1

什么需要这么多时间？

要考虑协方差矩阵的变量 NumPy computes the singular value decomposition 之间的关系，这需要大部分时间（underlying GESDD 通常是 Θ(n³) 和 5000³ 已经有点）了。

如何加快速度？

在所有变量独立的最简单情况下，您可以使用random.normal：

from numpy.random import normal

sample = normal(means, deviations, len(means))

否则，如果您的协方差矩阵恰好是满秩（因此是正定的），则通常用 cholesky 代替 svd（仍然是 Θ(n³)，但使用较小的常数）：

from numpy.random import standard_normal
from scipy.linalg import cholesky

l = cholesky(covariances, check_finite=False, overwrite_a=True)
sample = means + l.dot(standard_normal(len(means)))

如果矩阵可能是单数的（有时是这种情况），那么要么包装SPSTRF，要么考虑使用scipy#6202。

Cholesky 可能会明显更快，但如果这还不够，那么您可以进一步研究是否无法解析地分解矩阵，或者尝试使用不同的基础库（例如 ACML、MKL 或cuSOLVER）。