Комментарии 49
Интересно было бы сравнить по скорости с библиотекой MKT от того же Intel,
и например с Eigen еще.
На больших матрицах думаю будет небольшой проигрыш в пределах 5-10% (для 1 — потока). На малых размерах проигрыш будет в разы (причины этого указаны в статье).
А померять? Штрассен должен давать ускорение в 10-100 раз на больших матрицах и в 1000+ на ОЧЕНЬ больших и хорошие библиотеки обязаны его поддерживать в теории
UPD: вижу ваш ответ ниже. Странно.
UPD: вижу ваш ответ ниже. Странно.
А почему даже после введения правильного выравнивания используются инструкции для чтения/записи по невыровненному адресу? Разве от них не будет оверхеда? И вообще, мне кажется лучшим решением будет просто ввести требование по выравниванию непосредственно для входных данных.
Я уже достаточно давно разработываю проект по оптимизации различных алгоритмов (в основном по компьютерному зрению) при помощи SIMD инструкций. В частности я всегда старался выравнивать данные и использовать варианты с выровненными инструкциями (например, _mm_load_ps вместо _mm_loadu_ps), где только можно. Однако в последних поколениях процессоров Интел скорость работы этих инструкций одинаковая (при условии обращения к выровненному адресу конечно). Более того компилятор часто вместо _mm_load_ps подставляет ассемблерную инструкцию соответсвующую _mm_loadu_ps.
Ясно, спасибо за информацию, и за статью.
Быстрое гугление говорит, что разница между _mm_load_p и _mm_loadu_p на выравненных данных исчезла начиная с Nehalem (2008).
Я уже достаточно давно разработываю проект по оптимизации различных алгоритмов (в основном по компьютерному зрению) при помощи SIMD инструкций.
Обёртку на python ещё не завезли? :)
Очень информативно, а вы пробовали сравнить результаты для реализации "в лоб" у gcc и clang (-O3 -march=haswell -ffast-math)?
Я по правде говоря ожидал что оптимизации из первого шага и так должны быть сделаны компилятором (loop unswitch/rotate, const propogation?)
Вы не пробовали провести что-то подобное для матриц в столбцовом формате?
GPU просто созданы для быстрого перемножения матриц.
Зачем пытаться их обогнать на CPU? Ведь все равно не выйдет.
Зачем пытаться их обогнать на CPU? Ведь все равно не выйдет.
Я бы упомянул симпатичную математическую абстракцию часто используемую в подобных в подобных оптимизациях. Для матричного умножения справедливо что если мы поделим матрицу на блоки меньшего размера, то матрицы можно перемножать перемножая блоки
Это можно использовать либо явно (tiled версия aka блочное перемножение матриц) либо в рекурсивных схемах (где возникают всякие вкусные субкубические варианты типа алгоритма Штрассена). У Вас как я понимаю по сути tiled версия сильно заточенная под AVX — ширина блоков выбирается под ширину регистров, длина — под размер L1.
Это можно использовать либо явно (tiled версия aka блочное перемножение матриц) либо в рекурсивных схемах (где возникают всякие вкусные субкубические варианты типа алгоритма Штрассена). У Вас как я понимаю по сути tiled версия сильно заточенная под AVX — ширина блоков выбирается под ширину регистров, длина — под размер L1.
На сколько я видел в разных реализациях, нигде блочное перемножение не используется. Видимо потому, что это приводит к дополнительным копированиям памяти. А как следует из статьи выигрыш от локализации данных в кэше — до 10 раз. Выигрыш от подобных блочных алгоритмов будет проявляться только на матрицах циклопических размеров. А такие уже будет целесообразно считать на GPU.
Причём важна форма блоков! Типичная задача — перемножить большие матрицы в медленной памяти, цикл такой: разбиваем на блоки, пересылаем блоки аргументов во внутреннюю память, перемножаем, выгружаем блок результата. Как минимизировать пересылки? Оказывается, чем крупнее и квадратнее блок, тем лучше, поскольку объём работы по перемножению матриц IxK и KxJ пропорционален I*J*K, а объём пересылок — I*K + K*J + I*J.
Вспомнилось, как я году эдак в 2013 матрицы перемножал. Правда, там в ячейках были не float, а int64 (точнее, в A*B=C в A и B — инты, а в C должно было инт64 влезать). Размеры матричек 5000х5000. Подробнее о задаче на fastcomputing.org (вебархив, ибо проект ныне покойный). А вот табличка рекордов. В 40 секунд тогда упихал. Подробности реализации сейчас не вспомню, но вроде там Виноград-Штрассен, который где-то то ли на 200х200, то ли на 1000х1000 переключается уже на перемножение в лоб. Без Штрассена было тормознее, насколько помню. В перемножении в лоб тоже всякие SSE, помню разве что только один хак: мы умножаем не строчку A на строчку B (B транспонирована, конечно), а пару строчек A на пару строчек B, тем самым сокращая число чтений, потом все идет чисто на регистрах процессора, и на выходе получаем 4 значения для этих самых пар строчек.
MKL конечно платный, не всем подходит, но есть же MKL-DNN открытый.
Для мелких лучше использовать libxsmm. Он делает JIT компиляцию микрокернела под вашу платформу и бесплатный.
Для мелких лучше использовать libxsmm. Он делает JIT компиляцию микрокернела под вашу платформу и бесплатный.
Наверное, стоит добавить, что когда матриц много, то очень часто имеет смысл сначала определить порядок умножения матриц.
Очень полезная и хорошая статья. Насколько я знаю, MKL всегда работает на числе потоков, равном числу физических ядер. Насколько результаты статьи зависят от числа использованных виртуальных/физических потоков?
Цитата из введения: С целью ограничить объем изложения, я ограничился описанием однопоточного алгоритма для обычных процессоров. Тема многопоточности и алгоритмов для графических ускорителей явно заслуживает отдельной статьи.
MKL работает:
а) Если выставлен MKL_NUM_THREADS — то количество потоков ограничено этой переменной;
б) Если выставлен OMP_NUM_THREADS и не выставлена а) то ограничивается этой переменной;
По умолчанию считается что количество потоков = количеству ядер. Т.е если ОБЕ переменные не выставлены — будет работать на всем до чего дотянется.
а) Если выставлен MKL_NUM_THREADS — то количество потоков ограничено этой переменной;
б) Если выставлен OMP_NUM_THREADS и не выставлена а) то ограничивается этой переменной;
По умолчанию считается что количество потоков = количеству ядер. Т.е если ОБЕ переменные не выставлены — будет работать на всем до чего дотянется.
Т.е. для использования всех виртуальных ядер с учетом гипертрейдинга достаточно установить MKL_NUM_THREADS или OMP_NUM_THREADS в удвоенное количество физических ядер? Заранее спасибо.
Попробовал на dgemm с включенным гипертрейдингом. Уменьшить число потоков через set MKL_NUM_THREADS получается, а увеличить больше, чем число физических ядер — нет. Быстродействие упрямо показывает мах 50% загруженности cpu. Аналогично mkl_set_num_threads не дает загрузить число ядер больше физического. Функция mkl_get_max_threads возвращает значение равное числу физических ядер, а не виртуальных. Помнится, разработчики из Нижнего Новгорода подтверждали этот факт.
Тут решение уже за MKL. То, что вы видите загрузку в 50% — к реальности может(и скорее всего так и есть) не иметь никакого отношения, так как каждое из виртуальных ядер может полностью загрузить реальное. Т.е. с точки зрения операционки, одно ядро загружено полностью, на втором вообще ничего не запущено — загрузка 50%, а на деле единственное физическое ядро имеет полностью загруженные конвееры FPU.
MKL библиотека от производителя процессора, я сильно сомневаюсь чтобы они специально занижали производительность своих продуктов.
А в случае когда два потока работают на одном физическом ядре они разделяют как блоки FPU, так и кеш, что скажется на производительности скорее негативно, чем позитивно.
MKL библиотека от производителя процессора, я сильно сомневаюсь чтобы они специально занижали производительность своих продуктов.
А в случае когда два потока работают на одном физическом ядре они разделяют как блоки FPU, так и кеш, что скажется на производительности скорее негативно, чем позитивно.
Использование гипертридинга в нагруженных кернелах не приведет ни к чему хорошему.
50% загрузка видимо видна в диспетчере процессов Windows, там да. Логическое ядро отображается как нормальное, в связи с чем получается такое недопонимание.
50% загрузка видимо видна в диспетчере процессов Windows, там да. Логическое ядро отображается как нормальное, в связи с чем получается такое недопонимание.
Статья читается как какая-то магия. Я, как джавист из энтерпрайза, смог бы самостоятельно придти только к первому варианту решения, и если бы после этого встал вопрос повышения производительности, то только и смог бы придумать что «попробовать в несколько потоков выполнять» или «раскидать вычисления по облаку».
Делал подобное лет 10 назад, использовал практически все те же оптимизации, правда оптимизацией под L2/L3 уже не заморачивался, прирост не особо большой от этого получался. Результатами до сих пор пользуюсь — библиотека получилась небольшой, не надо что-то стороннее тянуть большое.
Но все-таки умножение матриц это достаточно простая вещь. Вот сделать то же самое для обращения матрицы (LUP разложения) и распараллелить обращение было посложнее.
Но все-таки умножение матриц это достаточно простая вещь. Вот сделать то же самое для обращения матрицы (LUP разложения) и распараллелить обращение было посложнее.
А что если перед умножением транспонировать матрицу B? Тогда по идее можно получить ускорение за счет того, что элементы матрицы B, которые участвуют в умножении на одну строку матрицы A станут соседями в памяти и будут лучше размещаться в кеше.
К сожалению нет. На первый взгляд действительно все проще: данные A и B лежат одинаково — только считай взаимное скалярное произведение их строчек.
Однако: максимальный размер микроядра получится 3x4, что дает нам (3 + 4)/(3*4) = ~0.58 загрузок на одну fma. Напомню, что при классической схеме с окном 6x16 получается (6 + 16)/(6*16) = ~0.23 загрузок на одну fma. Т.е. предложенная вами схема почти в 2.5 раза более требовательна к пропускной способности памяти. В принципе мои внутренние тесты это подтверждают.
Однако: максимальный размер микроядра получится 3x4, что дает нам (3 + 4)/(3*4) = ~0.58 загрузок на одну fma. Напомню, что при классической схеме с окном 6x16 получается (6 + 16)/(6*16) = ~0.23 загрузок на одну fma. Т.е. предложенная вами схема почти в 2.5 раза более требовательна к пропускной способности памяти. В принципе мои внутренние тесты это подтверждают.
Спасибо, классная статья. Не ожидал, что можно так сильно оптимизировать.
В начале статьи Вы критикуете малопонятность опенсурсных реализаций, ссылаясь в большей степени на ассемблер. Могу покритиковать Ваш код, так как он использует захардкоженные значения :) Думаю, именованные константы мне бы были более понятны, чем всякие 6, 24 и т. д…
В начале статьи Вы критикуете малопонятность опенсурсных реализаций, ссылаясь в большей степени на ассемблер. Могу покритиковать Ваш код, так как он использует захардкоженные значения :) Думаю, именованные константы мне бы были более понятны, чем всякие 6, 24 и т. д…
Как раз имею Core i5 4-го поколения (должно иметь приличные FMA), правда использовал для Matlab и для FFT (вроде как там по такой же логике польза от FMA) на комплексных числах. Вроде как прирост 4670 против Cel. 1005M достаточно большой был.
Ещё раньше использовал Athlon II X3 445, но для немного других задач. Там можно было явно записать решение задачи
dE/dt = f(E)*E + a*(laplas)E
в форме
dE/E = [f(E) — a*k2] dt.
То есть это обычный метод с функцией fft. При этом 3-ядерник грузился процентов на 60, но брат думает, что это тратилось второе ядро только на вывод графики.
Ещё раньше использовал Athlon II X3 445, но для немного других задач. Там можно было явно записать решение задачи
dE/dt = f(E)*E + a*(laplas)E
в форме
dE/E = [f(E) — a*k2] dt.
То есть это обычный метод с функцией fft. При этом 3-ядерник грузился процентов на 60, но брат думает, что это тратилось второе ядро только на вывод графики.
Мощно!
Почему-то в сети мало информации по умножению матриц, искал варианты реализации на Swift.
Почему-то в сети мало информации по умножению матриц, искал варианты реализации на Swift.
AlphaZero, говорят, вычислила новые алгоритмы умножения матриц, еще быстрее чем ранее, но подробностей не гуглится, только новость.
https://www.nature.com/articles/s41586-022-05172-4
Оригинальная статья.
Если коротко: они представили алгоритм умножения матриц как тензор (а):
На b изображен алгоритм Штрассена.
На с изображено разложение тензора для алгоритма Штрассена.
Далее с помощью нейросети стали искать разложения для различных размерностей матриц. Сам поиск представлен в виде игры:
На нулевом шаге есть тензор из примера а.
На каждом шаге нейросеть генерирует возможное разложение для текущего тензора, которое собирается в тензор и вычитается из текущего тензора.
Алгоритм работает пока не получится нулевой тензор, либо пока не будет достигнут лимит ходов.
Из всех предложенных нейросетью разложений собирается финальное разложение оригинального тензора.
Так, например, для умножения 4х4 на 4х4 в алгоритме Штрассена нужно 49 умножений, а они нашли вариант за 47 умножений.
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Умножение матриц: эффективная реализация шаг за шагом