根据您的应用程序，如果足以找到Bs @ Bs.T ~ B，则可以使用Cholesky分解。如果没有，你可以根据特征值分解得到平方根。

import numpy as np;
import scipy.linalg
A = np.random.randn(1500, 1500)

%%time
Bs = scipy.linalg.sqrtm(B)

Wall time: 4.4 s-我们的基准

%%time
Bs = scipy.linalg.cholesky(B)

Wall time: 52 ms Cholesky更快

D, V = scipy.linalg.eigh(B)
Bs = (V * np.sqrt(D)) @ V.T

Wall time: 1.62 s快了两倍多（它探索了对称性）

使用pytorch

Pytorch支持一些线性代数函数，它们可以跨多个CPU进行向量化

import torch.linalg
B_cpu = torch.tensor(B, device='cpu')

平方根使用eigh（12逻辑/ 6物理CPU）

%%time
D, V = torch.linalg.eigh(B_cpu)
Bs = (V * torch.sqrt(D)) @ V.T

Wall time: 400 ms
或乔列斯基分解

Bs = torch.linalg.cholesky(B_cpu)

Wall time: 27 ms

为什么是新的答案？

另一个答案是Cholesky分解，它只适用于正定矩阵。虽然该限制对于样本协方差矩阵非常合适，但该问题指定了一个正半定矩阵，这将导致Scipy/Numpy中的LinAlgError。eig和eigh都适用于半正定，但正如另一个答案所示，这些都很慢。有一个与Cholesky相关的半定分解，LDLt分解。如果D是对角的，我们可以像特征分解一样简单地取元素的平方根。然而，在实践中，scipy使用快速算法产生块对角D。

解决方案

因为LDL算法产生的D最多有2x2个块，所以我们可以从eigh升级到eigh_tridiagonal，这会快得多。此外，保证D具有与原始矩阵相同的（半）确定性，保证非负特征值。

定时

下面的代码将eigh与ldl + eigh_tridiagonal的组合进行比较

import numpy as np
from scipy.linalg import ldl, eigh, eigh_tridiagonal, cholesky
import time

ldl_time = 0
eig_time = 0
tri_time = 0
n_trials = 20
size = 400
for trial in range(n_trials):
    arr = np.random.normal(size=(size,size))
    arr = arr @ arr.T
    a = time.time()
    vals, vecs = eigh(arr)
    eig_time += time.time() - a
    a = time.time()
    l, d, p = ldl(arr)
    ldl_time += time.time() - a
    a = time.time()
    w, v = eigh_tridiagonal(np.diag(d), np.diag(d, 1))
    tri_time += time.time() - a
    a = time.time()
    L = cholesky(arr)
    chol_time += time.time() - a
print(f"LDL alone took {ldl_time:.3f} seconds")
print(f"    Additional factorization of D took {tri_time:.3f} seconds")
print(f"On the other hand, eigendecomposition took {eig_time:.3f} seconds")

LDL alone took 0.444 seconds
    Additional factorization of D took 0.041 seconds
On the other hand, eigendecomposition took 2.145 seconds

我们可以通过以下方式重建原始矩阵：

np.linalg.norm(l @ (v * w) @ v.T @ l.T - arr)

3.4043016350272527e-12