在一段Python代码中,我需要在某个时刻分别将两个2x2矩阵的大列表相乘.在代码中,这两个列表都是形状为(n,2,2)的数字array.另一(n,2,2)数组中的预期结果,其中矩阵1是第一列表的矩阵1与第二列表的矩阵1之间的乘法的结果,依此类推.

经过一些分析后,我发现矩阵乘法是性能瓶颈.出于好奇,我试着"显式"地编写矩阵乘法.下面是一个带有测量运行时的代码示例.

import timeit
import numpy as np

def explicit_2x2_matrices_multiplication(
    mats_a: np.ndarray, mats_b: np.ndarray
) -> np.ndarray:
    matrices_multiplied = np.empty_like(mats_b)
    for i in range(2):
        for j in range(2):
            matrices_multiplied[:, i, j] = (
                mats_a[:, i, 0] * mats_b[:, 0, j] + mats_a[:, i, 1] * mats_b[:, 1, j]
            )

    return matrices_multiplied


matrices_a = np.random.random((1000, 2, 2))
matrices_b = np.random.random((1000, 2, 2))

assert np.allclose( # Checking that the explicit version is correct 
    matrices_a @ matrices_b,
    explicit_2x2_matrices_multiplication(matrices_a, matrices_b),
)

print(  # 1.1814142999992328 seconds
    timeit.timeit(lambda: matrices_a @ matrices_b, number=10000)
)
print(  # 1.1954495010013488 seconds
    timeit.timeit(lambda: np.matmul(matrices_a, matrices_b), number=10000)
)
print(  # 2.2304022700009227 seconds
    timeit.timeit(lambda: np.einsum('lij,ljk->lik', matrices_a, matrices_b), number=10000)
)
print(  # 0.19581600800120214 seconds
    timeit.timeit(
        lambda: explicit_2x2_matrices_multiplication(matrices_a, matrices_b),
        number=10000,
    )
)

如在代码中测试的,该函数产生与常规矩阵__matmul__结果相同的结果.然而,不同的是速度:在我的机器上,显式表达式的速度最多快10倍.

这对我来说是一个相当令人惊讶的结果.我原本预计NumPy表达式会更快,或者至少与更长的Python版本相当,而不是像我们在这里看到的那样慢一个数量级.我很想知道为什么业绩差异如此之大.

我运行的是NumPy版本1.25和Python版本3.10.6.

推荐答案

TL;DR:.问题中提供的所有方法都非常低效.事实上,Numpy显然是次优的,没有办法只使用Numpy高效地计算这一点.尽管如此,还是有比问题中提供的更快的解决方案.


解释和更快的实施

Numpy代码利用强大通用iterators将给定计算(如矩阵乘法)应用于多个数组切片.这样的迭代器对于实现broadcasting以及生成相对简单的einsum的实现是有用的.问题是,当迭代次数很大而数组很小时,它们也相当昂贵.这正是在您的用例中发生的事情.虽然可以通过优化Numpy代码来减少Numpy迭代器的开销,但在这个特定用例中,没有办法将开销减少到可以忽略的时间.事实上,每个矩阵只需要执行12次浮点运算.相对较新的x86-64处理器可以使用标量FMA单位在不到10纳秒的时间内计算每个矩阵.事实上,人们可以使用SIMD指令来计算每个矩阵,只需几纳秒.

首先,我们可以通过自己对第一个轴上的向量进行矩阵乘法来消除内部Numpy迭代器的开销.这就是explicit_2x2_matrices_multiplication所做的事情!

虽然explicit_2x2_matrices_multiplication应该要快得多,但它仍然不是最优的:它执行非连续的内存读取,创建几个无用的临时数组,并且每次Numy调用都会带来很小的开销.更快的解决方案是编写Numba/Cython代码,这样底层编译器就可以生成针对2x2矩阵进行优化的非常好的指令序列.在这种情况下,优秀的编译器甚至可以生成SIMD指令,这对于Numpy来说是不可能的.以下是生成的代码:

import numba as nb
@nb.njit('(float64[:,:,::1], float64[:,:,::1])')
def compute_fastest(matrices_a, matrices_b):
    assert matrices_a.shape == matrices_b.shape
    assert matrices_a.shape[1] == 2 and matrices_a.shape[2] == 2

    n = matrices_a.shape[0]
    result_matrices = np.empty((n, 2, 2))
    for k in range(n):
        for i in range(2):
            for j in range(2):
                result_matrices[k,i,j] = matrices_a[k,i,0] * matrices_b[k,0,j] + matrices_a[k,i,1] * matrices_b[k,1,j]

    return result_matrices

性能结果

以下是我的机器在配备i5-9600KF CPU的1000x2x2矩阵上的性能结果:

Naive einsum:                           214   µs
matrices_a @ matrices_b:                102   µs
explicit_2x2_matrices_multiplication:    24   µs
compute_fastest:                          2.7 µs   <-----

讨论

Numba实现达到了4.5G Flop.每个矩阵的计算时间仅为12个CPU周期(2.7 ns)!我的机器在实践中能够达到300G Flop(理论上是432GFlop),但只能达到50GFlop的单核和12.5GFlop的标量代码(理论上是18GFlop).操作的粒度太小,多个线程没有用处(生成线程的开销至少是几微秒).此外,SIMD码很难饱和De FMA单元,因为输入数据布局需要SIMD混洗,因此50G触发器实际上是一个乐观的上限.因此,我们可以有把握地说,the Numba implementation is a pretty efficient implementation.不过,多亏了SIMD instructions,才能写出更快的代码(我预计实际上会有大约x2的速度提升).也就是说,使用帮助编译器生成快速SIMD代码的SIMD内部函数编写本机代码真的很不容易(更不用说最终的代码将是难看的和难以维护的).因此,SIMD实现可能不值得花这么大力气.

Python相关问答推荐

DataFrame groupby函数从列返回数组而不是值

根据不同列的值在收件箱中移动数据

Gekko:Spring-Mass系统的参数识别

抓取rotowire MLB球员新闻并使用Python形成表格

如何根据参数推断对象的返回类型?

将图像拖到另一个图像

如何在虚拟Python环境中运行Python程序?

如何从在虚拟Python环境中运行的脚本中运行需要宿主Python环境的Shell脚本?

DataFrames与NaN的条件乘法

无法在Docker内部运行Python的Matlab SDK模块,但本地没有问题

为一个组的每个子组绘制,

如何使用Pandas DataFrame按日期和项目汇总计数作为列标题

lityter不让我输入左边的方括号,'

在Python中计算连续天数

Pandas在rame中在组内洗牌行,保持相对组的顺序不变,

为什么后跟inplace方法的`.rename(Columns={';b';:';b';},Copy=False)`没有更新原始数据帧?

如何在Django模板中显示串行化器错误

Scipy差分进化:如何传递矩阵作为参数进行优化?

用0填充没有覆盖范围的垃圾箱

如何在不不断遇到ChromeDriver版本错误的情况下使用Selify?