为什么matrix与numpy相比，乘法速度要快于Python中的ctypes？
我试图找出matrix乘法的最快方法，并尝试了3种不同的方法：
- 纯Python的实现：这里没有什么惊喜。
-
numpy.dot(a, b)实现使用numpy.dot(a, b) - 在Python中使用
ctypes模块与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的实现，因为它很明显。
我对Numpy不是太熟悉，但是源码在Github上。 部分dot产品是在https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/arraytypes.c.src中实现的，我假设它被转换为每个特定的C实现数据types。; 例如：
/**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浏览器的几分钟内，我无法findmatrix的来源，但可能FLOAT_dot结果matrix中的每个元素使用一次FLOAT_dot调用。 这意味着这个函数中的循环对应于你最内层的循环。
它们之间的一个区别在于，“步幅” – input中连续元素之间的差异 – 在调用函数之前明确计算一次。 在你的情况下，没有步幅，每次计算每个input的偏移，例如a[i * n + k] 。 我会期望一个好的编译器将其优化到类似于Numpy的步伐，但也许不能certificate这个步骤是一个常量（或者没有被优化）。
在调用这个函数的高级代码中，Numpy也可能在caching效果上做了一些聪明的事情。 一个常见的技巧是考虑每行是连续的，还是每列 – 并尝试先遍历每个连续的部分。 对于每个点积来说，似乎很难做到完美最佳，一个inputmatrix必须被行和另一个列所遍历（除非它们碰巧以不同的主要顺序存储）。 但是至less可以为结果元素做到这一点。
Numpy还包含代码，用于从不同的基本实现中select特定操作（包括“点”）的实现。 例如它可以使用BLAS库。 从上面的讨论，听起来像使用CBLAS。 这是从Fortran翻译成C的。我认为在你的testing中使用的实现将是在这里find的： 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有更多的select来优化如何尽可能多地获得指令stream水线。
NumPy使用高度优化的，经过仔细调整的BLAS方法进行matrix乘法（另见： ATLAS ）。 这种情况下的具体function是GEMM（用于通用matrix乘法）。 您可以通过searchdgemm.f （它在Netlib中）来查找原始文件。
顺便说一句，优化超越了编译器优化。 在上面，菲利普提到了Coppersmith-Winograd。 如果我没有记错的话，这是在ATLAS中大多数matrix乘法情况下使用的algorithm（尽pipe一个评论者指出它可能是Strassen的algorithm）。
换句话说，你的matmultalgorithm是微不足道的实现。 有更快的方法来做同样的事情。
用来实现某种function的语言本身就是衡量性能的一个不好的方法。 通常，使用更合适的algorithm是决定因素。
在你的情况，你正在使用天真的方法matrix乘法在学校教，这是O（n ^ 3）。 但是，对于某些种类的matrix，例如平方matrix，备用matrix等，可以做得更好。
查看Coppersmith-Winogradalgorithm （O（n ^ 2.3737）中的方阵乘法），以获得快速matrix乘法的良好起点。 另请参阅“参考资料”一节，其中列出了一些更快的方法。
为了获得令人惊叹的性能增益的更加朴实的例子，试着写一个快速strlen()并将其与glibc实现进行比较。 如果你不能打败它，请阅读glibc的strlen()源代码，它有相当好的评论。
Numpy也是高度优化的代码。 “ 美码 ”一书中有部分文章。
ctypes必须经历从C到Python的dynamic转换，并且返回，这增加了一些开销。 在Numpy中，大多数matrix操作完全在内部完成。
写NumPy的人显然知道他们在做什么。
有很多方法可以优化matrix乘法。 例如，遍历matrix会影响内存访问模式，这会影响性能。
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; } } }
Fortran在数字代码afaik中速度优势的最常见原因是，该语言使得检测别名更容易 – 编译器可以告诉被乘的matrix不共享相同的内存，这可以帮助改进caching（不需要确保将结果立即写回到“共享”存储器中）。 这就是为什么C99引入了限制 。
然而，在这种情况下，我不知道是否numpy代码正在设法使用一些特殊的指令 ，C代码不是（因为差异似乎特别大）。
