Почему MATLAB так быстр в умножении матриц?

Я делаю некоторые тесты с CUDA, C ++, C #, Java и использую MATLAB для проверки и генерации матрицы. Когда я выполняю умножение матриц с помощью MATLAB, 2048x2048и даже большие матрицы умножаются практически мгновенно.

             1024x1024   2048x2048   4096x4096
             ---------   ---------   ---------
CUDA C (ms)      43.11      391.05     3407.99
C++ (ms)       6137.10    64369.29   551390.93
C# (ms)       10509.00   300684.00  2527250.00
Java (ms)      9149.90    92562.28   838357.94
MATLAB (ms)      75.01      423.10     3133.90

Только CUDA конкурентоспособна, но я думал, что по крайней мере C ++ будет несколько ближе, а не в 60 раз медленнее. Я также не знаю, что думать о результатах C #. Алгоритм такой же , как C ++ и Java, но гигантский скачок 2048от 1024.

Как MATLAB так быстро выполняет умножение матриц?

Код C ++:

float temp = 0;
timer.start();
for(int j = 0; j < rozmer; j++)
{
    for (int k = 0; k < rozmer; k++)
    {
        temp = 0;
        for (int m = 0; m < rozmer; m++)
        {
            temp = temp + matice1[j][m] * matice2[m][k];
        }
        matice3[j][k] = temp;
    }
}
timer.stop();

Вот мои результаты с использованием MATLAB R2011a + Parallel Computing Toolbox на машине с Tesla C2070:

>> A = rand(1024); gA = gpuArray(A);
% warm up by executing the operations a couple of times, and then:
>> tic, C = A * A; toc
Elapsed time is 0.075396 seconds.
>> tic, gC = gA * gA; toc
Elapsed time is 0.008621 seconds.

MATLAB использует высоко оптимизированные библиотеки для умножения матриц, поэтому умножение матриц в MATLAB так быстро. gpuArrayВерсия использует МАГМА .

Обновление с использованием R2014a на машине с Tesla K20c, а также новые timeitи gputimeitфункции:

>> A = rand(1024); gA = gpuArray(A);
>> timeit(@()A*A)
ans =
    0.0324
>> gputimeit(@()gA*gA)
ans =
    0.0022

Обновление с использованием R2018b на машине WIN64 с 16 физическими ядрами и Tesla V100:

>> timeit(@()A*A)
ans =
    0.0229
>> gputimeit(@()gA*gA)
ans =
   4.8019e-04

(Примечание: в какой-то момент (я забыл, когда именно) gpuArrayпереключился с MAGMA на cuBLAS - MAGMA все еще используется для некоторых gpuArrayопераций)

Этот вопрос повторяется и на него следует ответить более четко, чем «MATLAB использует высоко оптимизированные библиотеки» или «MATLAB использует MKL» один раз в переполнении стека.

История:

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

Я не специалист по истории, но, видимо, тогда все просто переписали его версию на Фортране с помощью простых циклов. Затем пришла некоторая стандартизация с идентификацией «ядер» (базовых процедур), которые необходимы для решения большинства задач линейной алгебры. Эти основные операции были затем стандартизированы в спецификации под названием: «Основные подпрограммы линейной алгебры (BLAS)». Инженеры могли бы тогда назвать эти стандартные, хорошо проверенные подпрограммы BLAS в своем коде, делая их работу намного проще.

BLAS:

BLAS эволюционировал от уровня 1 (первая версия, которая определяла скалярно-векторные и векторно-векторные операции) до уровня 2 (операции вектор-матрица) до уровня 3 (операции матрица-матрица) и предоставлял все больше и больше «ядер», настолько стандартизированных и другие основные операции линейной алгебры. Оригинальные реализации FORTRAN 77 по-прежнему доступны на веб-сайте Netlib .

На пути к лучшей производительности:

Таким образом, за прошедшие годы (особенно между выпусками уровня 1 и уровня 2: в начале 80-х годов) аппаратное обеспечение изменилось с появлением векторных операций и иерархий кэша. Эти изменения позволили существенно повысить производительность подпрограмм BLAS. Различные поставщики затем пришли с их реализацией подпрограмм BLAS, которые были все более и более эффективными.

