Алгоритм умножения матриц

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

Непосредственное применение математического определения умножения матриц дает алгоритм, требующий времени порядка n. ³ полевые операции для умножения двух матриц размера n × n над этим полем ( Θ( n ³ ) в большой записи O ). Лучшие асимптотические границы времени, необходимого для умножения матриц, были известны со времен алгоритма Штрассена в 1960-х годах, но оптимальное время (то есть вычислительная сложность умножения матриц ) остается неизвестным. По состоянию на апрель 2024 г. , наилучшая заявленная оценка асимптотической сложности алгоритма умножения матриц равна O( n ^2.371552 ) время, данное Уильямсом , Сюй, Сюй и Чжоу. Это улучшает оценку O( n ^2.3728596 ) время, данное Алманом и Уильямсом. Однако этот алгоритм является галактическим из-за больших констант и не может быть реализован практически.

Итерационный алгоритм

Определение умножения матриц таково: если C = AB для размера n × m матрицы A и размера m × p матрицы B , то C является матрицей размера n × p с элементами

$c_{ij}=\sum _{k=1}^{m}a_{ik}b_{kj}.$

Исходя из этого, можно построить простой алгоритм, который циклически перебирает индексы i от 1 до n и j от 1 до p , вычисляя вышеуказанное с помощью вложенного цикла:

• Входные данные: матрицы A и B.
• Пусть C — новая матрица соответствующего размера.
• Для i от 1 до n :
  • Для j от 1 до p :
    • Пусть сумма = 0
    • Для k от 1 до m :
      • Установить сумму ← sum + A _ik × B _kj
    • Установить C _ij ← am
• Возврат С

Этот алгоритм требует времени Θ( nmp ) (в асимптотических обозначениях ). Обычное упрощение для анализа алгоритма состоит в том, чтобы предположить, что все входные данные представляют собой квадратные матрицы размера n × n , и в этом случае время работы равно Θ( n ³ ) , т. е. кубической по размеру размерности.

Поведение кэша

Иллюстрация порядка строк и столбцов

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

В частности, в идеализированном случае полностью ассоциативного кэша, состоящего из M байтов и b байтов на строку кэша (т.е. ⁠ M / b ⁠ строк кэша), приведенный выше алгоритм неоптимален для A и B, хранящихся в порядке строк. Когда n > ⁠ M / b ⁠ , каждая итерация внутреннего цикла (одновременный просмотр строки A и столбца B к элементу B. ) вызывает промах в кэше при доступе Это означает, что алгоритм несет Θ( n ³ ) в худшем случае кэш промахивается. По состоянию на 2010 год , скорость памяти по сравнению со скоростью процессоров такова, что промахи в кэше, а не фактические вычисления, доминируют над временем работы для матриц значительного размера.

Оптимальным вариантом итерационного алгоритма для A и B в макете с разбивкой по строкам является мозаичная версия, в которой матрица неявно делится на квадратные плитки размером √ M на √ M :

• Входные данные: матрицы A и B.
• Пусть C — новая матрица соответствующего размера.
• Выберите размер плитки T = Θ( √ M )
• Для I от 1 до n с шагом T :
  • Для J от 1 до p с шагом T :
    • Для K от 1 до m с шагом T :
      • Умножим A _{I : I + T , K : K + T} и B _{K : K + T , J : J + T} на C _{I : I + T , J : J + T} , то есть:
      • Для i от I до min( I + T , n ) :
        • Для j от J до min( J + T , p ) :
          • Пусть сумма = 0
          • Для k от K до min( K + T , m ) :
            • Установить сумму ← sum + A _ik × B _kj
          • Установить C _ij ← C _ij + сумма
• Возврат С

В идеализированной модели кэша этот алгоритм требует только Θ( ⁠ n ³ / b √ M ⁠ ) промахи в кэше; делитель b √ M на современных машинах составляет несколько порядков, так что фактические вычисления доминируют во времени выполнения, а не промахи в кэше.

Алгоритм «разделяй и властвуй»

Альтернативой итеративному алгоритму является «разделяй и властвуй» алгоритм умножения матриц . Это зависит от разделения блоков

$C={\begin{pmatrix}C_{11}&C_{12}\\C_{21}&C_{22}\\\end{pmatrix}},\,A={\begin{pmatrix}A_{11}&A_{12}\\A_{21}&A_{22}\\\end{pmatrix}},\,B={\begin{pmatrix}B_{11}&B_{12}\\B_{21}&B_{22}\\\end{pmatrix}},$

