Dlaczego mnożenie macierzy jest szybsze przy użyciu numpy niż w przypadku ctypes w Pythonie?

Question

starałem się dowiedzieć najszybszy sposób to zrobić mnożenie macierzy i próbowałem 3 różne sposoby:

• Czysta wdrożenie Pythona: tutaj żadnych niespodzianek.
• realizacja Numpy pomocą numpy.dot(a, b)
• łączenia się z C, stosując ctypes moduł w Pythonie.

Jest to kod C, która przekształca się w udostępnionej biblioteki:

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

a kod Pythona, który wywołuje go:

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 

musiałbym założyć, że wersja użyciu C byłoby szybciej ... i bym się zgubił! Poniżej jest mój punkt odniesienia, który wydaje się wskazywać, że albo zrobił to nieprawidłowo, albo że numpy jest głupio szybki:

Chciałbym zrozumieć, dlaczego wersja numpy jest szybszy niż w wersji ctypes I” Nawet nie mówię o czystej implementacji Pythona, ponieważ jest to oczywiste.

Answer 1

Nie jestem zbyt zaznajomiony z Numpy, ale źródło jest na Github. Część produktów dot jest zaimplementowana w https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/arraytypes.c.src, która, jak przypuszczam, jest tłumaczona na konkretne implementacje C dla każdego typu danych. Na przykład:

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

Wydaje się obliczać jednowymiarowe produkty punktowe, tj. Na wektorach. W ciągu kilku minut przeglądania Github nie mogłem znaleźć źródła dla macierzy, ale możliwe jest, że używa jednego połączenia do FLOAT_dot dla każdego elementu w macierzy wyników. Oznacza to, że pętla w tej funkcji odpowiada twojej wewnętrznej pętli.

Jedną z różnic między nimi jest to, że "krok" - różnica między kolejnymi elementami na wejściach - jest jednoznacznie obliczany jeden raz przed wywołaniem funkcji. W twoim przypadku nie ma kroku, a przesunięcie każdego wejścia jest obliczane za każdym razem, np. a[i * n + k]. Spodziewałbym się, że dobry kompilator zoptymalizuje to do czegoś podobnego do kroku Numpy, ale może nie może udowodnić, że krok jest stały (lub nie jest zoptymalizowany).

Numpy może również robić coś inteligentnego dzięki efektom pamięci podręcznej w kodzie wyższego poziomu, który wywołuje tę funkcję. Typową sztuczką jest myślenie o tym, czy każdy wiersz jest ciągły, czy też o każdą kolumnę - i spróbuj najpierw powtórzyć każdą z sąsiadujących części. Wydaje się trudne, aby być całkowicie optymalnym, dla każdego produktu punktowego, jedna matryca wejściowa musi być przesuwana przez wiersze, a druga przez kolumny (chyba, że ​​akurat były przechowywane w różnej kolejności). Ale może to zrobić przynajmniej dla elementów wynikowych.

Numpy zawiera również kod do wyboru implementacji niektórych operacji, w tym "kropka", z różnych podstawowych implementacji. Na przykład może korzystać z biblioteki BLAS. Z powyższej dyskusji brzmi jak CBLAS. Zostało to przetłumaczone z Fortran na C. Myślę, że implementacja użyta w twoim teście byłaby znaleziona tutaj: http://www.netlib.org/clapack/cblas/sdot.c.

Należy zauważyć, że ten program został napisany przez komputer do odczytu z innego komputera. Ale można zobaczyć na dole, że to za pomocą odwiniętą pętlę na przetwarzanie 5 elementów naraz:

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

Ten rozwijanie czynnikiem jest prawdopodobnie były zbierane po profilowania kilka. Ale jedną z jego zalet jest to, że więcej operacji arytmetycznych jest wykonywanych między punktami rozgałęzień, a kompilator i procesor mają większy wybór, jak optymalnie je zaplanować, aby uzyskać jak najwięcej potoków instrukcji.

