Pourquoi est-matrice de multiplication plus rapide avec numpy qu'avec ctypes en Python?

J'ai été à essayer de comprendre le moyen le plus rapide pour faire de multiplication de matrice et essayé de 3 façons différentes:

  • Pur python de mise en œuvre: pas de surprises ici.
  • Numpy mise en œuvre à l'aide de numpy.dot(a, b)
  • L'interfaçage avec C à l'aide de ctypes module en Python.

C'est le code C qui est transformée en une bibliothèque partagée:

#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 ;
}

Et le code Python qui l'appelle:

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

J'aurais pu parier que la version à l'aide de C aurait été plus rapide ... et j'aurais perdu ! Ci-dessous est ma référence en matière ce qui semble montrer que je soit fait de manière incorrecte, ou que numpy est bêtement rapide:

Pourquoi est-matrice de multiplication plus rapide avec numpy qu'avec ctypes en Python?

J'aimerais comprendre pourquoi le numpy version est plus rapide que la ctypes version, je ne parle même pas de la pure Python de mise en œuvre, car il est évident.

  • Belle question -, il s'avère np.(dot) est également plus rapide que le naïf mise en œuvre du GPU en C.
  • L'une des choses les plus importantes décisions de votre naïf C matmul lent, c'est l'accès à la mémoire de modèle. b[k * n + j]; à l'intérieur de la boucle interne (plus de k) a une foulée de n, de sorte qu'il touche une autre ligne de cache sur chaque accès. Et la boucle ne peut pas auto-vectorisation avec SSE/AVX. le Résoudre en transposant b à l'avant, d'un coût de O(n^2) et de la paie pour lui-même réduit les défauts de cache alors que vous ne O(n^3) charges de b. ce serait encore une implémentation naïve sans cache-blocage (aka boucle carrelage), si.
  • Puisque vous utilisez un int sum (pour une raison quelconque...), votre boucle, pourrait en fait vectoriser sans -ffast-math si la boucle interne a accès à deux séquentielle des tableaux. FP mathématiques n'est pas associatif, de sorte que les compilateurs ne pouvez pas re-ordre des opérations sans -ffast-math, mais de math entier est associative (et a une latence inférieure à celle FP outre, ce qui permet, si vous n'allez pas à optimiser votre boucle avec plusieurs accumulateurs ou d'autres latence cacher des trucs). float -> int les coûts de conversion la même chose qu'un FP add (en fait, à l'aide de la FP ajouter ALU sur les Processeurs Intel), il n'est donc pas la peine dans le code optimisé.