Уличный алгоритм

Не путать с алгоритмом Шёнхаге – Штрассена для умножения многочленов.

В линейной алгебре алгоритм Штрассена , названный в честь Фолькера Штрассена , представляет собой алгоритм умножения матриц . Он быстрее, чем стандартный алгоритм умножения матриц для больших матриц, с лучшей асимптотической сложностью , хотя наивный алгоритм часто лучше подходит для меньших матриц. Алгоритм Штрассена медленнее самых быстрых известных алгоритмов для чрезвычайно больших матриц, но такие галактические алгоритмы бесполезны на практике, поскольку они намного медленнее для матриц практического размера. Для небольших матриц существуют еще более быстрые алгоритмы.

Алгоритм Штрассена работает для любого кольца , например плюс/умножение, но не для всех полуколец , например мин-плюс или булевой алгебры , где все еще работает наивный алгоритм, и так называемого комбинаторного умножения матриц .

История

Фолькер Штрассен впервые опубликовал этот алгоритм в 1969 году и тем самым доказал, что $n^{3}$ общий алгоритм умножения матриц не был оптимальным. Публикация алгоритма Штрассена привела к увеличению количества исследований по умножению матриц, которые привели как к асимптотически нижним оценкам, так и к улучшению вычислительных верхних границ.

Алгоритм

В левом столбце визуализируются вычисления, необходимые для определения результата умножения матрицы 2х2 . Наивное умножение матриц требует одного умножения на каждую «1» левого столбца. Каждый из остальных столбцов (M1-M7) представляет одно из 7 умножений в алгоритме Штрассена. Сумма столбцов M1-M7 дает тот же результат, что и умножение полной матрицы слева.

Позволять $A$, $B$ быть двумя квадратными матрицами над кольцом ${\mathcal {R}}$, например матрицы, элементы которых являются целыми или действительными числами. Целью умножения матриц является вычисление произведения матрицы. $C=AB$. Следующее описание алгоритма предполагает, что все эти матрицы имеют размеры, являющиеся степенями двойки (т. е. $A,\,B,\,C\in \operatorname {Matr} _{2^{n}\times 2^{n}}({\mathcal {R}})$), но это лишь концептуально необходимо — если матрицы $A$, $B$ не относятся к типу $2^{n}\times 2^{n}$, «недостающие» строки и столбцы можно заполнить нулями, чтобы получить матрицы с размерами степеней двойки, хотя реальные реализации алгоритма на практике этого не делают.

Разделения алгоритма Штрассена $A$, $B$ и $C$ одинакового размера на блочные матрицы

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

с $A_{ij},B_{ij},C_{ij}\in \operatorname {Mat} _{2^{n-1}\times 2^{n-1}}({\mathcal {R}})$. Наивный алгоритм будет таким:

${\begin{bmatrix}C_{11}&C_{12}\\C_{21}&C_{22}\end{bmatrix}}={\begin{bmatrix}A_{11}{\color {red}\times }B_{11}+A_{12}{\color {red}\times }B_{21}\quad &A_{11}{\color {red}\times }B_{12}+A_{12}{\color {red}\times }B_{22}\\A_{21}{\color {red}\times }B_{11}+A_{22}{\color {red}\times }B_{21}\quad &A_{21}{\color {red}\times }B_{12}+A_{22}{\color {red}\times }B_{22}\end{bmatrix}}.$

Данная конструкция не уменьшает количество умножений: для вычисления $C_{ij}$ матриц, то же количество умножений, необходимое при использовании стандартного умножения матриц.

Вместо этого алгоритм Штрассена определяет новые значения:

${\begin{aligned}M_{1}&=(A_{11}+A_{22}){\color {red}\times }(B_{11}+B_{22});\\M_{2}&=(A_{21}+A_{22}){\color {red}\times }B_{11};\\M_{3}&=A_{11}{\color {red}\times }(B_{12}-B_{22});\\M_{4}&=A_{22}{\color {red}\times }(B_{21}-B_{11});\\M_{5}&=(A_{11}+A_{12}){\color {red}\times }B_{22};\\M_{6}&=(A_{21}-A_{11}){\color {red}\times }(B_{11}+B_{12});\\M_{7}&=(A_{12}-A_{22}){\color {red}\times }(B_{21}+B_{22}),\\\end{aligned}}$

используя всего 7 умножений (по одному на каждое $M_{k}$) вместо 8. Теперь мы можем выразить $C_{ij}$ с точки зрения $M_{k}$:

${\begin{bmatrix}C_{11}&C_{12}\\C_{21}&C_{22}\end{bmatrix}}={\begin{bmatrix}M_{1}+M_{4}-M_{5}+M_{7}\quad &M_{3}+M_{5}\\M_{2}+M_{4}\quad &M_{1}-M_{2}+M_{3}+M_{6}\end{bmatrix}}.$

Мы рекурсивно повторяем этот процесс деления до тех пор, пока подматрицы не вырождаются в числа (элементы кольца ${\mathcal {R}}$). Если, как упоминалось выше, исходная матрица имела размер, не являющийся степенью 2, то результирующее произведение будет иметь нулевые строки и столбцы, как и $A$ и $B$, и затем на этом этапе они будут удалены, чтобы получить (меньшую) матрицу $C$ мы очень хотели.

Практическая реализация алгоритма Штрассена переключается на стандартные методы умножения матриц для достаточно маленьких подматриц, для которых эти алгоритмы более эффективны. Конкретная точка пересечения, для которой алгоритм Штрассена более эффективен, зависит от конкретной реализации и аппаратного обеспечения. Ранее авторы подсчитали, что алгоритм Штрассена работает быстрее для матриц шириной от 32 до 128 для оптимизированных реализаций. Однако было замечено, что эта точка пересечения в последние годы увеличивается, и исследование 2010 года показало, что даже один шаг алгоритма Штрассена часто не приносит пользы для современных архитектур по сравнению с высокооптимизированным традиционным умножением, пока размеры матрицы не превысят 1000 или более, и даже для размеров матриц в несколько тысяч выигрыш обычно в лучшем случае незначителен (около 10% или меньше). Более недавнее исследование (2016 г.) выявило преимущества для матриц размером всего 512 и около 20%.

Форма Винограда

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

${\begin{bmatrix}a&b\\c&d\end{bmatrix}}{\begin{bmatrix}A&C\\B&D\end{bmatrix}}={\begin{bmatrix}t+b{\color {red}\times }B&w+v+(a+b-c-d){\color {red}\times }D\\w+u+d{\color {red}\times }(B+C-A-D)&w+u+v\end{bmatrix}}$

где $t=a{\color {red}\times }A,\;u=(c-a){\color {red}\times }(C-D),\;v=(c+d){\color {red}\times }(C-A),\;w=t+(c+d-a){\color {red}\times }(A+D-C)$.

Это уменьшает количество сложений и вычитаний матриц с 18 до 15. Количество умножений матриц по-прежнему равно 7, а асимптотическая сложность та же.

Асимптотическая сложность

Схема приведенного выше алгоритма показала, что можно обойтись всего 7 вместо традиционных 8 матричных умножений для подблоков матрицы. С другой стороны, приходится делать сложение и вычитание блоков, хотя по общей сложности это не имеет значения: Сложение матриц размера $N/2$ требуется только $(N/2)^{2}$ операции, тогда как умножение существенно дороже (традиционно $2(N/2)^{3}$ операции сложения или умножения).

Тогда вопрос заключается в том, сколько именно операций нужно для алгоритмов Штрассена и как это соотносится со стандартным умножением матриц, которое занимает примерно $2N^{3}$ (где $N=2^{n}$) арифметические операции, т.е. асимптотическая сложность $\Theta (N^{3})$.

