为什么矩阵乘法与numpy的速度比使用Python ctypes的？为什么矩阵乘法与numpy的速度

我试图找出做矩阵乘法的最快方法，并试图三种方式：

纯Python实现：这里没有惊喜。
numpy的实施使用numpy.dot(a, b)
用下，使用接口ctypes模块在Python。

这是将其转化入一个共享库中的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实现，因为它是一种显而易见的。

Answer 1:

我不是太熟悉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对如何最优地调度他们得到尽可能多的指令流水线尽可能更多的选择。

Answer 2:

NumPy的使用对于矩阵乘法高度优化，仔细调谐BLAS方法（也可参见： ATLAS ）。在这种情况下的特定功能是GEMM（用于通用矩阵乘法）。您可以通过搜索查找原始dgemm.f （它在入Netlib）。

优化，顺便说一句，超出编译器优化。以上，菲利普提到的铜匠，威诺格拉德。如果我没有记错，这是它在ATLAS用于矩阵乘法的大多数情况下（虽然评论者指出它可以被STRASSEN算法）的算法。

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

Answer 3:

用于实现某些功能的语言，是表现不好的措施本身。通常情况下，使用更合适的算法是决定因素。

在你的情况，你使用的简单方法是在学校，这是在为O（n ^ 3）教给矩阵乘法。然而，你可以做很多对某些种类的基质，例如方阵，备用矩阵等等更好。

看一看的铜匠，威诺格拉德算法对快速矩阵乘法一个很好的起点（方形矩阵乘法为O（n ^ 2.3737））。另见“参考”，其中列出了一些指针，以更快的速度的方法。

对于惊人的性能提升更朴实例如，尝试写一个快速strlen()并将其与glibc的实现。如果不设法打败它，阅读的glibc的strlen()源，它具有较好的意见。

Answer 4:

numpy的也高度优化代码。有关于这本书的部分内容的文章漂亮的代码。

该ctypes的要经过从C到Python和回动态翻译，增加了一些开销。在NumPy的大部分矩阵运算完成后完全内部给它。

Answer 5:

谁写的NumPy的家伙明明知道自己在做什么。

有许多方法来优化矩阵乘法。例如，在遍历矩阵顺序会影响内存访问模式，这会影响性能。
用好SSE的是另一种方式来优化，这可能NumPy的员工。
可能有更多的方式，这与NumPy的开发者都知道，我不知道。

顺便说一句，你编译optiomization 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;
        }
    }
}