который работает для всех квадратных матриц, размеры которых являются степенями двойки, т. е. формы равны 2. ^н × 2 ^н для некоторых н . Матричный продукт теперь

${\begin{pmatrix}C_{11}&C_{12}\\C_{21}&C_{22}\\\end{pmatrix}}={\begin{pmatrix}A_{11}&A_{12}\\A_{21}&A_{22}\\\end{pmatrix}}{\begin{pmatrix}B_{11}&B_{12}\\B_{21}&B_{22}\\\end{pmatrix}}={\begin{pmatrix}A_{11}B_{11}+A_{12}B_{21}&A_{11}B_{12}+A_{12}B_{22}\\A_{21}B_{11}+A_{22}B_{21}&A_{21}B_{12}+A_{22}B_{22}\\\end{pmatrix}}$

который состоит из восьми умножений пар подматриц с последующим шагом сложения. Алгоритм «разделяй и властвуй» рекурсивно вычисляет меньшие умножения , используя скалярное умножение c ₁₁ = a ₁₁ b ₁₁ в качестве базового случая.

Сложность этого алгоритма как функция от n определяется рекуррентностью

$T(1)=\Theta (1);$

$T(n)=8T(n/2)+\Theta (n^{2}),$

с учетом восьми рекурсивных вызовов матриц размера n /2 и Θ( n ² ) для поэлементного суммирования четырех пар результирующих матриц. Применение основной теоремы для рекурсий «разделяй и властвуй» показывает, что эта рекурсия имеет решение Θ( n ³ ) , то же, что и итерационный алгоритм.

Неквадратные матрицы

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

• Входные данные: матрицы A размера n × m , B размера m × p .
• Базовый случай: если max( n , m , p ) ниже некоторого порога, используйте развернутую версию итерационного алгоритма.
• Рекурсивные случаи:

• Если max( n , m , p ) = n , разделите A горизонтально:

$C={\begin{pmatrix}A_{1}\\A_{2}\end{pmatrix}}{B}={\begin{pmatrix}A_{1}B\\A_{2}B\end{pmatrix}}$

• В противном случае, если max( n , m , p ) = p , разделите B вертикально:

$C=A{\begin{pmatrix}B_{1}&B_{2}\end{pmatrix}}={\begin{pmatrix}AB_{1}&AB_{2}\end{pmatrix}}$

• В противном случае max( n , m , p ) = m . Разделите A вертикально и B горизонтально:

$C={\begin{pmatrix}A_{1}&A_{2}\end{pmatrix}}{\begin{pmatrix}B_{1}\\B_{2}\end{pmatrix}}=A_{1}B_{1}+A_{2}B_{2}$

Поведение кэша

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

Количество промахов кэша, возникающих в результате этого алгоритма, на машине с M строками идеального кэша, каждая размером b байт, ограничено

$\Theta \left(m+n+p+{\frac {mn+np+mp}{b}}+{\frac {mnp}{b{\sqrt {M}}}}\right)$

Субкубические алгоритмы

Дополнительная информация: Вычислительная сложность умножения матриц.

Улучшение оценок показателя ω с течением времени для вычислительной сложности умножения матриц $O(n^{\omega })$.

Существуют алгоритмы, которые обеспечивают лучшее время работы, чем простые алгоритмы. Первым был открыт алгоритм Штрассена , разработанный Фолькером Штрассеном в 1969 году и часто называемый «быстрым умножением матриц». Он основан на способе умножения двух матриц размера 2 × 2 , который требует всего 7 умножений (вместо обычных 8) за счет нескольких дополнительных операций сложения и вычитания. Применение этого рекурсивного подхода дает алгоритм с мультипликативной стоимостью $O(n^{\log _{2}7})\approx O(n^{2.807})$. Алгоритм Штрассена более сложен, а численная устойчивость снижается по сравнению с наивным алгоритмом. но это быстрее в случаях, когда n > 100 или около того и появляется в нескольких библиотеках, таких как BLAS . Это очень полезно для больших матриц в точных областях, таких как конечные поля , где численная стабильность не является проблемой.

Поскольку алгоритм Штрассена фактически используется в практическом численном программном обеспечении и системах компьютерной алгебры, улучшение констант, скрытых в большой нотации O, имеет свои преимущества. Ниже приведена таблица, в которой сравниваются ключевые аспекты улучшенной версии, основанной на рекурсивном умножении блочных матриц 2x2 посредством умножения 7-блочных матриц. По-прежнему $n$ дает размеры матрицы и $M$ обозначает размер памяти.

