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
2条答案
按热度按时间0s7z1bwu1#
根据您的应用程序,如果足以找到
Bs @ Bs.T ~ B
,则可以使用Cholesky分解。如果没有,你可以根据特征值分解得到平方根。Wall time: 4.4 s
-我们的基准Wall time: 52 ms
Cholesky更快Wall time: 1.62 s
快了两倍多(它探索了对称性)使用pytorch
Pytorch支持一些线性代数函数,它们可以跨多个CPU进行向量化
平方根使用eigh(12逻辑/ 6物理CPU)
Wall time: 400 ms
或乔列斯基分解
Wall time: 27 ms
y53ybaqx2#
为什么是新的答案?
另一个答案是Cholesky分解,它只适用于正定矩阵。虽然该限制对于样本协方差矩阵非常合适,但该问题指定了一个正半定矩阵,这将导致Scipy/Numpy中的
LinAlgError
。eig
和eigh
都适用于半正定,但正如另一个答案所示,这些都很慢。有一个与Cholesky相关的半定分解,LDLt分解。如果D是对角的,我们可以像特征分解一样简单地取元素的平方根。然而,在实践中,scipy
使用快速算法产生块对角D。解决方案
因为LDL算法产生的D最多有2x2个块,所以我们可以从
eigh
升级到eigh_tridiagonal
,这会快得多。此外,保证D具有与原始矩阵相同的(半)确定性,保证非负特征值。定时
下面的代码将
eigh
与ldl
+eigh_tridiagonal
的组合进行比较我们可以通过以下方式重建原始矩阵: