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

Нерешённые проблемы информатики: Насколько быстр может быть алгоритм умножения матриц?

Поскольку умножение матриц является центральной операцией во многих численных алгоритмах, много усилий было вложено в повышение эффективности алгоритма умножения матриц. Приложения алгоритма умножения матриц в вычислительных задачах найдены во многих областях, включая научные вычисления^[en] и распознавания образов, а также во вроде бы не имеющих отношение к матрицам задачах, таких как подсчёт путей через граф^[1]. Было разработано много различных алгоритмов для умножения матриц на оборудовании различного типа, включая параллельные и распределённые системы, где вычисления распределены на несколько процессоров (и, может быть, по сети).

Прямое использование математического определения умножения матриц даёт алгоритм, который работает за время порядка $n^{3}$ операций поля для умножения двух $n\times n$ матриц над этим полем (или $\Theta (n^{3})$ в нотации «O» большое). Улучшенные асимптотические границы по времени были известны с момента появления алгоритма Штрассена в 1960-х годах, но оптимальное время остаётся неизвестным (то есть, неизвестна сложность задачи). К концу 2020 года алгоритм умножения матриц с лучшей асимптотической сложностью, работающий за время $O(n^{2{,}3728596})$ , был дан Джозефом Алманом и Вирджинией Василевска Уильямс^[en]^[2], однако этот алгоритм галактического масштаба, то есть только для данных галактического размера, поскольку содержит огромные константы и не может быть реализован на практике.

Итеративный алгоритм править

По определению умножения матриц для $n\times m$ матрицы A и по $m\times p$ матрицы B произведением $C=AB$ является $n\times 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:
  - Положим $sum=0$
  - Для k от 1 до m:
    - Положим $sum\leftarrow sum+A_{ik}\times B_{kj}$
  - Положим $C_{ij}\leftarrow sum$
Возвращаем C

Этот алгоритм работает за время $\Theta (nmp)$ (в асимптотических обозначениях) ^[1]. Обычно для упрощения анализа алгоритма предполагается, что входными матрицами являются квадратные матрицы размера $n\times n$ , и в этом случае время работы составляет $\Theta (n^{3})$ , то есть время зависит кубически от размера матриц^[3].

Поведение кэша править

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

Три цикла при итеративном умножении матриц можно произвольным образом переставлять друг с другом без влияния на правильность алгоритма или асимптотическое время работы. Однако, порядок циклов может влиять на практические характеристики доступа памяти^[en] и на алгоритм использования кэша^[1]. Какой порядок вычисления лучше, зависит от того, как хранятся входные матрицы — в построчном порядке^[en], постолбцовом порядке, или в смешении этих порядков.

В частности, в идеальном случае полностью ассоциативного кэша, состоящего из M байт и b байт на строку кэша (то есть, M/b строк кэша), вышеприведённый алгоритм является подоптимальным для A и B, хранящихся построчно. Если $n>{\tfrac {M}{b}}$ , любая итерация внутреннего цикла (одновременный проход по строкам A и столбцам B) приводит к промахам кэша при выборке элемента матрицы B. Это означает, что при работе алгоритма будет $\Theta (n^{3})$ промахов кэша в худшем случае^[en]. К 2010 году для больших матриц скорость доступа к памяти являлась доминирующим фактором, определяющим время работы алгоритма, а не скорость процессора ^[4].

Оптимальным вариантом итерационного алгоритма для матриц A и B при построчном хранении является версия с разбиением на блоки^[en], где матрицы в неявном виде разбиты на квадратные части размера ${\sqrt {M}}\times {\sqrt {M}}$ ^[4]^[5]:

Вход: матрицы A и B
Пусть C будет новой матрицей нужного размера
Выбираем часть размера $T=\Theta ({\sqrt {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)$ :
        Положим sum = 0
        
        Для k от K до $\min(K+T,m)$ :
        Положим $sum\leftarrow sum+A_{ik}\times B_{kj}$
        
        Положим $C_{ij}\leftarrow C_{ij}+sum$
Возвращаем C

В случае идеального кэша алгоритм приводит только к $\Theta ({\tfrac {n^{3}}{b{\sqrt {M}}}})$ промахам. Делитель $b{\sqrt {M}}$ составляет несколько порядков величины для современных машин, так что вычисления доминируют во времени работы, а не промахи кэша^[4].

Алгоритм Разделяй-и-властвуй править

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

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^{n}\times 2^{n}$ для некоторого n. Произведение матриц тогда равно

{\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_{11}=a_{11}b_{11}$ как в базовом случае.

Сложность этого алгоритма как функция от n задаётся рекурсией^[3]

T(1)=\Theta (1)

;

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

,

принимающей во внимание восемь рекурсивных вызовов на матрицах размера n/2 и величину $\Theta (n^{2})$ для суммирования четырёх пар полученных матриц. Применяя основную теорему о рекуррентных соотношениях, получим, что эта рекурсия имеет решение $\Theta (n^{3})$ , ту же самую сложность, что и для итеративного алгоритма^[3].

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

Вариант этого алгоритма, который работает для матриц произвольного размера и на практике быстрее^[4], разбивает матрицы на две, а не на четыре подматрицы, что продемонстрировано ниже^[6]. Разбиение матрицы теперь означает разделение её на две одинаковые или близкие к одинаковым части, если размер нечётен.

Вход: матрица A размера $n\times m$ и матрица B размера $m\times 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}

Поведение кэша править

Количество промахов кэша рекурсивного умножения матриц та же самая, что и у версии aлгоритма с разбиением на блоки^[en], но, в отличие от этого алгоритма, рекурсивный алгоритм нечувствителен к кешированию^[en]^[6] — не нужно никакого настоечного параметра для получения оптимального поведения кэша и он работает хорошо в многозадачном окружении, где размер кэша меняется динамически, поскольку другие процессы тоже используют кэш^[4]. (Простой итеративный алгоритм нечувствителен к кэшированию также, но на практике много медленнее, если матрицы не подогнаны под алгоритм.)

Число промахов кэша при этом алгоритме на машинах с M линиями идеального кэша, каждая из которых имеет размер в b байт, ограничена величиной^[7]

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

Подкубические алгоритмы править

Улучшение оценки показателя

\omega

от времени вычислительной сложности умножения матриц

O(n^{\omega })

.

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

Открытым вопросом теоретической информатики является вопрос, насколько можно улучшить алгоритм Штрассена. Показатель умножения матриц $\omega$ — это наименьшее вещественное число, для которого произведение любых двух $n\times n$ матриц над полем может быть вычислено за $n^{\omega +o(1)}$ операций в поле. Текущая лучшая граница $\omega$ равна 2,3728596 (алгоритм Джошуа Алманаи Вирджинии Вассилевска^[en]^[2]. Этот алгоритм, подобно всем другим недавно разработанным алгоритмам в этом направлении исследований, является обобщением алгоритма Копперсмита — Винограда, который представили Дон Копперсмит и Шмуэль Виноград^[en] в 1990 году, и который имеет асимптотическую сложность $O(n^{2,376})$ . Концептуальная идея этих алгоритмов аналогична алгоритму Штрассена — способ разрабатывается для умножения двух $k\times k$ матриц за менее чем $k^{3}$ умножений и эта техника применяется рекурсивно. Однако, константы, спрятанные в нотации «O большое» так велики, что эти алгоритмы целесообразно применять только для матриц, которые слишком велики, чтобы их можно было обрабатывать на существующих компьютерах^[10]^[11].

Поскольку любой алгоритм умножения двух $n\times n$ матриц должен обработать все $2n^{2}$ элемента, имеется асимптотическая нижняя граница числа операций $\Omega (n^{2})$ . Рэн Раз доказал нижнюю границу в $\Omega (n^{2}log(n))$ ограниченных коэффициентов арифметических цепей над вещественными или комплексными числами^[12].

Кон с соавторами переложил методы, такие как алгоритмы Штрассена и Копперсмита — Винограда в полностью другой контекст теории групп путём использования троек подмножеств конечных групп, которые удовлетворяют свойством дизъютктности, называемое свойством тройного произведения^[en] (СТП, англ. triple product property, TPP). Они показали, что если семейства веночного произведения^[en] абелевых групп с симметричными группами образуют семейства троек со свойствами, аналогичными СТП, то существуют алгоритмы умножения матриц с фактически квадратичной сложностью^[13]^[14]. Большинство исследователей верят, что это верно вообще^[11]. Однако Алон, Шпилька и Хрис Уманс недавно показали, что некоторые из гипотез о быстром умножении матриц несовместимы с другой правдоподобной гипотезой, гипотезе о подсолнухе^[en]^[15].

Алгоритм Фрейвалдса — это простой алгоритм Монте-Карло, который для заданных матриц $A,B$ и C проверяет, что $AB=C$ , за время $\Theta (n^{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 псевдокода^[16]:

Процедура multiply(C, A, B):

Базовый случай: если $n=1$ , положим $c_{11}\leftarrow a_{11}\times b_{11}$ (умножаем маленькие блочные матрицы).
В противном случае распределяем место для новой матрицы T размером $n\times n$ , затем:
- Разбиваем A на $A_{11},A_{12},A_{21},A_{22}$ .
- Разбиваем B на $B_{11},B_{12},B_{21},B_{22}$ .
- Разбиваем C на $C_{11},C_{12},C_{21},C_{22}$ .
- Разбиваем T на $T_{11},T_{12},T_{21},T_{22}$ .
- Распараллеливаем (Fork = вилка, то есть ответвляем процесс):
  - Fork $multiply(C_{11},A_{11},B_{11})$ .
  - Fork $multiply(C_{12},A_{11},B_{12})$ .
  - Fork $multiply(C_{21},A_{21},B_{11})$ .
  - Fork $multiply(C_{22},A_{21},B_{12})$ .
  - Fork $multiply(T_{11},A_{12},B_{21})$ .
  - Fork $multiply(T_{12},A_{12},B_{22})$ .
  - Fork $multiply(T_{21},A_{22},B_{21})$ .
  - Fork $multiply(T_{22},A_{22},B_{22})$ .
- Join (Join = объединение, ждём завершения разветвлённых процессов).
- add(C, T).
- Уничтожаем T.

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

Базовый случай: если $n=1$ , полагаем $c_{11}\leftarrow c_{11}+t_{11}$ (или делаем короткий цикл, возможно, с размоткой).
В противном случае:
- Разбиваем C на C₁₁, C₁₂, C₂₁, C₂₂.
- Разбиваем T на T₁₁, T₁₂, T₂₁, T₂₂.
- Распараллеливаем:
  - Fork $add(C_{11},T_{11})$ .
  - Fork $add(C_{12},T_{12})$ .
  - Fork $add(C_{21},T_{21})$ .
  - Fork $add(C_{22},T_{22})$ .
- Join.

Здесь, fork означает, что вычисления могут осуществляться параллельно к остальной части процедуры, а join означает ожидание завершения всех запущенных в параллельные ветки вычислений. Распараллеливание достигает своей цели лишь передачей указателей.

Этот алгоритм имеет длину критического пути^[en] в $\Theta (\log ^{2}n)$ шагов, что определяет требуемое время для идеальной машины с неограниченным числом процессоров. Поэтому алгоритм имеет максимальное возможное ускорения^[en] $\Theta (n^{3}/\log ^{2}n)$ на любом реальном компьютере. Алгоритм не имеет практического значения ввиду неустранимой цены передачи данных во временную матрицу T и из неё, но более практичный вариант, не использующий временных матриц, достигает ускорения $\Theta (n^{2})$ ^[16].

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

Алгоритмы без обмена данных и распределённые алгоритмы править

В современных архитектурах с иерархической памятью цена загрузки и выгрузки входной матрицы стремится к определяющей роли при обработке. На одной машине это количество данных, переносимых между RAM и кэшем, в то время как для распределённой памяти машин с несколькими узлами это величина, переносимая между узлами. В любом случае это называется полосой пропускания. Наивный алгоритм, где используются три вложенных цикла, использует полосу $\Omega (n^{3})$ .

Алгоритм Кэннона^[en], известный также как алгоритм 2D, — это алгоритм, предотвращающий обмен данными^[en], который превращает каждую входную матрицу в блочную матрицу, элементами которой являются подматрицы размера ${\sqrt {\tfrac {M}{3}}}\times {\sqrt {\tfrac {M}{3}}}$ , где M является размером быстрой памяти ^[17]. Затем используется наивный алгоритм над блоками матриц, вычисляющий произведение подматриц полностью в быстрой памяти. Это сокращает полосу частот канала связи до $O({\tfrac {n^{3}}{\sqrt {M}}})$ , что асимптотически оптимально (для алгоритмов, выполняющих $\Omega (n^{3})$ операций)^[18]^[19].

В распределённых вычислениях с p процессорами, организованными в двухмерную решётку ${\sqrt {p}}\times {\sqrt {p}}$ , каждая из результирующих подматриц может быть назначена процессору и произведение может быть вычислено каждым процессором с передачей $O({\tfrac {n^{2}}{\sqrt {p}}})$ слов, что асимптотически оптимально в предположении, что каждый узел сохраняет минимум $O({\tfrac {n^{2}}{p}})$ элементов^[19]. Это может быть улучшено алгоритмом 3D, который распределяет процессоры в трёхмерную кубическую решётку, путём назначения каждого произведения двух входных подматриц одному процессору. Полученные подматрицы затем генерируются путём работы над каждой строкой^[20]. Этот алгоритм передаёт $O({\tfrac {n^{2}}{p^{2/3}}})$ слов на процессор, что асимптотически оптимально^[19]. Однако, это требует репликации каждого элемента входной матрицы $p^{\tfrac {1}{3}}$ раз, а потому требует в $p^{\tfrac {1}{3}}$ больше памяти, чем нужно для хранения входных данных. Этот алгоритм можно комбинировать с алгоримом Штрассена для дальнейшего уменьшения времени работы^[20]. «2,5D» алгоритмы обеспечивают постоянный обмен между использованием памяти и полосой частот обмена^[21]. В современных системах распределённых вычислений, таких как MapReduce, были разработаны специализированные алгоритмы умножения^[22].

Алгоритмы для ячеистых топологий править

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

n\times n

матриц в ячеистой топологии с перекрёстными связями.

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

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

См. также править

Примечания править

↑ ¹ ² ³ ⁴ Skiena, 2008, с. 45–46, 401–403.
↑ ¹ ² Alman, Williams, 2020.
↑ ¹ ² ³ Кормен, Лейзерсон, Ривест, Штайн, 2005, с. 833-939.
↑ ¹ ² ³ ⁴ ⁵ Amarasinghe, Leiserson, 2010.
↑ Lam, Rothberg, Wolf, 1991.
↑ ¹ ² Prokop, 1999.
↑ Prokop, 1999, с. 13.
↑ Miller, 1975, с. 97–107.
↑ Press, Flannery, Teukolsky, Vetterling, 2007, с. 108.
↑ Iliopoulos, 1989, с. 658–669.
↑ ¹ ² Robinson, 2005.
↑ Raz, 2002, с. 144.
↑ Cohn, Kleinberg, Szegedy, Umans, 2005, с. 379–388.
↑ Cohn, Umans, 2003, с. 438–449.
↑ Alon, Shpilka, Umans, On Sunflowers and Matrix Multiplication Архивная копия от 9 декабря 2016 на Wayback Machine
↑ ¹ ² Randall, 1998, с. 54–57.
↑ Cannon, 1969.
↑ Hong, Kung, 1981, с. 326–333.
↑ ¹ ² ³ Irony, Toledo, Tiskin, 2004, с. 1017–1026.
↑ ¹ ² Agarwal, Balle, Gustavson, Joshi, Palkar, 1995, с. 575–582.
↑ Solomonik, Demmel, 2011, с. 90–109.
↑ Zadeh, Carlsson, 2013.
↑ Bae, Shinn, Takaoka, 2014, с. 2230–2240.
↑ Kak, 1988, с. 383–385.
↑ Kak, Subhash (2014) Efficiency of matrix multiplication on the cross-wired mesh array. https://arxiv.org/abs/1411.3273 Архивная копия от 23 марта 2019 на Wayback Machine
↑ Kak, 1988, с. 347–365.

Литература править

Saman Amarasinghe, Charles Leiserson. 6.172 Performance Engineering of Software Systems, Lecture 8. — Massachusetts Institute of Technology, 2010.
Monica S. Lam, Edward E. Rothberg, Michael E. Wolf. The Cache Performance and Optimizations of Blocked Algorithms // Int'l Conf. on Architectural Support for Programming Languages and Operating Systems (ASPLOS). — 1991.
Webb Miller. Computational complexity and numerical stability // SIAM News. — 1975. — Т. 4, вып. 2. — doi:10.1137/0204009.
William H. Press, Brian P. Flannery, Saul A. Teukolsky, William T. Vetterling. Numerical Recipes: The Art of Scientific Computing. — 3rd. — Cambridge University Press, 2007. — С. 108. — (Numerical Recipes). — ISBN 978-0-521-88068-8.
Lynn Elliot Cannon. A cellular computer to implement the Kalman Filter Algorithm. — Montana State University, 1969. — (Technical report, Ph.D. Thesis).
Ran Raz. On the complexity of matrix product // Proceedings of the Thirty-fourth Annual ACM Symposium on Theory of Computing. — 2002. — ISBN 1581134959. — doi:10.1145/509907.509932.
Henry Cohn, Robert Kleinberg, Balázs Szegedy, Chris Umans. Group-theoretic Algorithms for Matrix Multiplication // Proceedings of the 46th Annual Symposium on Foundations of Computer Science. — Pittsburgh, PA,: IEEE Computer Society, 2005. — С. pp. 379–388.
Henry Cohn, Chris Umans. A Group-theoretic Approach to Fast Matrix Multiplication // Proceedings of the 44th Annual IEEE Symposium on Foundations of Computer Science. — Cambridge, MA,: IEEE Computer Society, 2003. — С. pp. 438–449.
Josh Alman, Virginia Vassilevska Williams. A Refined Laser Method and Faster Matrix Multiplication // 32nd Annual ACM-SIAM Symposium on Discrete Algorithms (SODA 2021). — 2020.
Harald Prokop. Cache-Oblivious Algorithms. — MIT, 1999.
Steven Skiena. Sorting and Searching // The Algorithm Design Manual. — Springer, 2008. — ISBN 978-1-84800-069-8. — doi:10.1007/978-1-84800-070-4_4.
Costas S. Iliopoulos. Worst-case complexity bounds on algorithms for computing the canonical structure of finite abelian groups and the Hermite and Smith normal forms of an integer matrix // SIAM Journal on Computing. — 1989. — Т. 18, вып. 4. — С. 658–669. — doi:10.1137/0218045. Выдержка: «Алгоритм Копперсмита — Виноград непрактичен, поскольку содержит очень большие спрятанные константы в верхней границе числа требуемых умножений.

Sara Robinson. Toward an Optimal Algorithm for Matrix Multiplication // SIAM News. — 2005. — Т. 38, вып. 9.
Keith H. Randall. Cilk: Efficient Multithreaded Computing. — Massachusetts Institute of Technology, 1998. — (Ph.D.).
Hong J. W., Kung H. T. I/O complexity: The red-blue pebble game // STOC '81: Proceedings of the Thirteenth Annual ACM Symposium on Theory of Computing. — 1981.
Dror Irony, Sivan Toledo, Alexander Tiskin. Communication lower bounds for distributed-memory matrix multiplication // J. Parallel Distrib. Comput.. — 2004. — Сентябрь (т. 64, вып. 9). — С. 1017–1026. — doi:10.1016/j.jpdc.2004.03.021.
Agarwal R.C., Balle S. M., Gustavson F. G., Joshi M., Palkar P. A three-dimensional approach to parallel matrix multiplication // IBM J. Res. Dev.. — 1995. — Сентябрь (т. 39, вып. 5). — doi:10.1147/rd.395.0575.
Edgar Solomonik, James Demmel. Communication-optimal parallel 2.5D matrix multiplication and LU factorization algorithms // Proceedings of the 17th International Conference on Parallel Processing. — 2011. — Т. Part II.
Reza Bosagh Zadeh, Gunnar Carlsson. Dimension Independent Matrix Square Using MapReduce. — 2013. — Bibcode: 2013arXiv1304.1467B. — arXiv:1304.1467.
Bae S.E., Shinn T.-W., Takaoka T. A faster parallel algorithm for matrix multiplication on a mesh array // Procedia Computer Science. — 2014. — Т. 29. — doi:10.1016/j.procs.2014.05.208.
Томас Кормен, Чарльз Лейзерсон, Рональд Ривест, Клиффорд Штайн. Алгоритм умножения матриц Штрассена // Алгоритмы: построение и анализ. — Москва, Санко-Петербург, Киев: Вильямс, 2005. — ISBN 5-8459-0857-4.
Kak S. A two-layered mesh array for matrix multiplication // Parallel Computing. — 1988. — Т. 6, вып. 3. — doi:10.1016/0167-8191(88)90078-6.
Kak S. Multilayered array computing // Information Sciences. — 1988. — Т. 45, вып. 3. — doi:10.1016/0020-0255(88)90010-2.

Литература для дальнейшего чтения править

Alfredo Buttari, Julien Langou, Jakub Kurzak, Jack Dongarra. A class of parallel tiled linear algebra algorithms for multicore architectures // Parallel Computing. — 2009. — Т. 35. — С. 38–53. — doi:10.1016/j.parco.2008.10.002. — arXiv:0709.1272.
Kazushige Goto, Robert A. van de Geijn. Anatomy of high-performance matrix multiplication // ACM Transactions on Mathematical Software. — 2008. — Т. 34, вып. 3. — С. 1–25. — doi:10.1145/1356052.1356053.
Field G. Van Zee, Robert A. van de Geijn. BLIS: A Framework for Rapidly Instantiating BLAS Functionality // ACM Transactions on Mathematical Software. — 2015. — Т. 41, вып. 3. — С. 1–33. — doi:10.1145/2764454.
Как оптимизировать GEMM

[_3f1fd1dd340e2ef7-1] ¹ ² ³ ⁴ Skiena, 2008, с. 45–46, 401–403.

[_466e9e163ae2e30e-2] ¹ ² Alman, Williams, 2020.

[_b1191991bfc62747-3] ¹ ² ³ Кормен, Лейзерсон, Ривест, Штайн, 2005, с. 833-939.

[_2de4b0d319ba8cec-4] ¹ ² ³ ⁴ ⁵ Amarasinghe, Leiserson, 2010.

[_1dd448be08e899ec-5] Lam, Rothberg, Wolf, 1991.

[_1ccd63528c82c7be-6] ¹ ² Prokop, 1999.

[_f56bf8d6078b9462-7] Prokop, 1999, с. 13.

[_5b1ea7ac4a274f53-8] Miller, 1975, с. 97–107.

[_e2f287f99646e752-9] Press, Flannery, Teukolsky, Vetterling, 2007, с. 108.

[_edfe97a5ffe11a7c-10] Iliopoulos, 1989, с. 658–669.

[_31d21af877738f22-11] ¹ ² Robinson, 2005.

[_932bf69c4ca08efb-12] Raz, 2002, с. 144.

[_0a9701e6dfd28bad-13] Cohn, Kleinberg, Szegedy, Umans, 2005, с. 379–388.

[_eabd2afbee0b1ea7-14] Cohn, Umans, 2003, с. 438–449.

[15] Alon, Shpilka, Umans, On Sunflowers and Matrix Multiplication Архивная копия от 9 декабря 2016 на Wayback Machine

[_ec66f900021f25ec-16] ¹ ² Randall, 1998, с. 54–57.

[_2b79fbb4d27472e5-17] Cannon, 1969.

[_94f7756b5c80d9e0-18] Hong, Kung, 1981, с. 326–333.

[_175d905b03da02cc-19] ¹ ² ³ Irony, Toledo, Tiskin, 2004, с. 1017–1026.

[_4c4c3422ffe90365-20] ¹ ² Agarwal, Balle, Gustavson, Joshi, Palkar, 1995, с. 575–582.

[_70e1193193cf104c-21] Solomonik, Demmel, 2011, с. 90–109.

[_91de9980478a85a6-22] Zadeh, Carlsson, 2013.

[_94f95a80d1947cbe-23] Bae, Shinn, Takaoka, 2014, с. 2230–2240.

[_ee74851b0c30569d-24] Kak, 1988, с. 383–385.

[Kak-25] Kak, Subhash (2014) Efficiency of matrix multiplication on the cross-wired mesh array. https://arxiv.org/abs/1411.3273 Архивная копия от 23 марта 2019 на Wayback Machine

[_bd8990369b0ccc97-26] Kak, 1988, с. 347–365.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]