Answer 2

Faceci, którzy napisali NumPy, oczywiście wiedzą, co robią.

Istnieje wiele sposobów optymalizacji mnożenia macierzy. Na przykład kolejność przemierzania macierzy wpływa na wzorce dostępu do pamięci, które mają wpływ na wydajność.
Dobre wykorzystanie SSE to kolejny sposób na optymalizację, którą prawdopodobnie używa NumPy.
Może być więcej sposobów, o których wiedzą twórcy NumPy, a ja nie.

BTW, czy skompilowałeś swój kod C z optyfikacją?

Możesz wypróbować następującą optymalizację dla C. Działa ona równolegle, i przypuszczam, że NumPy robi coś wzdłuż tych samych linii.
UWAGA: Działa tylko w przypadku rozmiarów. Dzięki dodatkowej pracy możesz usunąć to ograniczenie i utrzymać wydajność.

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

Najczęstszą przyczyną podane dla prędkości korzyścią Fortran w kodzie numerycznym, AFAIK, jest to, że język sprawia, że ​​łatwiej wykryć aliasing - kompilator może powiedzieć, że matryce są mnożone nie podziela tą pamięć, co może pomóc w poprawie buforowania (nie ma potrzeby, aby pewne wyniki były natychmiast zapisywane w pamięci "współdzielonej"). Właśnie dlatego C99 wprowadziło restrict.

Jednak w tym przypadku, zastanawiam się, czy również kod numpy jest w stanie użyć jakiegoś special instructions, że kod C nie jest (ponieważ różnica wydaje się szczególnie duża).

Answer 4

Numpy to również wysoce zoptymalizowany kod. Istnieje esej na temat jego części w książce Beautiful Code.

The ctypes musi przejść przez dynamiczne tłumaczenie z C na Python iz powrotem, co dodaje trochę narzutów. W Numpy większość operacji macierzowych jest wykonywanych całkowicie wewnętrznie.

Answer 5

Język używany do implementacji określonej funkcjonalności sam w sobie jest złą miarą wydajności. Często decydującym czynnikiem jest zastosowanie bardziej odpowiedniego algorytmu.

W twoim przypadku używasz naiwnego podejścia do mnożenia macierzy, jak uczy się w szkole, która jest w O (n^3). Można jednak zrobić znacznie lepiej dla niektórych rodzajów matryc, np. macierze kwadratowe, macierze zapasowe i tak dalej.

Spójrz na Coppersmith–Winograd algorithm (mnożenie macierzy kwadratowych w O (n^2,3737)) dla dobrego punktu wyjścia do szybkiego mnożenia macierzy. Zobacz także sekcję "Referencje", która podaje niektóre wskaźniki do jeszcze szybszych metod.

Dla bardziej ziemskiego przykładu zadziwiającego przyrostu wydajności, spróbuj napisać szybki strlen() i porównaj go z implementacją glibc. Jeśli nie uda ci się go pokonać, przeczytaj źródło glibc strlen(), ma dość dobre komentarze.

Answer 6

NumPy wykorzystuje wysoce zoptymalizowaną, dokładnie dostrojoną metodę BLAS do multiplikacji macierzy (patrz także: ATLAS). Szczególną funkcją w tym przypadku jest GEMM (dla ogólnego mnożenia macierzy). Możesz wyszukać oryginał, wyszukując numer dgemm.f (w Netlib).

Optymalizacja, przy okazji, wykracza poza optymalizacje kompilatora. Powyżej Philip wspomniał o Coppersmithu – Winograd. Jeśli dobrze pamiętam, jest to algorytm używany w większości przypadków mnożenia macierzy w ATLAS (chociaż komentator zauważa, że ​​może to być algorytm Strassena).

Innymi słowy, Twój algorytm matmult jest banalną implementacją. Istnieją szybsze sposoby na zrobienie tego samego.