Алгоритм штрассена для умножения матриц

Алгоритм штрассена для умножения матриц

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

В отличие от традиционного алгоритма умножения матриц (по формуле c i , k = ∑ a i , j b j , k <displaystyle c_=sum a_b_> ), работающего за время Θ ( n log 2 ⁡ 8 ) = Θ ( n 3 ) <displaystyle Theta (n^<log _<2>8>)=Theta (n^<3>)> , алгоритм Штрассена умножает матрицы за время Θ ( n log 2 ⁡ 7 ) = O ( n 2.81 ) <displaystyle Theta (n^<log _<2>7>)=O(n^<2.81>)> , что даёт выигрыш на больших плотных матрицах начиная, примерно, от 64×64.

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

Содержание

Алгоритм [ править | править код ]

Пусть A, B — две квадратные матрицы над кольцом R. Матрица C получается по формуле:

C = A B A , B , C ∈ R 2 n × 2 n <displaystyle mathbf =mathbf mathbf qquad mathbf ,mathbf ,mathbf in R^<2^ imes 2^>>

Если размер умножаемых матриц n не является натуральной степенью двойки, мы дополняем исходные матрицы дополнительными нулевыми строками и столбцами. При этом мы получаем удобные для рекурсивного умножения размеры, но теряем в эффективности за счёт дополнительных ненужных умножений.

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

A = [ A 1 , 1 A 1 , 2 A 2 , 1 A 2 , 2 ] , B = [ B 1 , 1 B 1 , 2 B 2 , 1 B 2 , 2 ] , C = [ C 1 , 1 C 1 , 2 C 2 , 1 C 2 , 2 ] <displaystyle mathbf =<eginmathbf _<1,1>&mathbf _<1,2>\mathbf _<2,1>&mathbf _<2,2>end><mbox< , >>mathbf =<eginmathbf _<1,1>&mathbf _<1,2>\mathbf _<2,1>&mathbf _<2,2>end><mbox< , >>mathbf =<eginmathbf _<1,1>&mathbf _<1,2>\mathbf _<2,1>&mathbf _<2,2>end>>

A i , j , B i , j , C i , j ∈ R 2 n − 1 × 2 n − 1 <displaystyle mathbf _,mathbf _,mathbf _in R^<2^ imes 2^>>

C 1 , 1 = A 1 , 1 B 1 , 1 + A 1 , 2 B 2 , 1 <displaystyle mathbf _<1,1>=mathbf _<1,1>mathbf _<1,1>+mathbf _<1,2>mathbf _<2,1>> C 1 , 2 = A 1 , 1 B 1 , 2 + A 1 , 2 B 2 , 2 <displaystyle mathbf _<1,2>=mathbf _<1,1>mathbf _<1,2>+mathbf _<1,2>mathbf _<2,2>> C 2 , 1 = A 2 , 1 B 1 , 1 + A 2 , 2 B 2 , 1 <displaystyle mathbf _<2,1>=mathbf _<2,1>mathbf _<1,1>+mathbf _<2,2>mathbf _<2,1>> C 2 , 2 = A 2 , 1 B 1 , 2 + A 2 , 2 B 2 , 2 <displaystyle mathbf _<2,2>=mathbf _<2,1>mathbf _<1,2>+mathbf _<2,2>mathbf _<2,2>>

Однако с помощью этой процедуры нам не удалось уменьшить количество умножений. Как и в обычном методе, нам требуется 8 умножений.

Теперь определим новые элементы

P 1 := ( A 1 , 1 + A 2 , 2 ) ( B 1 , 1 + B 2 , 2 ) <displaystyle mathbf

_<1>:=(mathbf _<1,1>+mathbf _<2,2>)(mathbf _<1,1>+mathbf _<2,2>)> P 2 := ( A 2 , 1 + A 2 , 2 ) B 1 , 1 <displaystyle mathbf

_<2>:=(mathbf _<2,1>+mathbf _<2,2>)mathbf _<1,1>> P 3 := A 1 , 1 ( B 1 , 2 − B 2 , 2 ) <displaystyle mathbf

_<3>:=mathbf _<1,1>(mathbf _<1,2>-mathbf _<2,2>)> P 4 := A 2 , 2 ( B 2 , 1 − B 1 , 1 ) <displaystyle mathbf

_<4>:=mathbf _<2,2>(mathbf _<2,1>-mathbf _<1,1>)> P 5 := ( A 1 , 1 + A 1 , 2 ) B 2 , 2 <displaystyle mathbf

