为什么我教授的 LU 分解版本比我的快？ Python 麻木答案

【问题标题】：Why is my prof's version of LU decomposition faster than mine? Python numpy为什么我教授的 LU 分解版本比我的快？ Python 麻木
【发布时间】：2018-07-20 19:12:17
【问题描述】：

我正在我的大学学习数值分析课程。我们正在研究LU decomposition。在查看讲师的版本之前，我尝试实施我的版本。我认为我的速度相当快，但实际上比较它们，我的讲师版本更快，即使它使用循环！这是为什么呢？

讲师版

def LU_decomposition(A):
    """Perform LU decomposition using the Doolittle factorisation."""

    L = np.zeros_like(A)
    U = np.zeros_like(A)
    N = np.size(A, 0)

    for k in range(N):
        L[k, k] = 1
        U[k, k] = (A[k, k] - np.dot(L[k, :k], U[:k, k])) / L[k, k]
        for j in range(k+1, N):
            U[k, j] = (A[k, j] - np.dot(L[k, :k], U[:k, j])) / L[k, k]
        for i in range(k+1, N):
            L[i, k] = (A[i, k] - np.dot(L[i, :k], U[:k, k])) / U[k, k]

    return L, U

我的版本

def lu(A, non_zero = 1):
    '''
    Given a matrix A, factorizes it into two matrices L and U, where L is
    lower triangular and U is upper triangular. This method implements
    Doolittle's method which sets l_ii = 1, i.e. L is a unit triangular
    matrix.

    :param      A: Matrix to be factorized. NxN
    :type       A: numpy.array

    :param non_zero: Value to which l_ii is assigned to. Must be non_zero.
    :type  non_zero: non-zero float.

    :return: (L, U)
    '''
    # Check if the matrix is square
    if A.shape[0] != A.shape[1]:
        return 'Input argument is not a square matrix.'

    # Store the size of the matrix
    n = A.shape[0]

    # Instantiate two zero matrices NxN (L, U)
    L = np.zeros((n,n), dtype = float)
    U = np.zeros((n,n), dtype = float)

    # Start algorithm
    for k in range(n):
        # Specify non-zero value for l_kk (Doolittle's)
        L[k, k] = non_zero
        # Case k = 0 is trivial
        if k == 0:
            # Complete first row of U
            U[0, :] = A[0, :] / L[0, 0]
            # Complete first column of L
            L[:, 0] = A[:, 0] / U[0, 0]
        # Case k = n-1 is trivial
        elif k == n-1:
            # Obtain  u_nn
            U[-1, -1] = (A[-1, -1] - np.dot(L[-1, :], U[:, -1])) / L[-1, -1]

        else:
            # Obtain u_kk
            U[k, k] = (A[k, k] - np.dot(L[k, :], U[:, k])) / L[k, k]
            # Complete kth row of U
            U[k, k+1:] = (A[k, k+1:] - [np.dot(L[k, :], U[:, i]) for i in \
                         range(k+1, n)]) / L[k, k]
            # Complete kth column of L
            L[k+1:, k] = (A[k+1:, k] - [np.dot(L[i, :], U[:, k]) for i in \
                         range(k+1, n)]) / U[k, k]
    return L, U

基准测试

我使用了以下命令：

A = np.random.randint(1, 10, size = (4,4))
%timeit lu(A)
57.5 µs ± 2.67 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit LU_decomposition(A)
42.1 µs ± 776 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

还有，scipy 的版本怎么这么好？

scipy.linalg.lu(A)
6.47 µs ± 219 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

【问题讨论】：

标签： python arrays numpy matrix scipy

【解决方案1】：

您的代码在 Python 代码中有条件，而讲座版本没有。 numpy 库在本机代码中进行了高度优化，因此您可以采取任何措施将计算推送到 numpy 而不是 python 中，这将有助于使其更快。

Scipy 必须在其库中对此进行更优化的版本，因为它是如何执行此操作的单个调用，外部循环可能是优化的本机代码的一部分，而不是相对较慢的 python 代码。

您可以尝试使用 Cython 进行基准测试，看看更优化的 python 运行时有什么不同。

【讨论】：

是的，scipy.linalg.lu 可能正在包装一些用 Fortran 或 C 编写的 BLAS 例程。
@juanpa.arrivillaga：更有可能是LAPACK。
@juanpa.arrivillaga 你认为有没有办法让 python 代码和 C 或 Fortran 一样快？
@Euler_Salter 你试过 Cython 了吗？这应该是一个好的开始。
说实话我不知道 Cython 是什么，但我会调查一下！

【解决方案2】：

我认为您的速度较慢，因为您使用了中间数据结构：

使用[np.dot(L[k, :], U[:, i]) for i in range(k+1, n)] 创建一个python 列表
使用A[k, k+1:] - temp_list 创建一个numpy 数组
另一个临时 numpy 数组是用temp_ndarray / L[k, k] 创建的
最后，这个临时数组被复制到结果数组中

对于这些步骤中的每一步，CPU 都必须执行一个循环，即使您没有明确编写一个循环。 Numpy 将这些循环抽象出来，但它们仍然必须执行！当然，通常在 numpy 中使用 X 个隐式快速循环而不是 1 个 python 循环可以获得回报，但这仅适用于中等大小的数组。此外，列表推导式实际上只比常规 for 循环快一点。

scipy 速度更快，因为它是一种在低级编程语言中高度优化的实现（而 python 是一种非常高级的语言）。最后，这可能意味着您应该欣赏教授的代码，因为它的优雅和可读性而不是它的速度:)

【讨论】：