爲什麼矩陣乘numpy的速度要快於Python中的ctypes？

Question

我試圖找出做矩陣乘法的最快方法，並試圖三種方式：

• 純Python實現：這裏沒有驚喜。
• 使用numpy.dot(a, b)
• Numpy實現在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版本比​​版本快，我我甚至沒有談論純粹的Python實現，因爲它很明顯。

Answer 1

我對Numpy不太熟悉，但是源碼在Github上。 dot產品的一部分在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的傢伙顯然知道他們在做什麼。

有很多方法可以優化矩陣乘法。例如，遍歷矩陣的順序影響內存訪問模式，這會影響性能。
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; 
     } 
    } 
} 

Answer 3

以數字代碼的Fortran的速度優勢給出的最常見的原因，據我所知，是該語言可以更容易地檢測aliasing - 編譯器可以告訴大家，矩陣相乘不共享相同的內存，這可以幫助改進緩存（無需確保將結果立即寫回到「共享」內存中）。這就是爲什麼C99引入restrict。

但是，在這種情況下，我想知道numpy代碼是否正在設法使用C代碼沒有的一些special instructions（因爲差別似乎特別大）。

Answer 4

Numpy也是高度優化的代碼。在Beautiful Code這本書中有一篇關於它的一部分的文章。

ctypes必須經歷從C到Python的動態轉換，並且返回會增加一些開銷。在Numpy中，大多數矩陣操作完全在內部完成。

Answer 5

用於實現某種功能的語言本身就是衡量性能的一個壞指標。通常，使用更合適的算法是決定性因素。

在你的情況，你正在使用天真的方法矩陣乘法在學校教，這是在O（n^3）。然而，對某些類型的矩陣你可以做得更好，例如方矩陣，備用矩陣等等。

查看Coppersmith–Winograd algorithm（O（n^2.3737）中的方陣乘法），以獲得快速矩陣乘法的良好起點。另請參閱「參考資料」一節，其中列出了一些更快的方法。

爲驚人的性能提升更樸實例如，嘗試寫一個快速strlen()並將其與glibc的實現。如果你不能擊敗它，閱讀glibc的strlen()來源，它有相當好的評論。

Answer 6

NumPy使用高度優化的，經過仔細調整的BLAS方法進行矩陣乘法（另請參閱：ATLAS）。這種情況下的具體功能是GEMM（用於通用矩陣乘法）。您可以通過搜索dgemm.f（它位於Netlib中）查找原件。

順便說一句，優化不僅僅是編譯器優化。上面，菲利普提到Coppersmith – Winograd。如果我沒有記錯的話，這是ATLAS中大多數矩陣乘法情況下使用的算法（儘管評論者指出它可能是Strassen的算法）。

換句話說，你的matmult算法是微不足道的實現。有更快的方法來做同樣的事情。