_<5>:=(mathbf _<1,1>+mathbf _<1,2>)mathbf _<2,2>> P 6 := ( A 2 , 1 − A 1 , 1 ) ( B 1 , 1 + B 1 , 2 ) <displaystyle mathbf

_<6>:=(mathbf _<2,1>-mathbf _<1,1>)(mathbf _<1,1>+mathbf _<1,2>)> P 7 := ( A 1 , 2 − A 2 , 2 ) ( B 2 , 1 + B 2 , 2 ) <displaystyle mathbf

_<7>:=(mathbf _<1,2>-mathbf _<2,2>)(mathbf _<2,1>+mathbf _<2,2>)>

которые затем используются для выражения Ci, j. Таким образом, нам нужно всего 7 умножений на каждом этапе рекурсии. Элементы матрицы C выражаются из Pk по формулам

C 1 , 1 = P 1 + P 4 − P 5 + P 7 <displaystyle mathbf _<1,1>=mathbf

_<1>+mathbf

_<4>-mathbf

_<5>+mathbf

_<7>> C 1 , 2 = P 3 + P 5 <displaystyle mathbf _<1,2>=mathbf

_<3>+mathbf

_<5>> C 2 , 1 = P 2 + P 4 <displaystyle mathbf _<2,1>=mathbf

_<2>+mathbf

_<4>> C 2 , 2 = P 1 − P 2 + P 3 + P 6 <displaystyle mathbf _<2,2>=mathbf

_<1>-mathbf

_<2>+mathbf

_<3>+mathbf

_<6>>

Рекурсивный процесс продолжается n раз, до тех пор пока размер матриц Ci, j не станет достаточно малым, далее используют обычный метод умножения матриц. Это делают из-за того, что алгоритм Штрассена теряет эффективность по сравнению с обычным на малых матрицах в силу большего числа сложений. Оптимальный размер матриц для перехода к обычному умножению зависит от характеристик процессора и языка программирования и на практике лежит в пределах от 32 до 128.

Читайте также:  The crew оптимизация для слабых пк

Пример реализации на Фортране [ править | править код ]

Приведён пример реализации алгоритма Штрассена на Фортране. Предполагается, что все матрицы квадратные, их размер является степенью 2.

Вычисления промежуточных матриц P1P7 можно проводить параллельно при использовании таких библиотек как OpenMP или MPI, что позволяет значительно повысить скорость работы алгоритма на многопроцессорных машинах.

Дальнейшее развитие [ править | править код ]

Штрассен был первым, кто показал возможность умножения матриц более эффективным способом, чем стандартный. После публикации его работы в 1969 начались активные поиски более быстрого алгоритма. Самым асимптотически быстрым алгоритмом на сегодняшний день является алгоритм Копперсмита — Винограда, позволяющий перемножать матрицы за O ( n 2.376 ) <displaystyle <
m >(n^<2.376>)> операций [1] , предложенный в 1987 году и усовершенствованный в 2011 году до уровня O ( n 2.3727 ) <displaystyle <
m
>(n^<2.3727>)> [1] . Этот алгоритм не представляет практического интереса в силу астрономически большой константы в оценке арифметической сложности. Вопрос о предельной в асимптотике скорости перемножения матриц не решен. Существует гипотеза Штрассена о том, что для достаточно больших n существует алгоритм перемножения двух матриц размера n × n <displaystyle n imes n> за O ( n 2 + ε ) <displaystyle <
m
>(n^<2+varepsilon >)> операций, где ε <displaystyle varepsilon > сколь угодно малое наперед заданное положительное число. Эта гипотеза имеет сугубо теоретический интерес, так как размер матриц, для которых ε <displaystyle varepsilon > действительно мало, по всей видимости очень велик.

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

Алгоритм Винограда — Штрассена [ править | править код ]

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

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

Вычисляются промежуточные матрицы S1S8, P1P7, T1T2

S 1 := ( A 2 , 1 + A 2 , 2 ) <displaystyle mathbf _<1>:=(mathbf _<2,1>+mathbf _<2,2>)> S 2 := ( S 1 − A 1 , 1 ) <displaystyle mathbf _<2>:=(mathbf _<1>-mathbf _<1,1>)> S 3 := ( A 1 , 1 − A 2 , 1 ) <displaystyle mathbf _<3>:=(mathbf _<1,1>-mathbf _<2,1>)> S 4 := ( A 1 , 2 − S 2 ) <displaystyle mathbf _<4>:=(mathbf _<1,2>-mathbf _<2>)> S 5 := ( B 1 , 2 − B 1 , 1 ) <displaystyle mathbf _<5>:=(mathbf _<1,2>-mathbf _<1,1>)> S 6 := ( B 2 , 2 − S 5 ) <displaystyle mathbf _<6>:=(mathbf _<2,2>-mathbf _<5>)> S 7 := ( B 2 , 2 − B 1 , 2 ) <displaystyle mathbf _<7>:=(mathbf _<2,2>-mathbf _<1,2>)> S 8 := ( S 6 − B 2 , 1 ) <displaystyle mathbf _<8>:=(mathbf _<6>-mathbf _<2,1>)>

