np.linalg.solve 和 scipy.linalg.cho_solve 之间的性能差距

Question

我试图通过求解涉及 Cholesky 分解的两个线性方程组来找到 alpha。 scipy 有一个特殊的功能可以做到这一点。 scipy 和 numpy 之间存在显着的性能差距。我可以通过任何其他方式在 numpy 中实现与 scipy 一样好的性能吗？（假设不允许我使用scipy）。

import numpy as np
import scipy

def numpy_cho_solve(N,M):
    for seed in range(N):
        np.random.seed(seed)
        x = np.random.rand(M,1)
        y = np.random.rand(M,1)
        k = x@x.T + np.eye(M)# M*M
        L = np.linalg.cholesky(k)
        alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))
        
def scipy_cho_solve(N,M):
    for seed in range(N):
        np.random.seed(seed)
        x = np.random.rand(M,1)
        y = np.random.rand(M,1)
        k = x@x.T + np.eye(M)# M*M
        L = np.linalg.cholesky(k)
        alpha = scipy.linalg.cho_solve((L,True), y)

%timeit numpy_cho_solve(100,100)
%timeit scipy_cho_solve(100,100)

输出

317 ms ± 12.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
76.9 ms ± 3.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Answer 1

考虑到您只能使用 numpy，那么 np.linalg.solve 是求解线性方程的最佳函数，因为它给出了准确的结果。您可以使用numpy的inv和transpose函数，但准确的结果将是solve函数。

Answer 2

Cholesky 分解方法的正向和反向替换步骤非常快但不可向量化，因此 numpy 帮不上什么忙。你需要一个编译函数（如 scipy 实现） - 但如果你不能使用 scipy 我怀疑你可以使用 numba （通常用于制作 c 编译函数numpy).

np.linalg.solve 试图通过天真地应用 LU 替换来解决简单的前向替换步骤，因此它比专门构建的函数（甚至根本不使用 Cholesky）花费的时间长得多。

np.linalg.solve 和 scipy.linalg.cho_solve 之间的性能差距

Performance gap between np.linalg.solve and scipy.linalg.cho_solve

python

numpy

linear-algebra

scipy