Алгоритм Штрассена

Фолькер Штрассен опублікував свій алгоритм в 1969 році. Хоча його алгоритм лише трохи швидший, ніж стандартний алгоритм множення матриць, але він був першим, хто вказав на не оптимальність стандартного підходу. Його стаття надихнула на пошук більш швидких алгоритмів. Зокрема, знайдено більш складний алгоритм Копперсміта-Вінограда в 1987 році.

У лівій колонці представлення множення матриць 2x2. Наївне множення матриць вимагає одне множення для кожного «1» в лівій колонці. Кожен з решти стовпців являє собою єдине з 7 множень в алгоритмі, і сума стовпців дає повне матричне множення зліва.

Нехай A, B — дві квадратні матриці над кільцем R. Ми можемо обчислити матрицю C, як:

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

Якщо матриці A,B, не типу 2ⁿ x 2ⁿ заповнюємо відсутні рядки і стовпці нулями. При цьому ми отримуємо зручні для рекурсивного множення розміри, але втрачаємо ефективність за рахунок додаткових непотрібних множень.

Розділимо матриці A,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_{i, 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_{i, 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}}$

Ми повторюємо рекурсивний процес ділення n раз доти, поки розмір матриць C_i, _j не стане досить малим, далі використовують звичайний метод множення матриць.

Стандартне множення матриць займає приблизно 2N₃ (де N = 2ⁿ) арифметичних операцій (додавання і множення); асимптотична складність O (N³).

Число додавань і множень, необхідних в алгоритмі Штрассена може бути розрахована наступним чином: нехай f(n) число операцій для 2ⁿ × 2ⁿ матриці. Тоді, за рекурсивним застосуванням алгоритму Штрассена, ми бачимо, що f(n) = 7f(n-1) + L4ⁿ, з деякою постійною L, яка залежить від кількості доповнень, виконаних в кожному застосуванні алгоритму. Отже f(n) = (7 + o(1))ⁿ, тобто асимптотична складність для множення матриць розміру N = 2ⁿ використовуючи алгоритм Штрассена є

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

Зменшення кількості арифметичних операцій призводить до частково зменшеної числової стабільності,[1] і алгоритм також вимагає значно більше пам'яті, порівняно з наївним алгоритмом. Обидві початкові матриці, розміри яких повинні бути розширені до наступного ступеня двійки, в результаті чого зберігається до чотирьох разів більше елементів, та сім допоміжних матриць, кожна з яких містить в собі чверть елементів.

Ранг

Ранг є важливим поняттям в асимптотичній складності матричного множення. Ранг ${\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\}}$

Іншими словами, ранг білінійного відображення — це довжина його найкоротшого білінійного обчислення.[2] Існування алгоритму Штрассена показує, що ранг матриці 2×2 здійснює множення не більше ніж сім разів. Щоб переконатися в цьому, розглянемо цей алгоритм (поряд зі стандартним алгоритмом) в таблиці для обчислення матриць.

	Стандартний алгоритм				Алгоритм Штрассена
i	fi(a)	gi(b)	wi		fi(a)	gi(b)	wi
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}}$

Наведений вище опис стверджує, що матриці є квадратними, а розмір є ступенем двійки. Це обмеження дозволяє матриці бути розділеною навпіл рекурсивно, поки межа скалярного множення не буде досягнута. Обмеження спрощує пояснення і аналіз складності; і справді, перетворення матриці, як описано, призведе до збільшення часу обчислень і може легко усунути досить невелику економію часу, отриману за допомогою методу.

Штрассен був першим, хто показав можливість множення матриць більш ефективним способом, ніж стандартний. Після публікації його роботи в 1969, почалися активні пошуки більш швидкого алгоритму. Асимптотично найшвидшим алгоритмом на сьогоднішній день є алгоритм Копперсміта-Винограда, що дозволяє множити матриці за ${\displaystyle {\rm {O}}(n^{2.376})}$ операцій[3], запропонований 1987 року і вдосконалений 2011 року до рівня ${\displaystyle {\rm {O}}(n^{2.3727})}$[3]. Цей алгоритм не має практичного інтересу в оцінці арифметичної складності. Питання про граничну в асимптотиці швидкість множення матриць не вирішене. Існує гіпотеза Штрассена про те, що для достатньо великих n існує алгоритм множення двох матриць розміру ${\displaystyle n\times n}$ за ${\displaystyle {\rm {O}}(n^{2+\varepsilon })}$ операцій, де ${\displaystyle \varepsilon }$ як завгодно мале наперед задане позитивне число. Ця гіпотеза має суто теоретичний інтерес, оскільки розмір матриць для яких ${\displaystyle \varepsilon }$ дійсно дуже великий.

Питання про побудову найбільш швидкого та стійкого практичного алгоритму множення великих матриць також залишається невирішеним.

Існує модифікація алгоритму Штрассена, для якого вимагається 7 множень та 15 додавань (замість18для звичайного алгоритму Штрассена). Матриці A,B і C діляться на блокові підматриці.

Обчислюються проміжні матриціS1 —S8,P1 —P7,T1 —T2

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

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

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

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

${\displaystyle \mathbf {S} _{5}:=(\mathbf {B} _{1,2}-\mathbf {B} _{1,1})}$

${\displaystyle \mathbf {S} _{6}:=(\mathbf {B} _{2,2}-\mathbf {S} _{5})}$

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

${\displaystyle \mathbf {S} _{8}:=(\mathbf {S} _{6}-\mathbf {B} _{2,1})}$

${\displaystyle \mathbf {P} _{1}:=\mathbf {S} _{2}\mathbf {S} _{6}}$

${\displaystyle \mathbf {P} _{2}:=\mathbf {A} _{1,1}\mathbf {B} _{1,1}}$

${\displaystyle \mathbf {P} _{3}:=\mathbf {A} _{1,2}\mathbf {B} _{2,1}}$

${\displaystyle \mathbf {P} _{4}:=\mathbf {S} _{3}\mathbf {S} _{7}}$

${\displaystyle \mathbf {P} _{5}:=\mathbf {S} _{1}\mathbf {S} _{5}}$

${\displaystyle \mathbf {P} _{6}:=\mathbf {S} _{4}\mathbf {B} _{2,2}}$

${\displaystyle \mathbf {P} _{7}:=\mathbf {A} _{2,2}\mathbf {S} _{8}}$

${\displaystyle \mathbf {T} _{1}:=\mathbf {P} _{1}+\mathbf {P} _{2}}$

${\displaystyle \mathbf {T} _{2}:=\mathbf {T} _{1}+\mathbf {P} _{4}}$

Елементи матриці C обчислюються за формулами

${\displaystyle \mathbf {C} _{1,1}:=\mathbf {P} _{2}+\mathbf {P} _{3}}$

${\displaystyle \mathbf {C} _{1,2}:=\mathbf {T} _{1}+\mathbf {P} _{5}+\mathbf {P} _{6}}$

${\displaystyle \mathbf {C} _{2,1}:=\mathbf {T} _{2}-\mathbf {P} _{7}}$

${\displaystyle \mathbf {C} _{2,2}:=\mathbf {T} _{2}+\mathbf {P} _{5}}$

• Weisstein, Eric W. Strassen's Formulas(англ.) на сайті Wolfram MathWorld. (також включені формули для швидкого пошуку зворотньої матриці)
• Tyler J. Earnest, Strassen's Algorithm on the Cell Broadband Engine

• Simple Strassen's algorithms implementation in C (легко для освітніх цілей)

• Simple Strassen's algorithms implementation in Java