Прогресс в рекурсивном умножении матриц 2x2 в стиле Штрассена

Год	Ссылка	#умножения матрицы за шаг	#добавления матрицы за шаг	итоговые арифметические операции	общая сложность ввода-вывода
1969	улицы	7	18	$7n^{\log _{2}7}-6n^{2}$	$6\left({\frac {{\sqrt {3}}n}{\sqrt {M}}}\right)^{\log _{2}7}\cdot M-18n^{2}+3M$
1971	Виноград	7	15	$6n^{\log _{2}7}-5n^{2}$	$5\left({\frac {{\sqrt {3}}n}{\sqrt {M}}}\right)^{\log _{2}7}\cdot M-15n^{2}+3M$
2017	Карштадт, Шварц	7	12	$5n^{\log _{2}7}-4n^{2}+3n^{2}\log _{2}n$	$4\left({\frac {{\sqrt {3}}n}{\sqrt {M}}}\right)^{\log _{2}7}\cdot M-12n^{2}+3n^{2}\cdot \log _{2}\left({\frac {{\sqrt {2}}n}{\sqrt {M}}}\right)+5M$

Известно, что алгоритм Штрассена с шагом блочной матрицы 2x2 требует не менее 7 умножений блочной матрицы. В 1976 году Проберт показал, что такой алгоритм требует как минимум 15 сложений (включая вычитания), однако скрытым предположением было то, что блоки и блочная матрица 2x2 представлены в одном и том же базисе. Карштадт и Шварц рассчитали в разных базисах и обменяли 3 сложения на менее затратные преобразования базисов. Они также доказали, что при использовании разных базисов нельзя опускаться ниже 12 сложений за шаг. В последующей работе Бениамини и др. применил этот трюк с изменением базы к более общим разложениям, чем блочные матрицы 2x2, и улучшил ведущую константу для их времени выполнения.

остается открытым вопрос, В теоретической информатике насколько хорошо алгоритм Штрассена можно улучшить с точки зрения асимптотической сложности . Показатель умножения матрицы , обычно обозначаемый $\omega$, — наименьшее действительное число, для которого любое $n\times n$ матрицу над полем можно перемножить, используя $n^{\omega +o(1)}$ полевые операции. Текущий лучший связанный $\omega$ является $\omega <2.371552$Уильямс , Сюй , Сюй и Чжоу. Этот алгоритм, как и все другие недавние алгоритмы в этом направлении исследований, является обобщением алгоритма Копперсмита-Винограда, предложенного Доном Копперсмитом и Шмуэлем Виноградом в 1990 году. Концептуальная идея этих алгоритмов аналогична алгоритму Штрассена: разработан способ умножения двух k × k -матриц с числом менее k. ³ умножения, и этот метод применяется рекурсивно. Однако постоянный коэффициент, скрытый обозначением Big O, настолько велик, что эти алгоритмы имеют смысл только для матриц, которые слишком велики для обработки на современных компьютерах.

Алгоритм Фрейвалдса — это простой алгоритм Монте-Карло , который по матрицам A , B и C проверяется в Θ( n ² ) время, если AB = C .

АльфаТензор

В 2022 году DeepMind представила AlphaTensor — нейронную сеть , которая использовала аналогию с однопользовательской игрой для изобретения тысяч алгоритмов умножения матриц, включая некоторые из них, ранее обнаруженные людьми, а некоторые — нет. Операции были ограничены некоммутативным основным полем (нормальная арифметика) и конечным полем. $\mathbb {Z} /2\mathbb {Z}$ (арифметика по модулю 2). Лучший найденный «практический» (явное разложение тензора матричного умножения низкого ранга) работал за O(n ^2.778 ). Нахождение разложений низкого ранга таких тензоров (и не только) NP-трудно; оптимальное умножение даже для матриц 3х3 остается неизвестным даже в коммутативном поле. Для матриц 4x4 AlphaTensor неожиданно обнаружил решение с 47 шагами умножения, что является улучшением по сравнению с 49, необходимыми для алгоритма Штрассена 1969 года, хотя и ограниченным арифметикой по модулю 2. Аналогично, AlphaTensor решал матрицы 5x5 за 96 шагов вместо 98 шагов Штрассена. Основываясь на неожиданном открытии существования таких улучшений, другие исследователи быстро смогли найти аналогичный независимый алгоритм 4x4 и отдельно доработали 96-шаговый алгоритм 5x5 Deepmind до 95 шагов в арифметике mod 2 и до 97 шагов. в обычной арифметике. Некоторые алгоритмы никогда раньше не были обнаружены, например (4, 5, 5) был улучшен до 76 с 80 в обычной арифметике и арифметике по модулю 2.

