Алгоритм Штрассена - Strassen algorithm - Wikipedia

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

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

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

История

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

Алгоритм

Левый столбец представляет 2x2 матричное умножение. Наивное матричное умножение требует одного умножения на каждую единицу левого столбца. Каждый из других столбцов представляет собой одно из 7 умножений в алгоритме, а сумма столбцов дает полное матричное умножение слева.

Позволять А, B быть двумя квадратные матрицы через звенеть р. Мы хотим вычислить матричное произведение C в качестве

${\displaystyle \mathbf {C} =\mathbf {A} \mathbf {B} \qquad \mathbf {A} ,\mathbf {B} ,\mathbf {C} \in R^{2^{n}\times 2^{n}}}$

Если матрицы А, B не относятся к типу 2^п × 2^п заполняем недостающие строки и столбцы нулями.

Мы разделяем А, B и C в равные по размеру блочные матрицы

${\displaystyle \mathbf {A} ={\begin{bmatrix}\mathbf {A} _{1,1}&\mathbf {A} _{1,2}\\\mathbf {A} _{2,1}&\mathbf {A} _{2,2}\end{bmatrix}}{\mbox{ , }}\mathbf {B} ={\begin{bmatrix}\mathbf {B} _{1,1}&\mathbf {B} _{1,2}\\\mathbf {B} _{2,1}&\mathbf {B} _{2,2}\end{bmatrix}}{\mbox{ , }}\mathbf {C} ={\begin{bmatrix}\mathbf {C} _{1,1}&\mathbf {C} _{1,2}\\\mathbf {C} _{2,1}&\mathbf {C} _{2,2}\end{bmatrix}}}$

с

${\displaystyle \mathbf {A} _{i,j},\mathbf {B} _{i,j},\mathbf {C} _{i,j}\in R^{2^{n-1}\times 2^{n-1}}}$.

Наивный алгоритм был бы таким:

${\displaystyle \mathbf {C} _{1,1}=\mathbf {A} _{1,1}\mathbf {B} _{1,1}+\mathbf {A} _{1,2}\mathbf {B} _{2,1}}$

${\displaystyle \mathbf {C} _{1,2}=\mathbf {A} _{1,1}\mathbf {B} _{1,2}+\mathbf {A} _{1,2}\mathbf {B} _{2,2}}$

${\displaystyle \mathbf {C} _{2,1}=\mathbf {A} _{2,1}\mathbf {B} _{1,1}+\mathbf {A} _{2,2}\mathbf {B} _{2,1}}$

${\displaystyle \mathbf {C} _{2,2}=\mathbf {A} _{2,1}\mathbf {B} _{1,2}+\mathbf {A} _{2,2}\mathbf {B} _{2,2}}$

С помощью этой конструкции мы не уменьшили количество умножений. Нам все еще нужно 8 умножений, чтобы вычислить C_{я, j} матрицы, то же количество умножений, которое нам нужно при использовании стандартного умножения матриц.

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

${\displaystyle \mathbf {M} _{1}:=(\mathbf {A} _{1,1}+\mathbf {A} _{2,2})(\mathbf {B} _{1,1}+\mathbf {B} _{2,2})}$

${\displaystyle \mathbf {M} _{2}:=(\mathbf {A} _{2,1}+\mathbf {A} _{2,2})\mathbf {B} _{1,1}}$

${\displaystyle \mathbf {M} _{3}:=\mathbf {A} _{1,1}(\mathbf {B} _{1,2}-\mathbf {B} _{2,2})}$

${\displaystyle \mathbf {M} _{4}:=\mathbf {A} _{2,2}(\mathbf {B} _{2,1}-\mathbf {B} _{1,1})}$

${\displaystyle \mathbf {M} _{5}:=(\mathbf {A} _{1,1}+\mathbf {A} _{1,2})\mathbf {B} _{2,2}}$

${\displaystyle \mathbf {M} _{6}:=(\mathbf {A} _{2,1}-\mathbf {A} _{1,1})(\mathbf {B} _{1,1}+\mathbf {B} _{1,2})}$

${\displaystyle \mathbf {M} _{7}:=(\mathbf {A} _{1,2}-\mathbf {A} _{2,2})(\mathbf {B} _{2,1}+\mathbf {B} _{2,2})}$

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

${\displaystyle \mathbf {C} _{1,1}=\mathbf {M} _{1}+\mathbf {M} _{4}-\mathbf {M} _{5}+\mathbf {M} _{7}}$

${\displaystyle \mathbf {C} _{1,2}=\mathbf {M} _{3}+\mathbf {M} _{5}}$

${\displaystyle \mathbf {C} _{2,1}=\mathbf {M} _{2}+\mathbf {M} _{4}}$

${\displaystyle \mathbf {C} _{2,2}=\mathbf {M} _{1}-\mathbf {M} _{2}+\mathbf {M} _{3}+\mathbf {M} _{6}}$

