Почему умножение матриц быстрее в numpy, чем в ctypes в Python?

Я пытался найти самый быстрый способ умножения матриц и пробовал 3 разных способа:

Pure python implementation: no surprises here. Numpy implementation using numpy.dot(a, b) Interfacing with C using ctypes module in Python.

Это код C, который преобразуется в общую библиотеку:

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

И код Python, который вызывает его:

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

Я бы поспорил, что версия, использующая C, была бы быстрее ... и я проиграл! Ниже мой тест, который, кажется, показывает, что я сделал это неправильно, или чтоnumpy тупо быстро

benchmark

Я хотел бы понять, почемуnumpy версия быстрее, чемctypes версия, я даже не говорю о чистой реализации Python, так как это отчасти очевидно.

 Peter Cordes28 авг. 2017 г., 08:19
Одна из самых больших вещей, делающих вашу наивную матрешку медленной, - это модель доступа к памяти.b[k * n + j]; внутри внутренней петлиk) имеет успехn, так что он касается другой строки кэша при каждом доступе. И ваш цикл не может автоматически векторизоваться с помощью SSE / AVX.Solve this by transposing b up-front, which costs O(n^2) time and pays for itself in reduced cache misses while you do O(n^3) loads from b.  Это все равно будет наивной реализацией без блокировки кеша (или циклического разбиения).
 user239802914 нояб. 2012 г., 18:10
Хороший вопрос - оказывается, np.dot () также быстрее, чем простая реализация на GPU в C.
 Peter Cordes28 авг. 2017 г., 08:21
Так как вы используетеint sum (по какой-то причине ...), ваш цикл может фактически векторизовать без-ffast-math если внутренний цикл обращался к двум последовательным массивам. Математика FP не ассоциативна, поэтому компиляторы не могут переупорядочивать операции без-ffast-math, но целочисленная математика является ассоциативной (и имеет меньшую задержку, чем сложение FP, что помогает, если вы не собираетесь оптимизировать цикл с несколькими аккумуляторами или другими вещами, скрывающими задержку).float - & GT;int конверсия стоит примерно столько же, сколько и FPadd (фактически используя FP добавить ALU на процессорах Intel), так что в оптимизированном коде это не стоит.

Ответы на вопрос(6)

по которой преимущество Fortran в числовом коде, afaik, заключается в том, что язык облегчает его обнаружение.наложения спектров - компилятор может сказать, что умножаемые матрицы не используют одну и ту же память, что может помочь улучшить кэширование (нет необходимости быть уверенным, что результаты немедленно записываются обратно в «разделяемую» память). Вот почему C99 представилограничивать.

Тем не менее, в этом случае, я задаюсь вопросом, умудряется ли также код NumPy использовать некоторыеСпециальные инструкции что код C не является (поскольку разница кажется особенно большой).

Красивый код.

Ctypes должен пройти динамический перевод с C на Python и обратно, что добавляет некоторые накладные расходы. В Numpy большинство матричных операций выполняются полностью внутренне.

 20 авг. 2012 г., 04:44
Numpy не является оптимизированным кодом. Этоmakes use of оптимизированный код, например, ATLAS.
 20 авг. 2012 г., 06:59
@mohawkjohn Это оба.
 Charles Menguy04 мая 2012 г., 06:41
+1 за мою следующую прочитанную книгу!

Есть много способов оптимизировать матричное умножение. Например, порядок обхода матрицы влияет на шаблоны доступа к памяти, которые влияют на производительность.
Хорошее использование SSE - это еще один способ оптимизации, который, вероятно, использует NumPy.
Могут быть и другие способы, о которых разработчики NumPy знают, а я - нет.

Кстати, вы скомпилировали свой код C с optiomization?

Вы можете попробовать следующую оптимизацию для 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;
        }
    }
}
 04 мая 2012 г., 06:41
Хорошая реализация умножения побьет любой уровень оптимизации. Я предполагаю, что никакая оптимизация вообще не будет значительно хуже.
 Charles Menguy04 мая 2012 г., 06:40
Да, я пробовал с разными уровнями оптимизации при компиляции, но это не сильно изменило результат по сравнению с NumPy
 22 нояб. 2012 г., 11:26
Этот ответ делает много предположений о том, что делает Numpy. Тем не менее, он вряд ли выполняет какие-либо из них из коробки, вместо этого перенося работу в библиотеку BLAS, когда она доступна. Производительность умножения матриц сильно зависит от реализации BLAS.

тщательно настроенный метод BLAS для умножения матриц (см. Также:АТЛАС). Специфической функцией в этом случае является GEMM (для общего умножения матриц). Вы можете посмотреть оригинал, выполнив поискdgemm.f (это в Netlib).

Кстати, оптимизация выходит за рамки оптимизации компилятора. Выше Филипп упомянул Копперсмит - Виноград. Если я правильно помню, это алгоритм, который используется для большинства случаев умножения матриц в ATLAS (хотя комментатор отмечает, что это может быть алгоритм Штрассена).

Другими словами, вашmatmult Алгоритм является тривиальной реализацией. Есть более быстрые способы сделать то же самое.

 27 февр. 2013 г., 18:40
