Question Pourquoi la multiplication matricielle est-elle plus rapide avec numpy qu'avec ctypes en Python?


J'essayais de trouver le moyen le plus rapide de multiplier la matrice et j'ai essayé trois façons différentes:

  • Implémentation pure en python: pas de surprise ici.
  • Implémentation numpy utilisant numpy.dot(a, b)
  • Interfaçage avec C en utilisant ctypes module en Python.

C'est le code C qui est transformé en 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 parié que la version utilisant C aurait été plus rapide ... et j'aurais perdu! Ci-dessous est mon benchmark qui semble montrer que je l'ai soit mal fait, soit que numpy est bêtement rapide:

benchmark

Je voudrais comprendre pourquoi le numpy la version est plus rapide que la ctypes version, je ne parle même pas de la pure implémentation de Python car elle est assez évidente.


43
2018-05-04 03:36


origine


Réponses:


Je ne suis pas trop familier avec Numpy, mais la source est sur Github. Une partie des produits scalaires est implémentée dans https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/arraytypes.c.src, ce que je suppose est traduit en implémentations C spécifiques pour chaque type de données. Par exemple:

/**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**/

Cela semble calculer des produits scalaires à une dimension, c'est-à-dire sur des vecteurs. Dans mes quelques minutes de navigation sur Github, je n’ai pas pu trouver la source des matrices, mais il est possible qu’elle utilise un appel à FLOAT_dot pour chaque élément de la matrice de résultat. Cela signifie que la boucle dans cette fonction correspond à votre boucle la plus interne.

Une différence entre eux est que le "stride" - la différence entre les éléments successifs dans les entrées - est explicitement calculé une fois avant d'appeler la fonction. Dans votre cas, il n'y a pas de foulée et le décalage de chaque entrée est calculé à chaque fois, par ex. a[i * n + k]. Je m'attendrais à un bon compilateur pour l'optimiser à un niveau similaire à la foulée de Numpy, mais peut-être ne peut-il pas prouver que le pas est une constante (ou ce n'est pas optimisé).

Numpy peut également faire quelque chose d'intelligent avec des effets de cache dans le code de niveau supérieur qui appelle cette fonction. Une astuce courante consiste à réfléchir à la contiguïté de chaque ligne ou à chaque colonne - et à essayer d'aborder chaque partie contiguë. Il semble difficile d'être parfaitement optimal, pour chaque produit scalaire, une matrice d'entrée doit être parcourue par des lignes et l'autre par des colonnes (sauf si elles sont stockées dans un ordre majeur différent). Mais cela peut au moins le faire pour les éléments de résultat.

Numpy contient également du code pour choisir l'implémentation de certaines opérations, y compris "dot", à partir de différentes implémentations de base. Par exemple, il peut utiliser un BLAS bibliothèque. D'après la discussion ci-dessus, CBLAS est utilisé. Cela a été traduit de Fortran dans C. Je pense que l'implémentation utilisée dans votre test serait celle trouvée ici: http://www.netlib.org/clapack/cblas/sdot.c.

Notez que ce programme a été écrit par une machine pour une autre machine à lire. Mais vous pouvez voir au bas qu'il utilise une boucle déroulée pour traiter 5 éléments à la fois:

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

Ce facteur de déroulement a probablement été choisi après en avoir profilé plusieurs. Mais un avantage théorique en est que plus d'opérations arithmétiques sont effectuées entre chaque point de branchement et que le compilateur et le CPU ont plus de choix sur la manière de les planifier de manière optimale pour obtenir autant d'instructions en pipeline que possible.


19
2018-05-04 05:01



NumPy utilise une méthode BLAS optimisée et optimisée pour la multiplication de matrices (voir aussi: ATLAS). La fonction spécifique dans ce cas est GEMM (pour la multiplication matricielle générique). Vous pouvez rechercher l'original en recherchant dgemm.f (c'est dans Netlib).

L'optimisation, par ailleurs, va au-delà des optimisations du compilateur. Ci-dessus, Philip a mentionné Coppersmith-Winograd. Si je me souviens bien, cet algorithme est utilisé dans la plupart des cas de multiplication matricielle dans ATLAS (même si un commentateur note que cela pourrait être l’algorithme de Strassen).

En d'autres termes, votre matmult l'algorithme est l'implémentation triviale. Il existe des moyens plus rapides de faire la même chose.


26
2017-08-20 02:48



Le langage utilisé pour implémenter certaines fonctionnalités est une mauvaise mesure des performances. L'utilisation d'un algorithme plus approprié est souvent le facteur décisif.

Dans votre cas, vous utilisez l’approche naïve de la multiplication matricielle enseignée à l’école, qui est en O (n ^ 3). Cependant, vous pouvez faire beaucoup mieux pour certains types de matrices, par ex. matrices carrées, matrices de rechange et ainsi de suite.

Regardez le Algorithme Coppersmith-Winograd (multiplication par matrice carrée en O (n ^ 2.3737)) pour un bon point de départ pour la multiplication matricielle rapide. Voir aussi la section "Références", qui répertorie certains pointeurs vers des méthodes encore plus rapides.


Pour un exemple plus réaliste de gains de performances étonnants, essayez d’écrire rapidement strlen() et le comparer à l'implémentation glibc. Si vous n'arrivez pas à le battre, lisez les glibc strlen() source, il a assez bons commentaires.


10
2018-05-04 06:30



Numpy est également un code hautement optimisé. Il y a un essai sur certaines parties du livre Beau code.

Les ctypes doivent passer par une traduction dynamique de C en Python et inversement. Dans Numpy, la plupart des opérations matricielles sont entièrement internes.


2
2018-05-04 04:33



Les gars qui ont écrit NumPy savent évidemment ce qu'ils font.

Il existe de nombreuses façons d’optimiser la multiplication de la matrice. Par exemple, l'ordre dans lequel vous traversez la matrice affecte les modèles d'accès à la mémoire, qui affectent les performances.
Une bonne utilisation de SSE est un autre moyen d'optimiser, que NumPy emploie probablement.
Il y a peut-être d'autres moyens que les développeurs de NumPy connaissent et que je ne connais pas.

BTW, avez-vous compilé votre code C avec l'optiomisation?

Vous pouvez essayer l'optimisation suivante pour C. Il fonctionne en parallèle, et je suppose que NumPy fait quelque chose dans le même sens.
REMARQUE: Ne fonctionne que pour des tailles égales. Avec un travail supplémentaire, vous pouvez supprimer cette limitation et conserver l'amélioration des performances.

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

2
2018-05-04 04:26



La raison la plus souvent invoquée pour expliquer l’avantage de Fortran en termes de code numérique est que le langage facilite la détection. aliasing - le compilateur peut dire que les matrices en train de se multiplier ne partagent pas la même mémoire, ce qui peut aider à améliorer la mise en cache (inutile de s’assurer que les résultats sont réécrits immédiatement dans la mémoire "partagée"). C'est pourquoi C99 a introduit restreindre.

Cependant, dans ce cas, je me demande si le code numpy parvient également à en utiliser instructions spéciales que le code C n'est pas (car la différence semble particulièrement grande).


1
2018-05-04 04:31