Я не знаю всех исторических реализаций (тогда я не родился и не был ребенком), но в начале 2000-х вышли две наиболее заметные из них: Intel MKL и GotoBLAS. Ваш Matlab использует Intel MKL, который является очень хорошим, оптимизированным BLAS, и это объясняет великолепную производительность, которую вы видите.

Технические подробности по умножению матриц:

Так почему же Matlab (MKL) так быстр dgemm(общее умножение матрицы на матрицу с двойной точностью)? Проще говоря: потому что он использует векторизацию и хорошее кэширование данных. В более сложных терминах: см. Статью Джонатана Мура.

По сути, когда вы выполняете умножение в предоставленном вами коде C ++, вы совсем не дружественны к кешу. Поскольку я подозреваю, что вы создали массив указателей на массивы строк, ваш доступ во внутреннем цикле к k-му столбцу "matice2": matice2[m][k]очень медленный. Действительно, когда вы matice2[0][k]получаете доступ , вы должны получить k-й элемент массива 0 вашей матрицы. Затем на следующей итерации вы должны получить доступ matice2[1][k], который является k-м элементом другого массива (массив 1). Затем на следующей итерации вы получаете доступ к еще одному массиву и т. Д. Так как вся матрица matice2не может поместиться в старшие кеши (это 8*1024*1024большие байты), программа должна извлечь нужный элемент из основной памяти, потеряв много время.

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

timer.start();
float temp = 0;
//transpose matice2
for (int p = 0; p < rozmer; p++)
{
    for (int q = 0; q < rozmer; q++)
    {
        tempmat[p][q] = matice2[q][p];
    }
}
for(int j = 0; j < rozmer; j++)
{
    for (int k = 0; k < rozmer; k++)
    {
        temp = 0;
        for (int m = 0; m < rozmer; m++)
        {
            temp = temp + matice1[j][m] * tempmat[k][m];
        }
        matice3[j][k] = temp;
    }
}
timer.stop();

Таким образом, вы можете видеть, как простота размещения кэша значительно увеличила производительность вашего кода. Теперь реальные dgemmреализации используют это на очень обширном уровне: они выполняют умножение на блоки матрицы, определяемые размером TLB (буфер с переводом в сторону, короче говоря: что может эффективно кэшироваться), чтобы они передавались в процессор именно тот объем данных, который он может обработать. Другой аспект - векторизация, они используют векторизованные инструкции процессора для оптимальной пропускной способности команд, чего вы не можете сделать из своего кроссплатформенного кода C ++.

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

Вот почему . MATLAB не выполняет наивное матричное умножение, зацикливаясь на каждом отдельном элементе, как вы делали это в своем коде C ++.

Конечно, я предполагаю, что вы просто использовали C=A*Bвместо того, чтобы писать функцию умножения самостоятельно.

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

Вы также можете посмотреть статью Гото и Ван Де Гейна «Анатомия высокопроизводительного матричного умножения» по адресу http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.140.1785&rep=rep1&type=pdf.

Ответ заключается в том, что библиотеки LAPACK и BLAS делают MATLAB ослепительно быстрым при матричных операциях, а не каким-либо частным кодом, созданным людьми в MATLAB.

Используйте библиотеки LAPACK и / или BLAS в своем коде C ++ для матричных операций, и вы должны получить производительность, аналогичную MATLAB. Эти библиотеки должны быть свободно доступны в любой современной системе, а их части были разработаны в течение десятилетий в научных кругах. Обратите внимание, что существует несколько реализаций, включая некоторые с закрытым исходным кодом, такие как Intel MKL .

Обсуждение того, как BLAS достигает высокой производительности , доступно здесь.

Кстати, я испытываю серьезную боль при вызове библиотек LAPACK непосредственно из c (но оно того стоит). Вы должны прочитать документацию ОЧЕНЬ точно.

Делая матричное умножение, вы используете наивный метод умножения, который занимает время O(n^3).

Существует алгоритм умножения матриц, который принимает O(n^2.4). Это означает, что для n=2000вашего алгоритма требуется ~ 100 раз больше вычислений, чем для лучшего алгоритма.
Вы действительно должны проверить страницу википедии для умножения матриц для получения дополнительной информации об эффективных способах ее реализации.

Я полагаю, что в зависимости от вашей версии Matlab она уже использует ваш графический процессор.