Количество сложений и умножений, необходимых в алгоритме Штрассена, можно рассчитать следующим образом: пусть $f(n)$ быть числом операций за $2^{n}\times 2^{n}$ матрица. Затем, рекурсивно применив алгоритм Штрассена, мы видим, что $f(n)=7f(n-1)+l4^{n}$, для некоторой константы $l$ это зависит от количества сложений, выполняемых при каждом применении алгоритма. Следовательно $f(n)=(7+o(1))^{n}$, т. е. асимптотическая сложность умножения матриц размера $N=2^{n}$ использование алгоритма Штрассена $O([7+o(1)]^{n})=O(N^{\log _{2}7+o(1)})\approx O(N^{2.8074})$. Однако за уменьшение количества арифметических операций приходится несколько снижать числовую стабильность . и алгоритм также требует значительно больше памяти по сравнению с наивным алгоритмом. Обе исходные матрицы должны иметь размеры, расширенные до следующей степени 2, что приводит к хранению в четыре раза большего количества элементов, а каждая из семи вспомогательных матриц содержит четверть элементов расширенных.

Алгоритм Штрассена следует сравнить с «наивным» способом умножения матриц, который потребует 8 вместо 7 умножений подблоков. Тогда это привело бы к сложности, которую можно ожидать от стандартного подхода: $O(8^{n})=O(N^{\log _{2}8})=O(N^{3})$. Сравнение этих двух алгоритмов показывает, что асимптотически алгоритм Штрассена быстрее: существует размер $N_{\text{threshold}}$ так что матрицы большего размера более эффективно умножаются с помощью алгоритма Штрассена, чем «традиционным» способом. Однако из асимптотического утверждения не следует, что алгоритм Штрассена всегда быстрее даже для маленьких матриц, и на практике это на самом деле не так: для маленьких матриц стоимость дополнительных добавлений матричных блоков перевешивает экономию на количестве умножения. Существуют и другие факторы, не учтенные в приведенном выше анализе, например, разница в стоимости современного оборудования между загрузкой данных из памяти в процессоры и стоимостью фактического выполнения операций с этими данными. Вследствие такого рода соображений алгоритм Штрассена обычно используется только с «большими» матрицами. Этот вид эффекта еще более выражен при использовании альтернативных алгоритмов, таких как алгоритм Копперсмита и Винограда : хотя асимптотически даже быстрее, точка пересечения $N_{\text{threshold}}$ настолько велик, что этот алгоритм обычно не используется для матриц, встречающихся на практике.

Ранг или билинейная сложность

Билинейная сложность или ранг билинейного отображения является важным понятием асимптотической сложности умножения матриц. Ранг билинейного отображения $\phi :\mathbf {A} \times \mathbf {B} \rightarrow \mathbf {C}$ над полем F определяется как (что-то вроде злоупотребления обозначениями )

$R(\phi /\mathbf {F} )=\min \left\{r\left|\exists f_{i}\in \mathbf {A} ^{*},g_{i}\in \mathbf {B} ^{*},w_{i}\in \mathbf {C} ,\forall \mathbf {a} \in \mathbf {A} ,\mathbf {b} \in \mathbf {B} ,\phi (\mathbf {a} ,\mathbf {b} )=\sum _{i=1}^{r}f_{i}(\mathbf {a} )g_{i}(\mathbf {b} )w_{i}\right.\right\}$

Другими словами, ранг билинейного отображения — это длина его кратчайшего билинейного вычисления. Существование алгоритма Штрассена показывает, что ранг $2\times 2$ умножение матрицы не более семи. Чтобы убедиться в этом, давайте представим этот алгоритм (наряду со стандартным алгоритмом) как такое билинейное вычисление. В случае матриц двойственные пространства A * и B * состоят из отображений в поле F , индуцированных скалярным двойным скалярным произведением (т.е. в данном случае суммой всех элементов произведения Адамара ).

	Стандартный алгоритм				Уличный алгоритм
