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

Question

我试图通过求解两个涉及乔列斯基分解的线性方程组来寻找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

Cholesky分解方法的前向和后向替换步骤非常快，但不能向量化，因此numpy不能提供太多帮助。你需要一个编译过的函数(作为scipy的实现)--但是如果你不能使用scipy，我怀疑你能不能使用numba (它通常用于为numpy生成c编译的函数)。

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

Answer 2

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