Подтвердить что ты не робот

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

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

  • Чистая реализация python: здесь нет сюрпризов.
  • Реализация Numpy с использованием numpy.dot(a, b)
  • Взаимодействие с C с использованием модуля ctypes в Python.

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

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

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

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

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

benchmark

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

4b9b3361

Ответ 1

Я не слишком хорошо знаком с 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
@[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 также содержит код для выбора реализации определенных операций, включая "точку", из разных базовых реализаций. Например, он может использовать библиотеку 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);
}

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

Ответ 2

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

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

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

Ответ 3

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

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

Посмотрите на алгоритм Coppersmith-Winograd (умножение квадратной матрицы в O (n ^ 2.3737)) для хорошей отправной точки на быстрых матричное умножение. Также см. Раздел "Ссылки", в котором перечислены некоторые указатели на более быстрые методы.


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

Ответ 4

Ребята, которые писали NumPy, очевидно, знают, что они делают.

Существует множество способов оптимизации умножения матрицы. Например, упорядочить перемещение матрицы влияет на шаблоны доступа к памяти, которые влияют на производительность.
Хорошее использование 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;
        }
    }
}

Ответ 5

Numpy также является высоко оптимизированным кодом. В книге "Красивый код" есть эссе о его частях.

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

Ответ 6

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

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