Кстати,np.show_config() показывает, с чем связано lapack / blas.
 20 апр. 2013 г., 03:53
Вы и Филипп делаете правильную точку (проблема в том, что реализация OP медленная), но я бы предположил, что NumPy использует алгоритм Штрассена или какой-то другой вариант, а не Coppersmith-Winograd, который имеет такие большие константы, что он обычно не полезно на практике.

используемый для реализации определенной функциональности, сам по себе является плохой мерой производительности. Часто использование более подходящего алгоритма является решающим фактором.

В вашем случае вы используете наивный подход к умножению матриц, как преподается в школе, который находится в O (n ^ 3). Тем не менее, вы можете сделать намного лучше для определенных видов матриц, например, квадратные матрицы, запасные матрицы и так далее.

Посмотрите наCoppersmith & # x2013; алгоритм Винограда (умножение квадратной матрицы в O (n ^ 2.3737)) для хорошей начальной точки быстрого умножения матриц. Также см. Раздел «Ссылки», в котором перечислены некоторые ссылки на еще более быстрые методы.


Для более естественного примера удивительного прироста производительности, попробуйте написать быстрыйstrlen() и сравните это с реализацией glibc. Если вам не удается победить, прочитайте glibcstrlen() Источник, он имеет довольно хорошие комментарии.

 20 авг. 2012 г., 04:44
Numpy не является оптимизированным кодом. Этоmakes use of оптимизированный код, например, ATLAS.
 20 авг. 2012 г., 06:59
@mohawkjohn Это оба.
 22 нояб. 2012 г., 11:26
Этот ответ делает много предположений о том, что делает Numpy. Тем не менее, он вряд ли выполняет какие-либо из них из коробки, вместо этого перенося работу в библиотеку BLAS, когда она доступна. Производительность умножения матриц сильно зависит от реализации BLAS.
 Charles Menguy04 мая 2012 г., 06:41
+1 за мою следующую прочитанную книгу!
 20 мая 2015 г., 03:19
+1 За использование биг-ой нотации и анализа (я всегда помню, наивный метод n ^ 3 против Штрассена alg с составляет около n ^ 2.8). Опять же, хороший способ проверить скорость alg - это не о языке.
 Charles Menguy04 мая 2012 г., 06:40
Да, я пробовал с разными уровнями оптимизации при компиляции, но это не сильно изменило результат по сравнению с NumPy
 04 мая 2012 г., 06:41
Хорошая реализация умножения побьет любой уровень оптимизации. Я предполагаю, что никакая оптимизация вообще не будет значительно хуже.
 28 авг. 2017 г., 08:08
Вероятно, в этом случае более важным является то, что наивная матрица OP не блокируется кешем и даже не транспонирует один из входов. Он циклически перебирает строки в одной матрице и столбцы в другой, когда они оба располагаются в основном порядке строк, поэтому он получает огромные потери в кеше. (Транспонирование - это O (n ^ 2), работа с фронтом, чтобы сделать вектор столбца строки *, точечные продукты делают последовательный доступ, что также позволяет им автоматически векторизовать с SSE / AVX / что угодно, если вы используете-ffast-math.)
Решение Вопроса

но источник на Github. Часть точечных продуктов реализована в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
@[email protected]_dot(char *ip1, npy_intp is1, char *ip2, npy_intp is2, char *op, npy_intp n,
           void *NPY_UNUSED(ignore))
{
    @[email protected] tmp = (@[email protected])0;
    npy_intp i;

    for (i = 0; i < n; i++, ip1 += is1, ip2 += is2) {
        tmp += (@[email protected])(*((@[email protected] *)ip1)) *
               (@[email protected])(*((@[email protected] *)ip2));
    }
    *((@[email protected] *)op) = (@[email protected]) 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);
}

Этот фактор развертывания, вероятно, был выбран после профилирования нескольких. Но одно теоретическое преимущество состоит в том, что между каждой точкой ветвления выполняется больше арифметических операций, и у компилятора и ЦП больше выбора в том, как оптимально планировать их, чтобы получить как можно больше конвейерной обработки команд.

 20 авг. 2012 г., 04:47
Я предполагаю, что ни один из этих алгоритмов на самом деле не используется для чисел с плавающей запятой, двойных чисел, одиночного комплекса или двойного комплекса. NumPy требует ATLAS, который имеет свои собственные версииdaxpy а такжеdgemm, Есть версии для поплавка и сложные; для целых чисел и тому подобного NumPy, вероятно, использует шаблон C, который вы связали.
 04 мая 2012 г., 07:57
Я снова был неправ, похоже, рутины в Numpy под/linalg/blas_lite.c называются. первыйdaxpy_ является развернутым внутренним циклом для точечных продуктов на числах с плавающей запятой, и основан на коде долгого времени назад. Проверьте комментарий там:"constant times a vector plus a vector. uses unrolled loops for increments equal to one. jack dongarra, linpack, 3/11/78. modified 12/3/93, array(1) declarations changed to array(*)"

Ваш ответ на вопрос