Еще одна вещь; Matlab отслеживает многие свойства вашей матрицы; по его диагонали, по Герметику и т. д., и специализирует свои алгоритмы на ее основе. Может быть, это специализация, основанная на нулевой матрице, которую вы передаете, или что-то в этом роде? Может быть, это кеширование повторных вызовов функций, которое портит ваши настройки? Возможно, это оптимизирует повторяющиеся неиспользованные матричные продукты?

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

MATLAB использует высоко оптимизированную реализацию LAPACK от Intel, известную как Intel Math Kernel Library (Intel MKL), в частности, функцию dgemm . Скорость Эта библиотека использует преимущества процессорных функций, включая инструкции SIMD и многоядерные процессоры. Они не документируют, какой конкретный алгоритм они используют. Если бы вы вызывали Intel MKL из C ++, вы должны увидеть аналогичную производительность.

Я не уверен, какую библиотеку MATLAB использует для умножения GPU, но, вероятно, что-то вроде nVidia CUBLAS .

Общий ответ на вопрос «Почему Matlab работает быстрее, чем в других программах», заключается в том, что в Matlab много встроенных оптимизированных функций.

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

Это можно интерпретировать двумя способами:

1) Обычный / теоретический способ: Matlab не значительно быстрее, вы просто делаете неверный тест

2) Реалистичный способ: для этого Matlab быстрее на практике, потому что языки как c ++ слишком легко используются неэффективными способами.

Резкий контраст связан не только с удивительной оптимизацией Matlab (как уже говорилось во многих других ответах), но и с тем, как вы формулировали матрицу как объект.

Кажется, вы сделали матрицу списком списков? Список списков содержит указатели на списки, которые затем содержат ваши матричные элементы. Местоположения содержащихся списков назначаются произвольно. Поскольку вы перебираете свой первый индекс (номер строки?), Время доступа к памяти очень значительно. Для сравнения, почему бы вам не попробовать реализовать матрицу в виде одного списка / вектора, используя следующий метод?

#include <vector>

struct matrix {
    matrix(int x, int y) : n_row(x), n_col(y), M(x * y) {}
    int n_row;
    int n_col;
    std::vector<double> M;
    double &operator()(int i, int j);
};

И

double &matrix::operator()(int i, int j) {
    return M[n_col * i + j];
}

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

Я прошу вас рассчитать время, чтобы результат был сопоставим с тем, что вы имели ранее (на той же машине). Сравнение покажет, насколько значительным может быть время доступа к памяти!

Это медленно в C ++, потому что вы не используете многопоточность. По сути, если A = BC, где все они являются матрицами, первая строка A может быть вычислена независимо от 2-й строки и т. Д. Если A, B и C все n на n матриц, вы можете ускорить умножение на фактор n ^ 2, как

a_ {i, j} = sum_ {k} b_ {i, k} c_ {k, j}