Мы повторяем этот процесс деления п раз (рекурсивно), пока подматрицы вырождаются в числа (элементы кольца р). Полученный продукт будет дополнен нулями, как и А и B, и следует удалить соответствующие строки и столбцы.

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

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

Стандартное матричное умножение занимает примерно 2N³ (куда N = 2^п) арифметические операции (сложение и умножение); асимптотическая сложность Θ (N³).

Количество сложений и умножений, требуемых в алгоритме Штрассена, можно рассчитать следующим образом: пусть ж(п) быть количеством операций для 2^п × 2^п матрица. Затем, рекурсивно применяя алгоритм Штрассена, мы видим, что ж(п) = 7ж(п−1) + ℓ4^п, для некоторой постоянной ℓ это зависит от количества добавлений, выполняемых при каждом применении алгоритма. Следовательно ж(п) = (7 + o (1))^п, т.е. асимптотическая сложность умножения матриц размера N = 2^п с использованием алгоритма Штрассена

${\displaystyle O([7+o(1)]^{n})=O(N^{\log _{2}7+o(1)})\approx O(N^{2.8074})}$.

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

«Наивный» способ выполнения матричного умножения потребовал бы 8 вместо 7 умножений подблоков. Тогда это вызовет сложность, которую можно ожидать от стандартного подхода:

${\displaystyle O(8^{log_{2}n})=O(N^{\log _{2}8})=O(N^{3})}$.

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

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

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

Другими словами, ранг билинейного отображения - это длина его кратчайшего билинейного вычисления.^[5] Существование алгоритма Штрассена показывает, что ранг умножения матриц 2 × 2 не превышает семи. Чтобы убедиться в этом, представим этот алгоритм (наряду со стандартным алгоритмом) как такое билинейное вычисление. В случае матриц двойные пространства А* и B* состоят из карт в поле F индуцированный скаляром произведение с двумя точками, (т.е. в данном случае сумма всех записей Произведение Адамара.)

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

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

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

Алгоритм Штрассена кеш не обращает внимания. Анализ ее тайник алгоритм поведения показал, что

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

промахи кеша во время его выполнения, предполагая идеализированный размер кеша M (т.е. с ${\displaystyle {\frac {M}{b}}}$ линии длины б).^[6]:13

Соображения по реализации

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

В хорошей реализации будет соблюдаться следующее:

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

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

• Изделие размером [2N Икс N] * [N х 10N] можно сделать как 20 отдельных [N Икс N] * [N Икс N] операции, упорядоченные по формированию результата;
• Изделие размером [N х 10N] * [10N Икс N] можно сделать как 10 отдельных [N Икс N] * [N Икс N] операций, суммированных для формирования результата.

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

На практике алгоритм Штрассена может быть реализован для достижения лучшей производительности, чем обычное умножение, даже для небольших матриц, для матриц, которые совсем не квадратные, и без необходимости рабочего пространства за пределами буферов, которые уже необходимы для высокопроизводительного обычного умножения.^[3]

Смотрите также

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

Рекомендации

1. Скиена, Стивен С. (1998), "§8.2.3 Умножение матриц", Руководство по разработке алгоритмов, Берлин, Нью-Йорк: Springer-Verlag, ISBN  978-0-387-94860-7.
2. Д'Альберто, Паоло; Николау, Александру (2005). Использование рекурсии для повышения производительности ATLAS (PDF). Шестой международный симпозиум. по высокопроизводительным вычислениям.
3. Хуанг, Цзяньюй; Смит, Тайлер; Генри, Грег; ван де Гейн, Роберт (2016). Перезагрузка алгоритма Штрассена. Международная конференция по высокопроизводительным вычислениям, сетям, хранению данных и анализу (SC'16).
4. Уэбб, Миллер (1975). «Вычислительная сложность и численная устойчивость». SIAM J. Comput.: 97–107.
5. Бюргиссер, Клаузен и Шокроллахи. Алгебраическая теория сложности. Springer-Verlag 1997.
6. Frigo, M .; Лейзерсон, К.Э.; Прокоп, Х.; Рамачандран, С. (1999). Алгоритмы без кеширования (PDF). Proc. IEEE Symp. по основам компьютерных наук (FOCS). С. 285–297.
7. Хайэм, Николас Дж. (1990). «Использование быстрого матричного умножения в BLAS уровня 3» (PDF). Транзакции ACM на математическом ПО. 16 (4): 352–368. Дои:10.1145/98267.98290. HDL:1813/6900. S2CID  5715053.

• Штрассен, Фолькер, Исключение Гаусса не оптимально, Нумер. Математика. 13, стр. 354-356, 1969 г.
• Томас Х. Кормен, Чарльз Э. Лейзерсон, Рональд Л. Ривест, и Клиффорд Штайн. Введение в алгоритмы, Второе издание. MIT Press и McGraw-Hill, 2001. ISBN  0-262-03293-7. Глава 28: Раздел 28.2: Алгоритм Штрассена для умножения матриц, стр. 735–741.

внешняя ссылка

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