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

Алгоритм Штрассена предназначен для быстрого умножения матриц. Он был разработан Фолькером Штрассеном в 1969 году и является обобщением метода умножения Карацубы на матрицы.

В отличие от традиционного алгоритма умножения матриц (по формуле $c_{ij}=\sum a_{ik}b_{kj}$), работающего за время $\Theta (n^{\log _{2}8})=\Theta (n^{3})$, алгоритм Штрассена умножает матрицы за время $\Theta (n^{\log _{2}7})=O(n^{2.81})$, что даёт выигрыш на больших плотных матрицах.

Несмотря на то, что алгоритм Штрассена является асимптотически не самым быстрым из существующих алгоритмов быстрого умножения матриц, он проще программируется и эффективнее при умножении матриц относительно малого размера, поэтому именно он чаще используется на практике.

Описание алгоритма

Для умножения (2x2)-матриц используются произведения сумм и разностей элементов. На иллюстрации каждое такое произведение взято в отдельную цветную рамку, а способ их компоновки в финальную матрицу обозначен исходящими лучами.

Если добавить к матрицам $A$ и $B$ одинаковые нулевые строки и столбцы, их произведение станет равно матрице $AB$ с теми же добавленными строками и столбцами. Поэтому можно рассматривать только матрицы размера $n=2^{k},\ k\in {\mathbb {N} }$, а другие случаи сводить к этому добавлением нулей, отчего $n$ может увеличиться лишь вдвое.

Пусть $A,B$ – матрицы размера $2^{k}\times 2^{k}$. Их можно представить как блочные матрицы размера $(2\times 2)$ из $(2^{k-1}\times 2^{k-1})$–матриц:

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

По принципу блочного умножения, матрица $AB$ выражается через их произведение

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

где в правой части происходит восемь умножений матриц размера $2^{k-1}\times 2^{k-1}$. Поскольку матрицы образуют кольцо, то для вычисления правой части годится любой алгоритм умножения $(2\times 2)$-матриц, использующий лишь сложения, вычитания и умножения. Штрассен предложил такой алгоритм с семью умножениями:

${\begin{aligned}D&=(A_{11}+A_{22})(B_{11}+B_{22});\\D_{1}&=(A_{12}-A_{22})(B_{21}+B_{22});\\D_{2}&=(A_{21}-A_{11})(B_{11}+B_{12});\\H_{1}&=(A_{11}+A_{12})B_{22};\\H_{2}&=(A_{21}+A_{22})B_{11};\\V_{1}&=A_{22}(B_{21}-B_{11});\\V_{2}&=A_{11}(B_{12}-B_{22});\\\end{aligned}}$

${\begin{aligned}AB&={\begin{pmatrix}D&0\\0&D\end{pmatrix}}+{\begin{pmatrix}D_{1}&0\\0&D_{2}\end{pmatrix}}+{\begin{pmatrix}-H_{1}&H_{1}\\H_{2}&-H_{2}\end{pmatrix}}+{\begin{pmatrix}V_{1}&V_{2}\\V_{1}&V_{2}\end{pmatrix}}\\&={\begin{pmatrix}D+D_{1}+V_{1}-H_{1}&V_{2}+H_{1}\\V_{1}+H_{2}&D+D_{2}+V_{2}-H_{2}\end{pmatrix}}\ .\end{aligned}}$

Каждое умножение можно совершать рекурсивно по той же процедуре, а сложение – тривиально, складывая $(2^{k-1})^{2}$ элементов. Тогда время работы алгоритма $T(n)$ оценивается через рекуррентное соотношение:

$T(n)=7T(n/2)+O(n^{2})=O(n^{\log _{2}7})\ .$

Пример реализации

Ниже приведён пример реализации алгоритма на языке Python с использованием библиотеки NumPy для быстрого взятия подматриц. Основная функция – strassen_mul. Предполагается, что все матрицы квадратны, представлены типом numpy.array и их размер является степенью 2.

При небольших размерах матрицы прямое умножение оказывается быстрее алгоритма Штрассена из-за большого числа сложений в последнем. Граница таких размеров зависит от соотношения времени сложения и умножения элементов и поэтому может варьироваться в зависимости от аппаратной среды. В коде за её назначение отвечает константа TRIVIAL_MULTIPLICATION_BOUND.

from itertools import product
import numpy as np

def split_to_2x2_blocks(matrix):
	return list(map(
		lambda row: np.hsplit(row, 2),
		np.vsplit(matrix, 2)
	))

def strassen_mul_2x2(lb, rb):
	d = strassen_mul(lb[0][0] + lb[1][1], rb[0][0] + rb[1][1])
	d_1 = strassen_mul(lb[0][1] - lb[1][1], rb[1][0] + rb[1][1])
	d_2 = strassen_mul(lb[1][0] - lb[0][0], rb[0][0] + rb[0][1])

	left = strassen_mul(lb[1][1], rb[1][0] - rb[0][0])
	right = strassen_mul(lb[0][0], rb[0][1] - rb[1][1])
	top = strassen_mul(lb[0][0] + lb[0][1], rb[1][1])
	bottom = strassen_mul(lb[1][0] + lb[1][1], rb[0][0])

	return [[d + d_1 + left - top, right + top],
			[left + bottom, d + d_2 + right - bottom]]

