python - 为什么 scipy.linalg 的方法在 9x9 矩阵上运行速度较慢？

Question

我正在比较解决线性系统的各种方法的运行时间，并且发现了一个奇怪的模式。我正在测试的解决方法是 la.solve() , la.inv() ，和la.lu_factor_solve() .

import scipy.linalg as la
import numpy as np
from time import time
from matplotlib import pyplot as plt

N = 20  # up to NxN matrix
T = 100  # run T times
inv_time, solve_time = [[] for _ in range(20)], [[] for _ in range(20)]
lu_factor_solve, lu_just_solve = [[] for _ in range(20)], [[] for _ in range(20)]
for _ in range(T):
    for n in range(1, N + 1):
        A = np.random.rand(n, n)
        b = np.random.rand(n, 1)
        np.dot(la.inv(A), b)  # the first time through is always slow,
        la.solve(A, b)        # so we run it once to get it out of the way
        start = time()
        np.dot(la.inv(A), b)
        end = time()
        inv_time[n - 1].append(end - start)

        start = time()
        la.solve(A, b)
        end = time()
        solve_time[n - 1].append(end - start)

        start = time()
        la.lu_solve(la.lu_factor(A), b)
        end = time()
        lu_factor_solve[n - 1].append(end - start)

        temp = la.lu_factor(A)
        start = time()
        la.lu_solve(temp, b)
        end = time()
        lu_just_solve[n - 1].append(end - start)
inv_time = np.mean(np.array(inv_time), axis=1)
solve_time = np.mean(np.array(solve_time), axis=1)
lu_factor_solve = np.mean(np.array(lu_factor_solve), axis=1)
lu_just_solve = np.mean(np.array(lu_just_solve), axis=1)
# do some plots
plt.plot(range(1, N + 1), inv_time, '-o', label='by inverse')
plt.plot(range(1, N + 1), solve_time, '-o', label='by la.solve()')
plt.plot(range(1, N + 1), lu_factor_solve, '-o', label='by lu factor solve')
plt.yscale('log')
plt.plot(range(1, N + 1), lu_just_solve, '-o', label='just la.lu_solve()')
plt.legend()
plt.show()

这些在随机矩阵上运行 A和列向量b由 A = np.random.rand(n, n) 制作和b = np.random.rand(n, 1)对于 n 的值从 1 到 20。我发现每次运行该程序时，它找到 9x9 矩阵的解的速度比其他大小的矩阵慢得多，如下图所示。红线显示仅执行 la.lu_solve() 所需的时间。以下是 T = 1 的结果图:

这是 T = 100 的结果:

这是否与 9x9 矩阵的固有特性有关，某些优化不适用于此大小的矩阵，或者与其他东西有关？

Answer 1

我已经用 perfplot 尝试了您的代码示例(我的一个小项目，本质上是 timeit 的包装)并且没有发现这样的特殊性。不过，可以认识到一级缓存已耗尽:

这是使用 NumPy 1.13.3 和 SciPy 1.0.0rc1。

仅运行一次测试时，代码示例中出现峰值的原因是它们运行得太快，以至于对机器其他部分的微小干扰可能会变得很大。这可以是任何东西，从浏览器中 JS 引擎的重复任务到从深度 sleep 状态唤醒的 CPU。事实上，您在 n=9 处看到这些峰值纯属巧合。当仅运行一次测试时，我自己可以重现 5 到 20 之间的随机峰值。

<小时/>

重现情节的代码:

import numpy
import perfplot
import scipy.linalg as la


def solve(data):
    A, b = data
    return la.solve(A, b)


def dot_inv(data):
    A, b = data
    return numpy.dot(la.inv(A), b)


def lu_solve(data):
    A, b = data
    return la.lu_solve(la.lu_factor(A), b)


perfplot.show(
    setup=lambda n: (numpy.random.rand(n, n), numpy.random.rand(n)),
    kernels=[solve, dot_inv, lu_solve],
    n_range=range(1, 40),
    )