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

Алгоритм Штрассена — алгоритм, який в 1969 році розробив Фолькер Штрассен. Призначений для швидкого множення матриць, є узагальненням методу множення Карацуби на матриці. Цей алгоритм дозволяє швидше за стандартний спосіб множити матриці. На практиці це відчутно на великих матрицях, але існує ще більш швидкий алгоритм Копперсміта — Вінограда для множення дуже великих матриць.

На відміну від традиційного алгоритму множення матриць (за формулою $c_{ik}=\sum a_{ij}b_{jk}$), який виконується за час $O(N^{3})=O(N^{\log _{2}8})$, алгоритм Штрассена множить матриці за час $O(N^{\log _{2}7+o(1)})\approx O(N^{2.8074})$, що дає виграш на великих щільних матрицях починаючи, приблизно, з матриць розміру 64 × 64.

Історія

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

Алгоритм

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

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

$\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 на рівні за розміром блочні матриці

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

де

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

тоді

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

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

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

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

При такій конструкції ми не зменшили кількість операцій множення. Нам, як і раніше, потрібно 8 множень для обчислення C_{i, j} матриці, як і в звичайному методі.

Зараз важлива частина. Визначаємо нові матриці

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

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

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

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

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

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

$\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, ось так:

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

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

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

$\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ⁿ використовуючи алгоритм Штрассена є

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

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

Ранг

Ранг є важливим поняттям в асимптотичній складності матричного множення. Ранг $\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×2 здійснює множення не більше ніж сім разів. Щоб переконатися в цьому, розглянемо цей алгоритм (поряд зі стандартним алгоритмом) в таблиці для обчислення матриць.

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

Питання практичної реалізації

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

Подальший розвиток

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$\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 [Архівовано 20 березня 2012 у Wayback Machine.]