def trivial_mul(left, right):
	height, mid_size = left.shape
	mid_size, width = right.shape

	result = np.zeros((height, width))
	for row, col, mid in product(*map(range, [height, width, mid_size])):
		result[row][col] += left[row][mid] * right[mid][col]

	return result

TRIVIAL_MULTIPLICATION_BOUND = 8

def strassen_mul(left, right):
	assert(left.shape == right.shape)
	assert(left.shape[0] == left.shape[1])

	if left.shape[0] <= TRIVIAL_MULTIPLICATION_BOUND:
		return trivial_mul(left, right)

	assert(left.shape[0] % 2 == 0)
	return np.block(
		strassen_mul_2x2(*map(split_to_2x2_blocks, [left, right]))
	)

Дальнейшее развитие

Штрассен был первым, кто показал возможность умножения матриц более эффективным способом, чем стандартный. После публикации его работы в 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$ делятся на блочные подматрицы как показано выше.

Вычисляются промежуточные элементы $S_{1},\,\ldots ,\,S_{8},\,P_{1},\,\ldots ,\,P_{7},\,T_{1},\,T_{2}$

${\begin{aligned}S_{1}&=(A_{21}+A_{22});\\S_{2}&=(S_{1}-A_{11});\\S_{3}&=(A_{11}-A_{21});\\S_{4}&=(A_{12}-S_{2});\\S_{5}&=(B_{12}-B_{11});\\S_{6}&=(B_{22}-S_{5});\\S_{7}&=(B_{22}-B_{12});\\S_{8}&=(S_{6}-B_{21});\\P_{1}&=S_{2}S_{6};\\P_{2}&=A_{11}B_{11};\\P_{3}&=A_{12}B_{21};\\P_{4}&=S_{3}S_{7};\\P_{5}&=S_{1}S_{5};\\P_{6}&=S_{4}B_{22};\\P_{7}&=A_{22}S_{8};\\T_{1}&=P_{1}+P_{2};\\T_{2}&=T_{1}+P_{4}.\end{aligned}}$

Элементы матрицы $C$ вычисляются следующим образом:

${\begin{pmatrix}C_{11}&C_{12}\\C_{21}&C_{22}\end{pmatrix}}={\begin{pmatrix}P_{2}+P_{3}&T_{1}+P_{5}+P_{6}\\T_{2}-P_{7}&T_{2}+P_{5}\end{pmatrix}}.$

Современное состояние вопроса

Алгоритм Штрассена является билинейным алгоритмом, его коэффициенты являются корнями кубической системы уравнений Брента. Для класса точных алгоритмов <2x2x2> это минимальная задача, решение которой позволяет уменьшить число умножений в кольце элементов матриц. Проблема поиска новых алгоритмов заключается в том, что система Брента нелинейная, числа неизвестных и уравнений (эти числа не совпадают) быстро растет с увеличением размера матриц и требуются решения только с большим количеством нулей.

B 2013 году после частичного преодоления этих проблем удалось найти первый практически реализуемый билинейный алгоритм умножения матриц, который асимптотически быстрее алгоритма Штрассена. Алгоритм Смирнова <3x3x6; 40> умножает матрицу 3X3 на матрицу 3x6, используя 40 умножений вместо 54. Его асимптотическая сложность $O(n^{\log _{54}64000})=O(n^{2.78})$. (Тензорное умножение алгоритма на себя с циклическим сдвигом аргументов приводит к алгоритму для квадратных матриц <54x54x54; 64000> с той же сложностью). Для реального ускорения умножения требуется существенная оптимизация — удаление множества дублирующих вычислений в линейных формах.

По состоянию на 2022 год это асимптотически наиболее быстрый практически реализуемый билинейный алгоритм для произвольного поля элементов матриц.

В 2022 году компанией DeepMind при помощи нейросетевого алгоритма AlphaTensor были найдены несколько новых алгоритмов перемножения матриц различных размерностей. Однако их скорость для произвольного поля далека от скорости известных лучших алгоритмов. Так для матриц 4X4 алгоритм Штрассена требует 49 умножений, а AlphaTensor нашёл алгоритм, требующий 47 умножений, но работает он только для поля $\mathbb {Z} /2\mathbb {Z}$.

Литература

• Strassen V. Gaussian Elimination is not Optimal (англ.) // Numerische Mathematik / F. Brezzi — Springer Science+Business Media, 1969. — Vol. 13, Iss. 4. — P. 354—356. — ISSN 0029-599X; 0945-3245 — doi:10.1007/BF02165411
• Левитин А. В. Глава 4. Метод декомпозиции: Умножение больших целых чисел и алгоритм умножения матриц Штрассена // Алгоритмы. Введение в разработку и анализ — М.: Вильямс, 2006. — С. 189—195. — 576 с. — ISBN 978-5-8459-0987-9
• Кормен, Томас Х., Лейзерсон, Чарльз И., Ривест, Рональд Л., Штайн, Клифорд. Глава 28. Работа с матрицами // Алгоритмы: построение и анализ = Introduction to Algorithms. — 2-e издание. — М.: «Вильямс», 2005. — С. 833 - 839. — ISBN 5-8459-0857-4.