60

我试图找出进行矩阵乘法的最快方法并尝试了 3 种不同的方法：

  • 纯 python 实现：这里没有惊喜。
  • 使用 Numpy 实现numpy.dot(a, b)
  • ctypes在 Python 中使用模块与 C交互。

这是转换为共享库的 C 代码：

#include <stdio.h>
#include <stdlib.h>

void matmult(float* a, float* b, float* c, int n) {
    int i = 0;
    int j = 0;
    int k = 0;

    /*float* c = malloc(nay * sizeof(float));*/

    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            int sub = 0;
            for (k = 0; k < n; k++) {
                sub = sub + a[i * n + k] * b[k * n + j];
            }
            c[i * n + j] = sub;
        }
    }
    return ;
}

以及调用它的 Python 代码：

def C_mat_mult(a, b):
    libmatmult = ctypes.CDLL("./matmult.so")

    dima = len(a) * len(a)
    dimb = len(b) * len(b)

    array_a = ctypes.c_float * dima
    array_b = ctypes.c_float * dimb
    array_c = ctypes.c_float * dima

    suma = array_a()
    sumb = array_b()
    sumc = array_c()

    inda = 0
    for i in range(0, len(a)):
        for j in range(0, len(a[i])):
            suma[inda] = a[i][j]
            inda = inda + 1
        indb = 0
    for i in range(0, len(b)):
        for j in range(0, len(b[i])):
            sumb[indb] = b[i][j]
            indb = indb + 1

    libmatmult.matmult(ctypes.byref(suma), ctypes.byref(sumb), ctypes.byref(sumc), 2);

    res = numpy.zeros([len(a), len(a)])
    indc = 0
    for i in range(0, len(sumc)):
        res[indc][i % len(a)] = sumc[i]
        if i % len(a) == len(a) - 1:
            indc = indc + 1

    return res

我敢打赌，使用 C 的版本会更快……我会输的！下面是我的基准，它似乎表明我要么做错了，要么numpy太快了：

基准

我想了解为什么numpy版本比版本快ctypes，我什至不谈论纯 Python 实现，因为它很明显。

4

6 回答 6

37

NumPy 使用高度优化、仔细调整的 BLAS 方法进行矩阵乘法（另请参见：ATLAS）。这种情况下的具体功能是 GEMM（用于通用矩阵乘法）。dgemm.f您可以通过搜索（它在 Netlib 中）来查找原件。

顺便说一句，优化超越了编译器优化。上面，菲利普提到了 Coppersmith-Winograd。如果我没记错的话，这是用于 ATLAS 中大多数矩阵乘法的算法（尽管评论者指出它可能是 Strassen 的算法）。

换句话说，您的matmult算法是简单的实现。有更快的方法来做同样的事情。

于 2012-08-20T02:48:41.520 回答
28

我对 Numpy 不太熟悉，但源代码在 Github 上。部分点积在https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/arraytypes.c.src中实现，我假设它被翻译成每个特定的 C 实现数据类型。例如：

/**begin repeat
 *
 * #name = BYTE, UBYTE, SHORT, USHORT, INT, UINT,
 * LONG, ULONG, LONGLONG, ULONGLONG,
 * FLOAT, DOUBLE, LONGDOUBLE,
 * DATETIME, TIMEDELTA#
 * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 * npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 * npy_float, npy_double, npy_longdouble,
 * npy_datetime, npy_timedelta#
 * #out = npy_long, npy_ulong, npy_long, npy_ulong, npy_long, npy_ulong,
 * npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 * npy_float, npy_double, npy_longdouble,
 * npy_datetime, npy_timedelta#
 */
static void
@name@_dot(char *ip1, npy_intp is1, char *ip2, npy_intp is2, char *op, npy_intp n,
           void *NPY_UNUSED(ignore))
{
    @out@ tmp = (@out@)0;
    npy_intp i;

    for (i = 0; i < n; i++, ip1 += is1, ip2 += is2) {
        tmp += (@out@)(*((@type@ *)ip1)) *
               (@out@)(*((@type@ *)ip2));
    }
    *((@type@ *)op) = (@type@) tmp;
}
/**end repeat**/

这似乎计算一维点积，即在向量上。在我浏览 Github 的几分钟里，我找不到矩阵的来源，但它可能FLOAT_dot对结果矩阵中的每个元素都使用了一次调用。这意味着此函数中的循环对应于您最内层的循环。

