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

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

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

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

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

По определению умножения матриц для ${\displaystyle n\times m}$ матрицы A и по ${\displaystyle m\times p}$ матрицы B произведением ${\displaystyle C=AB}$ является ${\displaystyle n\times p}$ матрица, состоящая из элементов

${\displaystyle 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:
    • Положим ${\displaystyle sum=0}$
    • Для k от 1 до m:
      • Положим ${\displaystyle sum\leftarrow sum+A_{ik}\times B_{kj}}$
    • Положим ${\displaystyle C_{ij}\leftarrow sum}$
• Возвращаем C

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

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

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

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

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

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

• Вход: матрицы A и B
• Пусть C будет новой матрицей нужного размера
• Выбираем часть размера ${\displaystyle T=\Theta ({\sqrt {M}})}$
• Для I от 1 до n для части T:
  • Для J от 1 до p для части T:
    • Для K от 1 до m для части T:
      • Умножаем ${\displaystyle A_{I:I+T,K:K+T}}$ и ${\displaystyle B_{K:K+T,J:J+T}}$, помещая результат в ${\displaystyle C_{I:I+T,J:J+T}}$, то есть:
      • Для i от I до ${\displaystyle \min(I+T,n)}$:
        • Для j от J до ${\displaystyle \min(J+T,p)}$:
          • Положим sum = 0
          • Для k от K до ${\displaystyle \min(K+T,m)}$:
            • Положим ${\displaystyle sum\leftarrow sum+A_{ik}\times B_{kj}}$
          • Положим ${\displaystyle C_{ij}\leftarrow C_{ij}+sum}$
• Возвращаем C

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

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

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

${\displaystyle 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}}}$,

что работает для всех квадратных матриц с размерностями, равными степеням двойки, то есть ${\displaystyle 2^{n}\times 2^{n}}$ для некоторого n. Произведение матриц тогда равно

${\displaystyle {\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}}}$

что составляет восемь умножений пар подматриц с последующим шагом сложения. Алгоритм разделяй-и-властвуй вычисляет элементы рекурсивно с помощью скалярного произведения ${\displaystyle c_{11}=a_{11}b_{11}}$ как в базовом случае.

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

${\displaystyle T(1)=\Theta (1)}$;

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

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

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

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

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

• Если ${\displaystyle \max(n,m,p)=n}$, разбиваем матрицу A горизонтально:

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

• В противном случае, если ${\displaystyle \max(n,m,p)=p}$, разбиваем матрицу B вертикально:

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

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

${\displaystyle 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лгоритма с разбиением на блоки^[англ.], но, в отличие от этого алгоритма, рекурсивный алгоритм нечувствителен к кешированию^{[англ.][6]} — не нужно никакого настоечного параметра для получения оптимального поведения кэша и он работает хорошо в многозадачном окружении, где размер кэша меняется динамически, поскольку другие процессы тоже используют кэш^[4]. (Простой итеративный алгоритм нечувствителен к кэшированию также, но на практике много медленнее, если матрицы не подогнаны под алгоритм.)

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

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

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

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

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

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

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

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

Алгоритм Фрейвалдса — это простой алгоритм Монте-Карло, который для заданных матриц ${\displaystyle A,B}$ и C проверяет, что ${\displaystyle AB=C}$, за время ${\displaystyle \Theta (n^{2})}$.

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

Параллельность с разделением памяти

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

${\displaystyle {\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):

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

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

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

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

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

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

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

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

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

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

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

Умножение матриц, выполненное в 2n-1 шагов для двух ${\displaystyle n\times n}$ матриц в ячеистой топологии с перекрёстными связями.

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

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

См. также

• Вычислительная сложность математических операций
• Алгоритм Валианта
• Задача о порядке перемножения матриц

Примечания

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