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

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

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

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

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

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

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

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

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

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

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

${\displaystyle 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, right = 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 начались активные поиски более быстрого алгоритма. Самым асимптотически быстрым алгоритмом на сегодняшний день является алгоритм Копперсмита — Винограда, позволяющий перемножать матрицы за ${\displaystyle {\rm {O}}(n^{2.376})}$ операций[1], предложенный в 1987 году и усовершенствованный в 2011 году до уровня ${\displaystyle {\rm {O}}(n^{2.3727})}$[1]. Этот алгоритм не представляет практического интереса в силу астрономически большой константы в оценке арифметической сложности. Вопрос о предельной в асимптотике скорости перемножения матриц не решен. Существует гипотеза Штрассена о том, что для достаточно больших ${\displaystyle n}$ существует алгоритм перемножения двух матриц размера ${\displaystyle n\times n}$ за ${\displaystyle {\rm {O}}(n^{2+\varepsilon })}$ операций, где ${\displaystyle \varepsilon }$ сколь угодно малое наперед заданное положительное число. Эта гипотеза имеет сугубо теоретический интерес, так как размер матриц, для которых ${\displaystyle \varepsilon }$ действительно мало, по всей видимости очень велик.

Вопрос о построении наиболее быстрого и устойчивого практического алгоритма умножения больших матриц также остается нерешённым.

Существует модификация алгоритма Штрассена, для которой требуется 7 умножений и 15 сложений (вместо 18 для обычного алгоритма Штрассена).

Матрицы ${\displaystyle A,\,B,\,C}$ делятся на блочные подматрицы как показано выше.

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

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

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

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

• Strassen V. Gaussian Elimination is not Optimal (англ.) // Numer. Math / 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.