Параллельные и распределенные алгоритмы

Параллелизм с общей памятью

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

${\begin{pmatrix}A_{11}B_{11}+A_{12}B_{21}&A_{11}B_{12}+A_{12}B_{22}\\A_{21}B_{11}+A_{22}B_{21}&A_{21}B_{12}+A_{22}B_{22}\\\end{pmatrix}}$

могут выполняться независимо друг от друга, как и четыре суммирования (хотя алгоритму необходимо «объединить» умножения перед выполнением суммирования). Используя полный параллелизм задачи, можно получить алгоритм, который можно выразить в в стиле fork-join псевдокоде :

Процедура умножения( C , A , B ) :

• Базовый случай: если n = 1 , установите c ₁₁ ← a ₁₁ × b ₁₁ (или умножьте небольшую блочную матрицу).
• В противном случае выделите место для новой матрицы T формы n × n , тогда:
  • Разделите A на A ₁₁ , A ₁₂ , A ₂₁ , A ₂₂ .
  • Разделите B на B11 _​ , B12 _​ , B21 _​ , B22 _​ .
  • Разделить C на C11 _​ , C12 _​ , C21 _​ , C22 _​ .
  • Разделить Т на Т ₁₁ , Т ₁₂ , Т ₂₁ , Т ₂₂ .
  • Параллельное выполнение:
    • Вилка умножения( C ₁₁ , A ₁₁ , B ₁₁ ) .
    • Вилочное умножение( C ₁₂ , A ₁₁ , B ₁₂ ) .
    • Вилка умножения( C ₂₁ , A ₂₁ , B ₁₁ ) .
    • Вилочное умножение( C ₂₂ , A ₂₁ , B ₁₂ ) .
    • Вилка умножается( T ₁₁ , A ₁₂ , B ₂₁ ) .
    • Вилка умножается( T ₁₂ , A ₁₂ , B ₂₂ ) .
    • Вилка умножения( T ₂₁ , A ₂₂ , B ₂₁ ) .
    • Вилка умножения( T ₂₂ , A ₂₂ , B ₂₂ ) .
  • Присоединяйтесь (дождитесь завершения параллельных разветвлений).
  • добавить( С , Т ) .
  • Освободить Т. ​

Процедура add( C , T ) добавляет T в C поэлементно:

• Базовый случай: если n = 1 , установите c ₁₁ ← c ₁₁ + t ₁₁ (или сделайте короткий цикл, возможно, развернутый).
• В противном случае:
  • Разделить C на C11 _​ , C12 _​ , C21 _​ , C22 _​ .
  • Разделить Т на Т ₁₁ , Т ₁₂ , Т ₂₁ , Т ₂₂ .
  • Параллельно:
    • Вилка добавить( C ₁₁ , T ₁₁ ) .
    • Вилка добавить( C ₁₂ , T ₁₂ ) .
    • Вилка добавить( C ₂₁ , T ₂₁ ) .
    • Вилка добавить( C ₂₂ , T ₂₂ ) .
  • Присоединиться .

Здесь fork — это ключевое слово, которое сигнализирует, что вычисление может выполняться параллельно с остальной частью вызова функции, в то время как соединение ожидает завершения всех ранее «разветвленных» вычислений. Раздел достигает своей цели только путем манипулирования указателями.

Этот алгоритм имеет пути критическую длину Θ(log ² n ) шагов, то есть столько времени требуется на идеальной машине с бесконечным числом процессоров; следовательно, он имеет максимально возможное ускорение Θ ( n ³ /бревно ² n ) на любом реальном компьютере. Алгоритм непрактичен из-за затрат на связь, присущих перемещению данных во временную матрицу T и обратно , но более практичный вариант обеспечивает Θ( n ² ) ускорение, без использования временной матрицы.

Блочное умножение матриц. алгоритме каждый процессор отвечает за одну подматрицу C. В 2D - В 3D-алгоритме каждая пара подматриц из A и B назначается одному процессору. перемножаемая