P 1 := S 2 S 6 <displaystyle mathbf

_<1>:=mathbf _<2>mathbf _<6>> P 2 := A 1 , 1 B 1 , 1 <displaystyle mathbf

_<2>:=mathbf _<1,1>mathbf _<1,1>> P 3 := A 1 , 2 B 2 , 1 <displaystyle mathbf

_<3>:=mathbf _<1,2>mathbf _<2,1>> P 4 := S 3 S 7 <displaystyle mathbf

_<4>:=mathbf _<3>mathbf _<7>> P 5 := S 1 S 5 <displaystyle mathbf

Читайте также:  Cameleo ru одноклассники официальный сайт

_<5>:=mathbf _<1>mathbf _<5>> P 6 := S 4 B 2 , 2 <displaystyle mathbf

_<6>:=mathbf _<4>mathbf _<2,2>> P 7 := A 2 , 2 S 8 <displaystyle mathbf

_<7>:=mathbf _<2,2>mathbf _<8>>

T 1 := P 1 + P 2 <displaystyle mathbf _<1>:=mathbf

_<1>+mathbf

_<2>> T 2 := T 1 + P 4 <displaystyle mathbf _<2>:=mathbf _<1>+mathbf

_<4>>

Элементы матрицы C вычисляются по формулам

C 1 , 1 := P 2 + P 3 <displaystyle mathbf _<1,1>:=mathbf

_<2>+mathbf

_<3>> C 1 , 2 := T 1 + P 5 + P 6 <displaystyle mathbf _<1,2>:=mathbf _<1>+mathbf

_<5>+mathbf

_<6>> C 2 , 1 := T 2 − P 7 <displaystyle mathbf _<2,1>:=mathbf _<2>-mathbf

_<7>> C 2 , 2 := T 2 + P 5 <displaystyle mathbf _<2,2>:=mathbf _<2>+mathbf

_<5>>

Метод матричного умножения Штрассена является типичным алгоритмом «разделяй и властвуй». Мы обсудили алгоритм Штрассена здесь . Однако давайте еще раз рассмотрим, что стоит за подходом «разделяй и властвуй», и реализуем его.

Предварительное условие: перед дальнейшим пониманием обязательно посмотрите этот пост.

Реализация

// Программа CPP для реализации матрицы Штрассена
// Алгоритм умножения
#include

using namespace std;

typedef long long lld;

/ * Алгоритм Штрассена для умножения матриц

Сложность: O (n ^ 2.808) * /