Если вы используете, скажем, Eigen [ http://eigen.tuxfamily.org/dox/GettingStarted.html ], многопоточность встроена, а количество потоков настраивается.

Потому что MATLAB - это язык программирования, изначально разработанный для числовой линейной алгебры (матричные манипуляции), в котором есть библиотеки, специально разработанные для умножения матриц. И теперь MATLAB также может использовать для этого графические процессоры (графические процессоры) .

И если мы посмотрим на ваши результаты вычислений:

             1024x1024   2048x2048   4096x4096
             ---------   ---------   ---------
CUDA C (ms)      43.11      391.05     3407.99
C++ (ms)       6137.10    64369.29   551390.93
C# (ms)       10509.00   300684.00  2527250.00
Java (ms)      9149.90    92562.28   838357.94
MATLAB (ms)      75.01      423.10     3133.90

тогда мы можем видеть, что не только MATLAB так быстр в умножении матриц: CUDA C (язык программирования от NVIDIA) дает лучшие результаты, чем MATLAB. CUDA C также имеет библиотеки, специально разработанные для умножения матриц, и использует графические процессоры.

Краткая история MATLAB

Клив Молер, председатель факультета компьютерных наук Университета Нью-Мексико, начал разработку MATLAB в конце 1970-х годов. Он разработал его, чтобы дать своим студентам доступ к LINPACK (программная библиотека для выполнения числовой линейной алгебры) и EISPACK(это библиотека программного обеспечения для численных расчетов линейной алгебры) без необходимости изучать Фортран. Вскоре он распространился на другие университеты и нашел сильную аудиторию в сообществе прикладной математики. Джек Литтл, инженер, познакомился с ним во время визита, который Мёлер совершил в Стэнфордском университете в 1983 году. Признавая его коммерческий потенциал, он присоединился к Мёлеру и Стиву Бангерту. Они переписали MATLAB на C и основали MathWorks в 1984 году, чтобы продолжить его развитие. Эти переписанные библиотеки были известны как JACKPAC. В 2000 году MATLAB был переписан для использования более нового набора библиотек для работы с матрицами, LAPACK (стандартная библиотека программного обеспечения для числовой линейной алгебры).

Источник

Что такое CUDA C

CUDA C также использует библиотеки, специально разработанные для умножения матриц, такие как OpenGL (Open Graphics Library). Он также использует GPU и Direct3D (на MS Windows).

Платформа CUDA предназначена для работы с такими языками программирования, как C, C ++ и Fortran. Такая доступность облегчает специалистам по параллельному программированию использование ресурсов графического процессора, в отличие от предыдущих API, таких как Direct3D и OpenGL , которые требовали дополнительных навыков в графическом программировании. Кроме того, CUDA поддерживает фреймворки программирования, такие как OpenACC и OpenCL .

Пример процесса обработки CUDA:

1. Копировать данные из основной памяти в память графического процессора
2. CPU запускает вычислительное ядро ​​на GPU
3. Ядра CUDA графического процессора выполняют ядро ​​параллельно
4. Скопируйте полученные данные из памяти графического процессора в основную память

Сравнение скорости выполнения CPU и GPU

Мы провели тест, в котором мы измерили количество времени, которое потребовалось для выполнения 50 временных шагов для размеров сетки 64, 128, 512, 1024 и 2048 на процессоре Intel Xeon X5650, а затем с помощью графического процессора NVIDIA Tesla C2050.

Для размера сетки 2048 алгоритм показывает уменьшение времени вычислений в 7,5 раз с более минуты на ЦП до менее 10 секунд на ГП. График масштаба журнала показывает, что ЦП на самом деле быстрее для небольших размеров сетки. Однако по мере развития и развития технологии решения для графических процессоров все в большей степени способны справляться с небольшими проблемами, и эта тенденция, как мы ожидаем, сохранится.

Источник

Из введения в Руководство по программированию CUDA C:

Вследствие ненасытного рыночного спроса на трехмерную графику высокой четкости в реальном времени программируемый графический процессор или графический процессор превратился в высокопараллельный многопоточный многоядерный процессор с огромной вычислительной мощностью и очень высокой пропускной способностью памяти, что иллюстрируется Figure 1и Figure 2.

Рисунок 1. Операции с плавающей запятой в секунду для CPU и GPU

Рисунок 2 . Пропускная способность памяти для процессора и графического процессора

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

Рисунок 3 . Графический процессор выделяет больше транзисторов для обработки данных

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

Параллельная обработка данных отображает элементы данных в потоки параллельной обработки. Многие приложения, обрабатывающие большие наборы данных, могут использовать модель параллельного программирования для ускорения вычислений. При 3D-рендеринге большие наборы пикселей и вершин отображаются в параллельные потоки. Аналогично, приложения обработки изображений и мультимедиа, такие как постобработка визуализированных изображений, кодирование и декодирование видео, масштабирование изображений, стереозрение и распознавание образов, могут отображать блоки изображений и пиксели в потоки параллельной обработки. Фактически, многие алгоритмы за пределами области визуализации и обработки изображений ускоряются параллельной обработкой данных, от общей обработки сигналов или физического моделирования до вычислительных финансов или вычислительной биологии.

Источник

Расширенное чтение

• GPU (графический процессор)
• MATLAB
• Руководство по программированию CUDA C
• Использование графических процессоров в MATLAB

• Основные подпрограммы линейной алгебры (BLAS)

• Анатомия высокопроизводительного матричного умножения , от Казушиге Гото и Роберта А. Ван Де Гейна

Некоторые интересные лица

Я написал C ++ матричное умножение, которое так же быстро, как и Matlab, но позаботилось об этом. (До этого Matlab использовал для этого графические процессоры).

Цитата из этого ответа .