Почему умножение матриц выполняется быстрее с помощью numpy, чем с ctypes в Python?

Я не слишком знаком с Numpy, но исходник находится на 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
@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. Это было переведено с Фортрана на 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);
}

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

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

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

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

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

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

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

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

Numpy также является высоко оптимизированным кодом. Об этом есть эссе в книге Beautiful Code.

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

Люди, написавшие NumPy, очевидно, знают, что делают.

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

Кстати, вы скомпилировали свой код 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;
        }
    }
}

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

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