它们之间的一个区别是“步幅”——输入中连续元素之间的差异——在调用函数之前显式计算一次。在您的情况下，没有步幅，并且每次都会计算每个输入的偏移量，例如a[i * n + k]. 我本来期望一个好的编译器将其优化为类似于 Numpy 步幅的东西，但也许它不能证明该步长是一个常数（或者它没有被优化）。

Numpy 也可能在调用此函数的高级代码中使用缓存效果做一些聪明的事情。一个常见的技巧是考虑每一行是连续的，还是每一列——并首先尝试迭代每个连续的部分。似乎很难做到完美，对于每个点积，一个输入矩阵必须按行遍历，另一个按列遍历（除非它们碰巧以不同的主要顺序存储）。但它至少可以为结果元素做到这一点。

Numpy 还包含从不同的基本实现中选择某些操作的实现的代码，包括“点”。例如，它可以使用BLAS库。从上面的讨论中，听起来像是使用了 CBLAS。这是从 Fortran 翻译成 C 的。我认为您的测试中使用的实现将是在这里找到的实现：http: //www.netlib.org/clapack/cblas/sdot.c。

请注意，该程序是由一台机器编写的，供另一台机器读取。但是您可以在底部看到它使用展开的循环一次处理 5 个元素：

for (i = mp1; i <= *n; i += 5) {
stemp = stemp + SX(i) * SY(i) + SX(i + 1) * SY(i + 1) + SX(i + 2) * 
    SY(i + 2) + SX(i + 3) * SY(i + 3) + SX(i + 4) * SY(i + 4);
}

这个展开因素很可能是在分析了几个之后才选择的。但是它的一个理论优势是每个分支点之间完成了更多的算术运算，并且编译器和 CPU 在如何优化调度它们以获得尽可能多的指令流水线方面有更多选择。

于 2012-05-04T05:01:05.363 回答
10

用于实现某种功能的语言本身就是一个不好的性能衡量标准。通常，使用更合适的算法是决定因素。

在您的情况下，您正在使用学校教授的矩阵乘法的简单方法，即 O(n^3)。但是，对于某些类型的矩阵，您可以做得更好，例如方阵、备用矩阵等。

查看Coppersmith–Winograd 算法（O(n^2.3737) 中的方阵乘法），作为快速矩阵乘法的良好起点。另请参阅“参考”部分，其中列出了一些指向更快方法的指针。


对于一个更实际的惊人性能提升示例，请尝试编写一个 faststrlen()并将其与 glibc 实现进行比较。如果您无法击败它，请阅读 glibc 的strlen()源代码，它有相当不错的评论。

于 2012-05-04T06:30:48.547 回答
5

编写 NumPy 的人显然知道他们在做什么。

有很多方法可以优化矩阵乘法。例如，您遍历矩阵的顺序会影响内存访问模式，从而影响性能。
充分利用 SSE 是另一种优化方法，NumPy 可能会采用这种方法。
可能还有更多方法，NumPy 的开发人员知道而我不知道。

顺便说一句，您是否使用优化编译了您的 C 代码？

您可以尝试对 C 进行以下优化。它确实可以并行工作，我想 NumPy 也做了同样的事情。
注意：仅适用于偶数尺寸。通过额外的工作，您可以消除此限制并保持性能改进。

for (i = 0; i < n; i++) {
        for (j = 0; j < n; j+=2) {
            int sub1 = 0, sub2 = 0;
            for (k = 0; k < n; k++) {
                sub1 = sub1 + a[i * n + k] * b[k * n + j];
                sub1 = sub1 + a[i * n + k] * b[k * n + j + 1];
            }
            c[i * n + j]     = sub;
            c[i * n + j + 1] = sub;
        }
    }
}
于 2012-05-04T04:26:52.903 回答
5

Numpy 也是高度优化的代码。Beautiful Code一书中有一篇关于它的部分内容的文章。

ctypes 必须经过从 C 到 Python 的动态转换并返回，这会增加一些开销。在 Numpy 中，大多数矩阵运算完全在其内部完成。

于 2012-05-04T04:33:13.567 回答
2

Fortran 在数字代码中的速度优势 afaik 给出的最常见原因是该语言更容易检测别名- 编译器可以判断被相乘的矩阵不共享相同的内存，这有助于改善缓存（不需要确保结果立即写回“共享”内存）。这就是 C99 引入限制的原因。

但是，在这种情况下，我想知道 numpy 代码是否也在设法使用 C 代码没有的一些特殊指令（因为差异似乎特别大）。

于 2012-05-04T04:31:36.193 回答