inline lld** MatrixMultiply(lld** a, lld** b, int n,

lld** c = new lld*[n];

for ( int i = 0; i

for ( int i = 0; i

for ( int j = 0; j

for ( int k = 0; k

inline lld** Strassen(lld** a, lld** b, int n,

if (n == 1 || l == 1 || m == 1)

return MatrixMultiply(a, b, n, l, m);

lld** c = new lld*[n];

for ( int i = 0; i

int adjN = (n >> 1) + (n & 1);

int adjL = (l >> 1) + (l & 1);

int adjM = (m >> 1) + (m & 1);

lld**** As = new lld***[2];

for ( int x = 0; x

for ( int y = 0; y

As[x][y] = new lld*[adjN];

for ( int i = 0; i

As[x][y][i] = new lld[adjL];

for ( int j = 0; j

lld**** Bs = new lld***[2];

for ( int x = 0; x

for ( int y = 0; y

Bs[x][y] = new lld*[adjN];

for ( int i = 0; i

Bs[x][y][i] = new lld[adjM];

for ( int j = 0; j

lld*** s = new lld**[10];

for ( int i = 0; i

s[i] = new lld*[adjL];

for ( int j = 0; j

s[i][j] = new lld[adjM];

for ( int k = 0; k

s[i] = new lld*[adjN];

for ( int j = 0; j

s[i][j] = new lld[adjL];

for ( int k = 0; k

s[i] = new lld*[adjN];

for ( int j = 0; j

s[i][j] = new lld[adjL];

for ( int k = 0; k

s[i] = new lld*[adjL];

for ( int j = 0; j

s[i][j] = new lld[adjM];

for ( int k = 0; k

s[i] = new lld*[adjN];

for ( int j = 0; j

s[i][j] = new lld[adjL];

for ( int k = 0; k

s[i] = new lld*[adjL];

for ( int j = 0; j

s[i][j] = new lld[adjM];

for ( int k = 0; k

s[i] = new lld*[adjN];

for ( int j = 0; j

s[i][j] = new lld[adjL];

for ( int k = 0; k

s[i] = new lld*[adjL];

for ( int j = 0; j

s[i][j] = new lld[adjM];

for ( int k = 0; k

s[i] = new lld*[adjN];

for ( int j = 0; j

s[i][j] = new lld[adjL];

for ( int k = 0; k

s[i] = new lld*[adjL];

for ( int j = 0; j

s[i][j] = new lld[adjM];

for ( int k = 0; k

lld*** p = new lld**[7];

p[0] = Strassen(As[0][0], s[0], adjN, adjL, adjM);

p[1] = Strassen(s[1], Bs[1][1], adjN, adjL, adjM);

p[2] = Strassen(s[2], Bs[0][0], adjN, adjL, adjM);

p[3] = Strassen(As[1][1], s[3], adjN, adjL, adjM);

Читайте также:  Javascript текст в число

p[4] = Strassen(s[4], s[5], adjN, adjL, adjM);

p[5] = Strassen(s[6], s[7], adjN, adjL, adjM);

p[6] = Strassen(s[8], s[9], adjN, adjL, adjM);

for ( int i = 0; i

for ( int j = 0; j

c[i][j + adjM] = p[0][i][j] + p[1][i][j];

c[i + adjN][j] = p[2][i][j] + p[3][i][j];

c[i + adjN][j + adjM] = p[4][i][j] + p[0][i][j] — p[2][i][j] — p[6][i][j];

Матрица — математический объект, эквивалентный двумерному массиву. Числа располагаются в матрице по строкам и столбцам. Числа располагаются в матрице по строкам и столбцам. Две матрицы одинакового размера можно поэлементно сложить или вычесть друг их друга. Если число столбцов в первой матрице совпадает с числом строк во второй, то эти две матрицы можно перемножить. У произведения будет столько же строк, сколько в первой матрице, и столько же столбцов, сколько во второй. При умножении матрицы размером 3 х 4 на матрицу размером 4 х 7 мы получаем матрицу размером 3 х 7. Умножение матриц некоммутативно: оба произведения AB и BA двух квадратных матриц одинакового размера можно вычислить, однако результаты, вообще говоря, будут отличаться друг от друга. (Отметим, что умножение чисел коммутативно, и произведения AB и BA двух чисел A и B равны.)
Для вычисления произведения двух матриц каждая строка первой почленно умножается на каждый столбец второй. Затем подсчитывается сумма таких произведений и записывается в соответствующую клетку результата. На рисунке 1 приведен пример умножения двух матриц.
Умножение матриц на рисунке 1 требует 24 умножений и 16 сложений. Вообще, стандартный алгоритм умножения матрицы размером a х b на матрицу размером b х c выполняет abc умножений и a(b — 1)c сложений.
Вот как выглядит стандартный алгоритм умножения матрицы G размером a x b на матрицу H размером b x c. Результат записывается в матрицу R размером a x c. На первый взгляд это минимальный объем работы, необходимый для перемножения двух матриц. Однако исследователям не удалось доказать минимальность, и в результате они обнаружили другие алгоритмы, умножающие матрицы более эффективно.

Рассмотрим случай четной общей размерности b. Результаты подсчета сложений и умножений сведены в следующую таблицу.

Предварительная обработка матрицы G Предварительная обработка матрицы H Вычисление матрицы R Всего

Алгоритм Штрассена работает с квадратными матрицами. На самом деле он настолько эффективен, что иногда разумно расширить матрицы до квадратных, и при этом он все равно дает выигрыш, превышающий расходы на введение дополнительных элементов.
Для умножения матриц алгоритм Штрассена использует семь формул. Эти формулы чрезвычайно неестественны и, к сожалению, в оригинальной статье Штрассена не объясняется, как он пришел к ним. Замечательно, что как сами формулы, так и их использование не требуют, чтобы умножение элементов матриц было коммутативным. Это означает, в частности, что сами эти элементы могут быть матрицами, а значит, алгоритм Штрассена можно применять рекурсивно. Вот формулы Штрассена:

Ссылка на основную публикацию
Adblock detector