我正在使用顺序蒙特卡罗方法来估计某种金融模型。然而，这种方法在计算上确实很昂贵，瓶颈归结为从多元高斯分布中采样。假设我使用的样本（粒子）数量为 5000，问题的维度为 5。这意味着我将有 5000 个形状为 (5,) 的不同均值向量，5000 个形状为 (5,5) 的不同协方差数组)，我需要从这 5000 对中抽取 5000 个样本。我该如何快速完成？到目前为止，我最快的解决方案是：
seeds = np.random.multivariate_normal(np.zeros(5), np.eye(5), size = 5000)
Ls = [np.linalg.cholesky(covs[i]) for i in range(5000)]
results = means + np.einsum('nij,njk->nik', Ls, seeds[:,:,np.newaxis])[:,:,0]
其中 np.einsum 是 numpy 中做爱因斯坦求和的函数，mean 和 covs 是上面提到的给定的 5000 对。
我认为在 for 循环中执行 cholesky 的步骤非常缓慢，但在这种情况下我找不到任何进行矢量化计算的方法。欢迎大家提出意见。