$i$	$f_{i}(\mathbf {a} )$	$g_{i}(\mathbf {b} )$	$w_{i}$		$f_{i}(\mathbf {a} )$	$g_{i}(\mathbf {b} )$	$w_{i}$
1	${\begin{bmatrix}1&0\\0&0\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}1&0\\0&0\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}1&0\\0&0\end{bmatrix}}$		${\begin{bmatrix}1&0\\0&1\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}1&0\\0&1\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}1&0\\0&1\end{bmatrix}}$
2	${\begin{bmatrix}0&1\\0&0\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}0&0\\1&0\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}1&0\\0&0\end{bmatrix}}$		${\begin{bmatrix}0&0\\1&1\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}1&0\\0&0\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}0&0\\1&-1\end{bmatrix}}$
3	${\begin{bmatrix}1&0\\0&0\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}0&1\\0&0\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}0&1\\0&0\end{bmatrix}}$		${\begin{bmatrix}1&0\\0&0\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}0&1\\0&-1\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}0&1\\0&1\end{bmatrix}}$
4	${\begin{bmatrix}0&1\\0&0\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}0&0\\0&1\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}0&1\\0&0\end{bmatrix}}$		${\begin{bmatrix}0&0\\0&1\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}-1&0\\1&0\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}1&0\\1&0\end{bmatrix}}$
5	${\begin{bmatrix}0&0\\1&0\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}1&0\\0&0\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}0&0\\1&0\end{bmatrix}}$		${\begin{bmatrix}1&1\\0&0\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}0&0\\0&1\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}-1&1\\0&0\end{bmatrix}}$
6	${\begin{bmatrix}0&0\\0&1\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}0&0\\1&0\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}0&0\\1&0\end{bmatrix}}$		${\begin{bmatrix}-1&0\\1&0\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}1&1\\0&0\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}0&0\\0&1\end{bmatrix}}$
7	${\begin{bmatrix}0&0\\1&0\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}0&1\\0&0\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}0&0\\0&1\end{bmatrix}}$		${\begin{bmatrix}0&1\\0&-1\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}0&0\\1&1\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}1&0\\0&0\end{bmatrix}}$
8	${\begin{bmatrix}0&0\\0&1\end{bmatrix}}:\mathbf {a}$	${\begin{bmatrix}0&0\\0&1\end{bmatrix}}:\mathbf {b}$	${\begin{bmatrix}0&0\\0&1\end{bmatrix}}$
	$\mathbf {a} \mathbf {b} =\sum _{i=1}^{8}f_{i}(\mathbf {a} )g_{i}(\mathbf {b} )w_{i}$				$\mathbf {a} \mathbf {b} =\sum _{i=1}^{7}f_{i}(\mathbf {a} )g_{i}(\mathbf {b} )w_{i}$

Можно показать, что общее количество элементарных умножений $L$ необходимое для умножения матриц тесно асимптотически связано с рангом $R$, то есть $L=\Theta (R)$или, более конкретно, поскольку константы известны, $R/2\leq L\leq R$. Одним из полезных свойств ранга является то, что он субмультипликативен для тензорных произведений , и это позволяет показать, что $2^{n}\times 2^{n}\times 2^{n}$ умножение матриц может быть выполнено не более чем за $7n$ элементарное умножение любого $n$. (Этот $n$-кратное тензорное произведение $2\times 2\times 2$ карта умножения матрицы сама на себя — $n$-я тензорная степень — реализуется на рекурсивном этапе показанного алгоритма.)

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

Алгоритм Штрассена не учитывает кэш . Анализ алгоритма поведения кэша показал, что он несет

$\Theta \left(1+{\frac {n^{2}}{b}}+{\frac {n^{\log _{2}7}}{b{\sqrt {M}}}}\right)$

кеш промахивается во время своего выполнения, предполагая, что кеш идеализированного размера $M$ (т.е. с $M/b$ линии длины $b$).

Рекомендации по реализации

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

Хорошая реализация будет наблюдать следующее:

• Нет необходимости и нежелательно использовать алгоритм Штрассена вплоть до предела скаляров. По сравнению с обычным умножением матриц этот алгоритм добавляет значительные $O(n^{2})$ нагрузка на сложение/вычитание; поэтому ниже определенного размера лучше использовать обычное умножение. Так, например, $1600\times 1600$ не нужно дополнять $2048\times 2048$, поскольку его можно разделить на $25\times 25$ Затем на этом уровне можно использовать матрицы и обычное умножение.
• Этот метод действительно можно применить к квадратным матрицам любой размерности. Если размер четный, они делятся пополам, как описано. Если размерность нечетная, сначала применяется заполнение нулями одной строки и одного столбца. Такое дополнение можно применять «на лету» и лениво, а лишние строки и столбцы отбрасываются по мере формирования результата. Например, предположим, что матрицы $199\times 199$. Их можно разделить так, чтобы верхняя левая часть была $100\times 100$ и правый нижний угол $99\times 99$. Везде, где этого требуют операции, размеры $99$ дополняются нулями до $100$ первый. Обратите внимание, например, что произведение $M_{2}$ используется только в нижней строке вывода, поэтому требуется только $99$ ряды высокие; и, следовательно, левый фактор $A_{21}+A_{22}$ использовать для его создания нужно только $99$ ряды высокие; соответственно, нет необходимости дополнять эту сумму до $100$ ряды; нужно только проложить $A_{22}$ к $100$ столбцы для соответствия $A_{21}$.

Более того, нет необходимости, чтобы матрицы были квадратными. Неквадратные матрицы можно разделить пополам, используя те же методы, в результате чего неквадратные матрицы будут меньшего размера. Если матрицы достаточно неквадратны, имеет смысл свести первоначальную операцию к более квадратным произведениям, используя простые методы, которые по сути являются $O(n^{2})$, например:

• Продукт размером $[2N\times N]\ast [N\times 10N]$ можно сделать как 20 отдельных $[N\times N]\ast [N\times N]$ операции, организованные для формирования результата;
• Продукт размером $[N\times 10N]\ast [10N\times N]$ можно сделать как 10 отдельных $[N\times N]\ast [N\times N]$ операции, суммированные для формирования результата.

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

На практике алгоритм Штрассена может быть реализован для достижения большей производительности, чем обычное умножение, даже для матриц такого маленького размера, как $500\times 500$, для матриц, которые совсем не являются квадратными и не требуют рабочего пространства за пределами буферов, которые уже необходимы для высокопроизводительного обычного умножения.

См. также

• Вычислительная сложность математических операций
• Элиминация Гаусса – Джордана
• Вычислительная сложность умножения матриц
• Кривая Z-порядка
• Алгоритм Карацубы для умножения n -значных целых чисел в $O(n^{\log _{2}3})$ вместо того, чтобы в $O(n^{2})$ время
  • Похожий алгоритм комплексного умножения умножает два комплексных числа, используя 3 действительных умножения вместо 4.
• Алгоритм Тума-Кука , более быстрое обобщение алгоритма Карацубы, которое позволяет рекурсивное разложение по принципу «разделяй и властвуй» на более чем два блока одновременно.

Ссылки

• Томас Х. Кормен , Чарльз Э. Лейзерсон , Рональд Л. Ривест и Клиффорд Стейн . Введение в алгоритмы , второе издание. MIT Press и McGraw-Hill, 2001. ISBN   0-262-03293-7 . Глава 28: Раздел 28.2: Алгоритм Штрассена для умножения матриц, стр. 735–741.
• Кнут, Дональд (1997). Искусство компьютерного программирования, получисловые алгоритмы . Том. II (3-е изд.). Аддисон-Уэсли. ISBN  0-201-89684-2 .

Внешние ссылки

• Вайсштейн, Эрик В. «Формулы Штрассена» . Математический мир . (также включает формулы для быстрого обращения матрицы )
• Тайлер Дж. Эрнест, Алгоритм Штрассена для механизма сотовой широкополосной связи