Алгоритмы, избегающие общения, и распределенные алгоритмы

В современных архитектурах с иерархической памятью стоимость загрузки и хранения элементов входной матрицы обычно превышает стоимость арифметических операций. На одной машине это объем данных, передаваемых между ОЗУ и кешем, а на многоузловой машине с распределенной памятью - это объем, передаваемый между узлами; в любом случае это называется полосой пропускания связи . Наивный алгоритм, использующий три вложенных цикла, использует Ω( n ³ ) пропускная способность связи.

Алгоритм Кэннона , также известный как 2D-алгоритм , представляет собой алгоритм, позволяющий избежать связи , который разбивает каждую входную матрицу на блочную матрицу, элементы которой являются подматрицами размера √ M /3 на √ M /3 , где M — размер быстрой памяти. Затем наивный алгоритм применяется к блочным матрицам, вычисляя продукты подматриц полностью в быстрой памяти. Это уменьшает пропускную способность связи до O ( n ³ / √ M ) , что является асимптотически оптимальным (для алгоритмов, выполняющих Ω( n ³ ) расчет).

В распределенной среде с p процессорами, расположенными в двумерной сетке √ p на √ p , каждому процессору может быть назначена одна подматрица результата, и произведение может быть вычислено с каждым процессором, передающим O ( n ² / √ p ) слов, что является асимптотически оптимальным при условии, что каждый узел хранит минимум O ( n ² / п ) элементы. Это можно улучшить с помощью 3D-алгоритма, который размещает процессоры в сетке 3D-куба, назначая каждый продукт двух входных подматриц одному процессору. Затем генерируются результирующие подматрицы путем сокращения каждой строки. Этот алгоритм передает O ( n ² / п ^2/3 ) слов на процессор, что асимптотически оптимально. Однако для этого требуется репликация каждого элемента входной матрицы p ^1/3 раз, поэтому требуется коэффициент p ^1/3 больше памяти, чем необходимо для хранения входных данных. Этот алгоритм можно комбинировать с алгоритмом Штрассена для дальнейшего сокращения времени выполнения. Алгоритмы «2.5D» обеспечивают постоянный компромисс между использованием памяти и пропускной способностью связи. В современных распределенных вычислительных средах, таких как MapReduce , разработаны специализированные алгоритмы умножения.

Алгоритмы для сеток

Умножение матриц выполняется за 2n-1 шагов для двух матриц n×n в сетке с перекрестными связями.

Существует множество алгоритмов умножения на сетках . Для умножения двух n × n на стандартной двумерной сетке с использованием алгоритма 2D Кэннона можно выполнить умножение за 3 n -2 шага, хотя при повторных вычислениях это число сокращается вдвое. Стандартный массив неэффективен, поскольку данные из двух матриц поступают не одновременно и их необходимо дополнять нулями.

Результат еще быстрее на двухслойной сетке с перекрестными проволоками, где всего 2 n -1 шагов. требуется Производительность еще больше улучшается при повторных вычислениях, что приводит к 100% эффективности. Массив с перекрестной сеткой можно рассматривать как частный случай непланарной (т.е. многослойной) структуры обработки.

См. также

• Вычислительная сложность математических операций
• Вычислительная сложность умножения матриц
• Алгоритм CYK § Алгоритм Валианта
• Умножение цепочки матриц
• Умножение разреженной матрицы на вектор
• Метод четырех русских

Дальнейшее чтение

• Буттари, Альфредо; Лангу, Жюльен; Курзак, Якуб; Донгарра, Джек (2009). «Класс параллельных алгоритмов тайловой линейной алгебры для многоядерных архитектур». Параллельные вычисления . 35 : 38–53. arXiv : 0709.1272 . дои : 10.1016/j.parco.2008.10.002 . S2CID   955 .
• Гото, Казусигэ; ван де Гейн, Роберт А. (2008). «Анатомия высокопроизводительного умножения матриц». Транзакции ACM в математическом программном обеспечении . 34 (3): 1–25. CiteSeerX   10.1.1.140.3583 . дои : 10.1145/1356052.1356053 . S2CID   9359223 .
• Ван Зи, Филд Г.; ван де Гейн, Роберт А. (2015). «BLIS: платформа для быстрого создания экземпляров функциональности BLAS». Транзакции ACM в математическом программном обеспечении . 41 (3): 1–33. дои : 10.1145/2764454 . S2CID   1242360 .
• Как оптимизировать GEMM