Участник:Y.ryabkova/Умножение матриц методом Штрассена: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
Строка 79: Строка 79:
 
== Вычислительное ядро алгоритма ==
 
== Вычислительное ядро алгоритма ==
  
Вычислительное ядро последовательного варианта алгоритма Штрассена на каждом рекурсивном шаге, кроме последнего, можно составить из вычисления матриц <math> S_{i}, T_{i}, \: i = \overline{1, 7} </math>, и блоков матрицы <math> C </math> из матриц <math> \alpha_{i}, \: i = \overline{1, 7} </math> по формулам, данным в математическом описании алгоритма. Все эти вычисления представляют из себя <math> 18 </math> операций сложения (вычитания) матриц порядка <math> k </math>, зависящего от номера шага <math> s_m </math> следующим образом: <math> k = \frac{n}{2^{s_m}} </math>. Операции выполняются поэлементно. Таким образом, вычислительное ядро алгоритма сводится к <math> 18 </math> операциям сложения (вычитания) матриц или к <math> 18 k^2 </math> операциям сложения (вычитания) чисел на каждом рекурсивном шаге.
+
В описанном алгоритме можно выделить два вычислительных ядра: первое соответствует операциям, выполняемым на первых <math> m </math> рекурсивных шагах, а второе - последнему рекурсивному шагу с номером <math> m + 1</math>.
 +
 
 +
Вычислительное ядро последовательного варианта алгоритма Штрассена на каждом рекурсивном шаге, кроме последнего, можно составить из вычисления матриц <math> S_{i}, T_{i}, \: i = \overline{1, 7} </math>, и блоков матрицы <math> C </math> из матриц <math> \alpha_{i}, \: i = \overline{1, 7} </math> по формулам, данным в математическом описании алгоритма. Все эти вычисления представляют из себя сложения (вычитания) матриц порядка <math> k </math>, зависящего от номера шага <math> s_m </math> следующим образом: <math> k = \frac{n}{2^{s_m}} </math>. Операции выполняются поэлементно. Таким образом, первое вычислительное ядро алгоритма сводится к операциям сложения (вычитания) чисел на каждом рекурсивном шаге.
 +
 
 +
Вычислительное ядро алгоритма на последнем, <math> m+1 </math>-ом рекурсивном шаге, состоит из операций перемножения двух матриц по классическому алгоритму кубической сложности, основную часть которого составляют умножения чисел.
 +
 
 +
Указанные операции составляют основную часть метода, причем с увеличением <math> m </math> доля операций первого вычислительного ядра будет расти, а доля операций второго - падать, и наоборот.
 +
 
  
 
== Макроструктура алгоритма ==
 
== Макроструктура алгоритма ==
Строка 85: Строка 92:
 
Алгоритм Штрассена состоит из <math> m+1 </math> рекурсивных шагов и может быть представлен в виде двух этапов.  
 
Алгоритм Штрассена состоит из <math> m+1 </math> рекурсивных шагов и может быть представлен в виде двух этапов.  
  
Первый этап соответствует первым <math> m </math> рекурсивным шагам, на каждом из которых можно выделить следующие три последовательные макрооперации:
+
'''Первый этап''' соответствует первым <math> m </math> рекурсивным шагам, на каждом из которых можно выделить следующие три последовательные макрооперации:
1. вычисление матриц <math> S_{i}, T_{i}, \: i = \overline{1, 7} </math>;
+
 
2. рекурсивный вызов для вычисления произведений <math> \alpha_{i} = S_{i} \cdot T_{i}, \: i = \overline{1, 7} </math>, таким образом при переходе на следующий рекурсивный шаг задача размера <math> k </math> сводится к семи задачам размера <math> \frac{k}{2} </math>;
+
* вычисление матриц <math> S_{i}, T_{i}, \: i = \overline{1, 7} </math>;
3. вычисление матрицы <math> C </math> из матриц <math> \alpha_{i} </math>.
+
 
 +
* рекурсивный вызов для вычисления произведений <math> \alpha_{i} = S_{i} \cdot T_{i}, \: i = \overline{1, 7} </math>, таким образом при переходе на следующий рекурсивный шаг задача размера <math> k </math> сводится к семи задачам размера <math> \frac{k}{2} </math>;
 +
 
 +
* вычисление матрицы <math> C </math> из матриц <math> \alpha_{i} </math>.
  
Второй этап соответствует последнему, <math> m+1 </math>-ому шагу алгоритма и выходу из рекурсии. Он представляет из себя вычисление произведения двух матриц матриц порядка <math> n_{m} = \frac{n}{2^{m}} </math> по классическому алгоритму: <br>
+
'''Второй этап''' соответствует последнему, <math> m+1 </math>-ому шагу алгоритма и выходу из рекурсии. Он представляет из себя вычисление произведения двух матриц матриц порядка <math> n_{m} = \frac{n}{2^{m}} </math> по классическому алгоритму: <br>
:<math> c_{ij} = \sum_{k = 1}^{n_{m}} a_{ik} b_{kj}, \: i,j = \overline{n_{m}} </math>.
+
:<math> c_{ij} = \sum_{k = 1}^{n_{m}} a_{ik} b_{kj}, \: i,j = \overline{1, n_{m}} </math>.
  
 
== Схема реализации последовательного алгоритма ==
 
== Схема реализации последовательного алгоритма ==
Строка 98: Строка 108:
  
 
<source lang="fortran">
 
<source lang="fortran">
        call Strassen(a, b, c, n, m)
 
  
 
Strassen(A, B, C, k, m) ! C = A * B
 
Strassen(A, B, C, k, m) ! C = A * B
           if (m == 0) call mul(A, B, C, k) ! классический алгоритм умножения матриц
+
           if (m == 0) call mul(A, B, C, k) ! переход к классическому алгоритму умножения матриц
 
           do i = 1, 7
 
           do i = 1, 7
 
               compute S[i], T[i] from A and B
 
               compute S[i], T[i] from A and B
Строка 112: Строка 121:
 
           end do
 
           end do
 
         end Strassen
 
         end Strassen
 +
 
</source>
 
</source>
  
 
Таким образом, на каждом рекурсивном шаге эта схема состоит из трёх последовательных этапов:
 
Таким образом, на каждом рекурсивном шаге эта схема состоит из трёх последовательных этапов:
  
1. Вычисление семи пар матриц <math> S_{i}, T_{i} </math>. Так как они определены независимо друг от друга, порядок их вычисления может быть произвольным, в том числе параллельным.
+
* Вычисление семи пар матриц <math> S_{i}, T_{i} </math>. Так как они определены независимо друг от друга и хранятся в разных областях памяти, порядок их вычисления может быть произвольным, в том числе параллельным.
  
2. Выполнение семи рекурсивных вызовов функции Strassen -- по одному для каждой тройки <math> \alpha_{i} = S_{i} \cdot T_{i}</math>. Порядок этих вызовов также может быть произвольным и параллельным.
+
* Выполнение семи рекурсивных вызовов функции ''Strassen'', по одному для каждой тройки <math> \alpha_{i} = S_{i} \cdot T_{i}</math>. Порядок этих вызовов также может быть произвольным и параллельным.
  
3. Вычисление матрицы <math> С </math> из полученных после рекурсивных вызовов матриц <math> \alpha_{i} </math>.
+
* Вычисление матрицы <math> С </math> из полученных после рекурсивных вызовов матриц <math> \alpha_{i} </math>.
  
Как видно, данная реализация требует большого количества дополнительной памяти, так как на каждом рекурсивном шаге требуется одновременно хранить семь пар <math>\frac{k}{2}\times\frac{k}{2}</math>-матриц <math> S_{i}, T_{i} </math> и семь результатов их произведений <math> \alpha_{i} </math>.
+
Как видно, данная реализация требует большого количества дополнительной памяти, так как на каждом рекурсивном шаге необходимо одновременно хранить семь пар <math>\frac{k}{2}\times\frac{k}{2}</math>-матриц <math> S_{i}, T_{i} </math> и семь результатов их умножения <math> \alpha_{i} </math>.
  
 
Есть и другой вариант реализации последовательного алгоритма, позволяющий значительно сократить количество требуемой памяти:
 
Есть и другой вариант реализации последовательного алгоритма, позволяющий значительно сократить количество требуемой памяти:
  
 
<source lang="fortran">
 
<source lang="fortran">
        call Strassen(a, b, c, n, m)
 
  
 
Strassen(A, B, C, k, m) ! C = A * B
 
Strassen(A, B, C, k, m) ! C = A * B
           if (m == 0) call mul(A, B, C, k) ! классический алгоритм умножения матриц
+
           if (m == 0) call mul(A, B, C, k) ! переход к классическому алгоритму умножения матриц
 
           do i = 1, 7
 
           do i = 1, 7
 
               compute S[i], T[i] from A and B
 
               compute S[i], T[i] from A and B
Строка 137: Строка 146:
 
           end do
 
           end do
 
         end Strassen
 
         end Strassen
 +
 
</source>
 
</source>
  
 
На каждом рекурсивном шаге схема этой реализации представляет собой цикл по <math>i = 1,\dotsc,7</math>, каждая итерация которого состоит из трёх последовательных этапов:
 
На каждом рекурсивном шаге схема этой реализации представляет собой цикл по <math>i = 1,\dotsc,7</math>, каждая итерация которого состоит из трёх последовательных этапов:
  
1. Вычисление <math>i</math>-ой пары матриц <math> S_{i}, T_{i} </math>.
+
* Вычисление <math>i</math>-ой пары матриц <math> S_{i}, T_{i} </math>.
  
2. Выполнение рекурсивного вызова функции Strassen для вычисления <math> \alpha_{i} = S_{i} \cdot T_{i} </math>.
+
* Выполнение рекурсивного вызова функции ''Strassen'' для вычисления <math> \alpha_{i} = S_{i} \cdot T_{i} </math>.
  
3. Изменение значений блоков матрицы <math> С </math> в соответствии с формулами алгоритма в зависимости от номера полученной после рекурсивного вызова матрицы <math> \alpha_{i} </math>.
+
* Обновление блоков матрицы <math> С </math> в соответствии с формулами алгоритма в зависимости от номера полученной после рекурсивного вызова матрицы <math> \alpha_{i} </math>.
  
Порядок итераций такого цикла может быть произвольным, но они должны выполняться строго последовательно. Такое требование позволяет выделить память всего для трёх <math>\frac{k}{2}\times\frac{k}{2}</math>-матриц на каждом рекурсивном шаге: для одной пары <math> S_{i}, T_{i} </math> и результата их произведения <math> \alpha_{i} </math>, и на каждой итерации цикла использовать одни и те же области памяти для хранения соответствующей тройки матриц.
+
Порядок итераций такого цикла может быть произвольным, но они должны выполняться строго последовательно. Такое требование позволяет выделить память всего для трёх <math>\frac{k}{2}\times\frac{k}{2}</math>-матриц на каждом рекурсивном шаге: для одной пары <math> S_{i}, T_{i} </math> и их произведения <math> \alpha_{i} </math>, и на каждой итерации цикла использовать одни и те же области памяти для хранения соответствующей тройки матриц.
  
 
== Последовательная сложность алгоритма ==
 
== Последовательная сложность алгоритма ==
 +
 +
Оба вышеописанных варианта реализации метода Штрассена требуют одинакового числа арифметических операций. В соответствии с макроструктурой алгоритма, рассмотрим каждый из его двух этапов отдельно.
 +
 +
'''Первый этап''' требует на первом рекурсивном шаге выполнения <math>18 \left( \frac{n}{2} \right)^2 </math> сложений (вычитаний). При переходе на следующий рекурсивный шаг, задача будет сведена к семи аналогичным задачам, но вдвое меньшей размерности. Отсюда получаем, что на <math>k</math>-ом шаге потребуется <math>18 \cdot 7^{k-1} \left( \frac{n}{2^{k}} \right)^2 </math> операций. Тогда всего на первых <math> m </math> шагах потребуется следующее число сложений (вычитаний):<br>
 +
:<math> 18 \left( \frac{n}{2} \right)^2 + 18 \cdot 7 \left( \frac{n}{4} \right)^2 + \dotsb + 18 \cdot 7^{m-1} \left( \frac{n}{2^{m}} \right)^2 = \frac{18}{4}n^{2}\left( 1 + \frac{7}{4} + \dotsb + \frac{7^{m-1}}{4^{m-1}}\right) = 6n^2\left(\frac{7^{m}}{4^{m}} - 1 \right) </math>.
 +
Получена сложность первого этапа алгоритма.
 +
 +
'''Второй этап''' включает в себя решение <math>7^m</math> задач перемножения матриц порядка <math> \frac{n}{2^{m}} </math> по классическому алгоритму, что требует <math>7^m \left(\frac{n}{2^{m}}\right)^3</math> умножений и <math>7^m \left(\frac{n}{2^{m}}\right)^2\left(\frac{n}{2^{m}}-1\right)</math> сложений, всего <math>O\left(\frac{n^3}{2^{3m}}\right)</math> арифметических операций.
 +
 +
'''Всего''' последовательный алгоритм Штрассена требует:
 +
* <math>6n^2\left(\frac{7^{m}}{4^{m}} - 1 \right) + 7^m\left(\frac{n}{2^{m}}\right)^2\left(\frac{n}{2^{m}}-1\right)</math> сложений;
 +
* <math>7^m\left(\frac{n}{2^{m}}\right)^3</math> умножений.
 +
Как видно, сложность алгоритма зависит не только от порядка матрицы <math> n </math>, но и от числа <math> m </math> рекурсивных шагов первого этапа.
 +
 +
Если <math> m = L</math>, где <math>n = 2^L</math>, то есть второй этап алгоритма сводится лишь к <math>7^L</math> умножениям двух чисел, то сложность всего последовательного алгоритма равна <br>
 +
:<math> 6n^2\left(\frac{7^{L}}{4^{L}} - 1 \right) + 7^L\left(\frac{n}{2^{L}}\right)^3 = 6n^2\left(\frac{7^{L}}{n^{2}} - 1 \right) + 7^L = 7 \cdot 7^{L} - 6n^2 = 7 n^{\log_2 7} - 6n^2 = O\left(n^{\log_2 7}\right) \approx O\left(n^{2.81}\right) </math>.
 +
Таким образом, получена оценка сложности последовательного алгоритма Штрассена с максимально возможным числом <math>m = L</math> рекурсивных шагов, соответствующая максимальной оценке сложности первого этапа алгоритма для произвольного <math>m</math>. Эта оценка асимптотически лучше оценки для классического метода перемножения матриц, требующего <math>O(n^3)</math> операций.
 +
 +
Стоит отметить роль правильного выбора параметра <math> m </math>. Из полученных оценок видно, что при малых <math> m </math> сложность алгоритма с уменьшением этого параметра будет стремиться к кубической сложности классического метода, однако при <math> m </math>, очень близких к <math> L </math>, алгоритм будет требовать большего числа сложений, чем классический метод, которое будет увеличиваться с ростом <math> m </math>. Отсюда следует, что существует некоторое оптимальное значение <math> m \in [1, L] </math>, которое минимизирует число арифметических операций всего алгоритма. Также при выборе <math> m </math> стоит учитывать роль накладных расходов на реализацию рекурсии при малых размерах матриц.
  
 
== Информационный граф ==
 
== Информационный граф ==

Версия 10:52, 30 ноября 2017

Основные авторы описания: Ю.В.Рябкова

1 Свойства и структура алгоритмов

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

Алгоритм Штрассена, описанный в 1969 году Фолькером Штрассеном, предназначен для умножения двух произвольных матриц при условии согласованности их размеров. В данном варианте алгоритма будет рассмотрен только частный случай умножения квадратных матриц порядка n над полем вещественных чисел \mathbb{R}.

В 1965 году Штрассен нашел способ умножения двух матриц размера 2 \times 2 с использованием семи умножений[1], в то время как классический алгоритм перемножения матриц требует восьми умножений. Стоит отметить, что данный способ не использует свойство коммутативности для операции перемножения чисел.

Теперь предположим, что n = 2^L, L \in \mathbb{N} (далее будет показано, что это предположение не нарушает общности рассуждений), и будем рассматривать матрицы A и B как блочные 2 \times 2-матрицы с размером блоков \frac{n}{2} \times \frac{n}{2} . Так как для блочных матриц правила умножения остаются теми же, как если бы на месте блока стояло число, а в методе Штрассена умножения 2 \times 2-матриц не используется коммутативность, этот метод подходит и для умножения блочных 2 \times 2-матриц. Таким образом, от задачи умножения 2 \times 2-матриц с семью умножениями можно перейти к задаче умножения n \times n-матриц, требующей не более чем O(7n^{\log_27}) операций, что асимптотически лучше классического способа умножения матриц сложностью O(n^3). На данный момент известны и более быстрые методы умножения матриц, однако алгоритм Штрассена относительно прост в реализации и эффективен на плотных матрицах не очень больших размеров, поэтому всё ещё широко используется.

1.2 Математическое описание алгоритма

Пусть A, B \in \mathbb{R}^{n \times n} — заданные плотные матрицы. Вычислим матрицу C \in \mathbb{R}^{n \times n} , такую, что C = AB .

Для удобства будем полагать, что n = 2^L, L \in \mathbb{N} . Это предположение не ограничивает общности рассуждений, так как если это свойство не выполнено и n \lt 2^L , можно рассмотреть квадратные матрицы \tilde{A} и \tilde{B} порядка 2^L следующего вида:

\tilde{A} = \begin{bmatrix} A & 0 \\ 0 & 0 \end{bmatrix}, \: \: \tilde{B} = \begin{bmatrix} B & 0 \\ 0 & 0 \end{bmatrix},

то есть матрицы, полученные из матриц A и B дополнением их нулевыми строками и столбцами в позициях i = \overline{n+1, 2^L} до размера 2^L \times 2^L . По правилу перемножения блочных матриц:

\tilde{C} = \tilde{A}\tilde{B} = \begin{bmatrix} AB & 0 \\ 0 & 0 \end{bmatrix} = \begin{bmatrix} C & 0 \\ 0 & 0 \end{bmatrix}.

Таким образом, решая задачу для матриц \tilde{A} и \tilde{B} , получим решение исходной задачи умножения матриц A и B как левый верхний n \times n-блок матрицы \tilde{C}. Аналогичным образом можно не ограничивая общности предположить, что n = s \cdot 2^L, \: L, s \in \mathbb{N} .

Итак, пусть n = 2^L, L \in \mathbb{N} . Рассмотрим матрицы A и B как блочные 2 \times 2 -матрицы с квадратными блоками порядка \frac{n}{2}:

A = \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}, \: \: B = \begin{bmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{bmatrix}.

Их произведение C = AB , согласно алгоритму Штрассена, вычисляется следующим образом:

\begin{align*} &S_{1} = A_{11} + A_{22}, & &T_{1} = B_{11} + B_{22}, \\ &S_{2} = A_{21} + A_{22}, & &T_{2} = B_{11}, \\ &S_{3} = A_{11}, & &T_{3} = B_{12} − B_{22}, \\ &S_{4} = A_{22}, & &T_{4} = B_{21} − B_{11}, \\ &S_{5} = A_{11} + A_{12}, & &T_{5} = B_{22}, \\ &S_{6} = A_{21} − A_{11}, & &T_{6} = B_{11} + B_{12}, \\ &S_{7} = A_{12} − A_{22}, & &T_{7} = B_{21} + B_{22}, \end{align*}
\alpha_{i} = S_{i} \cdot T_{i}, \: i = \overline{1,7},
C_{11} = \alpha_{1} + \alpha_{4} - \alpha_{5} + \alpha_{7},
C_{12} = \alpha_{3} + \alpha_{5},
C_{21} = \alpha_{2} + \alpha_{4},
C_{22} = \alpha_{1} + \alpha_{3} - \alpha_{2} + \alpha_{6},
C = \begin{bmatrix} C_{11} & C_{12} \\ C_{21} & C_{22} \end{bmatrix}.

Для вычисления матриц \alpha_{i}, i = \overline{1,7} , также требуется выполнить семь матричных умножений по формулам, данным выше, но уже для матриц порядка \frac{n}{2} . Таким образом, получаем рекурсивный алгоритм, на каждом шаге рекурсии требующий семи умножений и сводящий задачу размера k к задаче размера \frac{k}{2}. Всего потребуется \log_{2}n рекурсивных шагов, на последнем из которых задача будет сведена к перемножению двух чисел. Однако, на практике разумно продолжать рекурсивный процесс до тех пор, пока k не станет достаточно малым (или, в случае n = s \cdot 2^L , пока k кратно 2), а далее использовать классический способ умножения матриц. Именно такой вариант метода рассматривается в данном описании. Это связано с тем, что на малых матрицах алгоритм Штрассена теряет эффективность в сравнении с классическим алгоритмом умножения матриц, так как требует большого числа сложений.

Итак, пусть весь процесс перемножения матриц требует m+1 рекурсивных шагов, где первые m шагов выполняются в соответствии с алгоритмом Штрассена, а на последнем шаге используется классический алгоритм умножения матриц.

1.3 Вычислительное ядро алгоритма

В описанном алгоритме можно выделить два вычислительных ядра: первое соответствует операциям, выполняемым на первых m рекурсивных шагах, а второе - последнему рекурсивному шагу с номером m + 1.

Вычислительное ядро последовательного варианта алгоритма Штрассена на каждом рекурсивном шаге, кроме последнего, можно составить из вычисления матриц S_{i}, T_{i}, \: i = \overline{1, 7} , и блоков матрицы C из матриц \alpha_{i}, \: i = \overline{1, 7} по формулам, данным в математическом описании алгоритма. Все эти вычисления представляют из себя сложения (вычитания) матриц порядка k , зависящего от номера шага s_m следующим образом: k = \frac{n}{2^{s_m}} . Операции выполняются поэлементно. Таким образом, первое вычислительное ядро алгоритма сводится к операциям сложения (вычитания) чисел на каждом рекурсивном шаге.

Вычислительное ядро алгоритма на последнем, m+1 -ом рекурсивном шаге, состоит из операций перемножения двух матриц по классическому алгоритму кубической сложности, основную часть которого составляют умножения чисел.

Указанные операции составляют основную часть метода, причем с увеличением m доля операций первого вычислительного ядра будет расти, а доля операций второго - падать, и наоборот.


1.4 Макроструктура алгоритма

Алгоритм Штрассена состоит из m+1 рекурсивных шагов и может быть представлен в виде двух этапов.

Первый этап соответствует первым m рекурсивным шагам, на каждом из которых можно выделить следующие три последовательные макрооперации:

  • вычисление матриц S_{i}, T_{i}, \: i = \overline{1, 7} ;
  • рекурсивный вызов для вычисления произведений \alpha_{i} = S_{i} \cdot T_{i}, \: i = \overline{1, 7} , таким образом при переходе на следующий рекурсивный шаг задача размера k сводится к семи задачам размера \frac{k}{2} ;
  • вычисление матрицы C из матриц \alpha_{i} .

Второй этап соответствует последнему, m+1 -ому шагу алгоритма и выходу из рекурсии. Он представляет из себя вычисление произведения двух матриц матриц порядка n_{m} = \frac{n}{2^{m}} по классическому алгоритму:

c_{ij} = \sum_{k = 1}^{n_{m}} a_{ik} b_{kj}, \: i,j = \overline{1, n_{m}} .

1.5 Схема реализации последовательного алгоритма

Псевдокод простейшего варианта реализации последовательного алгоритма Штрассена, непосредственно вытекающего из его математического описания, может быть записан следующим образом:

	Strassen(A, B, C, k, m) ! C = A * B
           if (m == 0) call mul(A, B, C, k) ! переход к классическому алгоритму умножения матриц
           do i = 1, 7
               compute S[i], T[i] from A and B
           end do
           do i = 1, 7
               call Strassen(S[i], T[i], alpha[i], k/2, m-1)
           end do
           do i = 1, 7
               compute parts of C from alpha[i]
           end do
        end Strassen

Таким образом, на каждом рекурсивном шаге эта схема состоит из трёх последовательных этапов:

  • Вычисление семи пар матриц S_{i}, T_{i} . Так как они определены независимо друг от друга и хранятся в разных областях памяти, порядок их вычисления может быть произвольным, в том числе параллельным.
  • Выполнение семи рекурсивных вызовов функции Strassen, по одному для каждой тройки \alpha_{i} = S_{i} \cdot T_{i}. Порядок этих вызовов также может быть произвольным и параллельным.
  • Вычисление матрицы С из полученных после рекурсивных вызовов матриц \alpha_{i} .

Как видно, данная реализация требует большого количества дополнительной памяти, так как на каждом рекурсивном шаге необходимо одновременно хранить семь пар \frac{k}{2}\times\frac{k}{2}-матриц S_{i}, T_{i} и семь результатов их умножения \alpha_{i} .

Есть и другой вариант реализации последовательного алгоритма, позволяющий значительно сократить количество требуемой памяти:

	Strassen(A, B, C, k, m) ! C = A * B
           if (m == 0) call mul(A, B, C, k) ! переход к классическому алгоритму умножения матриц
           do i = 1, 7
               compute S[i], T[i] from A and B
               call Strassen(S[i], T[i], alpha[i], k/2, m-1)
               compute part of C from alpha[i]
           end do
        end Strassen

На каждом рекурсивном шаге схема этой реализации представляет собой цикл по i = 1,\dotsc,7, каждая итерация которого состоит из трёх последовательных этапов:

  • Вычисление i-ой пары матриц S_{i}, T_{i} .
  • Выполнение рекурсивного вызова функции Strassen для вычисления \alpha_{i} = S_{i} \cdot T_{i} .
  • Обновление блоков матрицы С в соответствии с формулами алгоритма в зависимости от номера полученной после рекурсивного вызова матрицы \alpha_{i} .

Порядок итераций такого цикла может быть произвольным, но они должны выполняться строго последовательно. Такое требование позволяет выделить память всего для трёх \frac{k}{2}\times\frac{k}{2}-матриц на каждом рекурсивном шаге: для одной пары S_{i}, T_{i} и их произведения \alpha_{i} , и на каждой итерации цикла использовать одни и те же области памяти для хранения соответствующей тройки матриц.

1.6 Последовательная сложность алгоритма

Оба вышеописанных варианта реализации метода Штрассена требуют одинакового числа арифметических операций. В соответствии с макроструктурой алгоритма, рассмотрим каждый из его двух этапов отдельно.

Первый этап требует на первом рекурсивном шаге выполнения 18 \left( \frac{n}{2} \right)^2 сложений (вычитаний). При переходе на следующий рекурсивный шаг, задача будет сведена к семи аналогичным задачам, но вдвое меньшей размерности. Отсюда получаем, что на k-ом шаге потребуется 18 \cdot 7^{k-1} \left( \frac{n}{2^{k}} \right)^2 операций. Тогда всего на первых m шагах потребуется следующее число сложений (вычитаний):

18 \left( \frac{n}{2} \right)^2 + 18 \cdot 7 \left( \frac{n}{4} \right)^2 + \dotsb + 18 \cdot 7^{m-1} \left( \frac{n}{2^{m}} \right)^2 = \frac{18}{4}n^{2}\left( 1 + \frac{7}{4} + \dotsb + \frac{7^{m-1}}{4^{m-1}}\right) = 6n^2\left(\frac{7^{m}}{4^{m}} - 1 \right) .

Получена сложность первого этапа алгоритма.

Второй этап включает в себя решение 7^m задач перемножения матриц порядка \frac{n}{2^{m}} по классическому алгоритму, что требует 7^m \left(\frac{n}{2^{m}}\right)^3 умножений и 7^m \left(\frac{n}{2^{m}}\right)^2\left(\frac{n}{2^{m}}-1\right) сложений, всего O\left(\frac{n^3}{2^{3m}}\right) арифметических операций.

Всего последовательный алгоритм Штрассена требует:

  • 6n^2\left(\frac{7^{m}}{4^{m}} - 1 \right) + 7^m\left(\frac{n}{2^{m}}\right)^2\left(\frac{n}{2^{m}}-1\right) сложений;
  • 7^m\left(\frac{n}{2^{m}}\right)^3 умножений.

Как видно, сложность алгоритма зависит не только от порядка матрицы n , но и от числа m рекурсивных шагов первого этапа.

Если m = L, где n = 2^L, то есть второй этап алгоритма сводится лишь к 7^L умножениям двух чисел, то сложность всего последовательного алгоритма равна

6n^2\left(\frac{7^{L}}{4^{L}} - 1 \right) + 7^L\left(\frac{n}{2^{L}}\right)^3 = 6n^2\left(\frac{7^{L}}{n^{2}} - 1 \right) + 7^L = 7 \cdot 7^{L} - 6n^2 = 7 n^{\log_2 7} - 6n^2 = O\left(n^{\log_2 7}\right) \approx O\left(n^{2.81}\right) .

Таким образом, получена оценка сложности последовательного алгоритма Штрассена с максимально возможным числом m = L рекурсивных шагов, соответствующая максимальной оценке сложности первого этапа алгоритма для произвольного m. Эта оценка асимптотически лучше оценки для классического метода перемножения матриц, требующего O(n^3) операций.

Стоит отметить роль правильного выбора параметра m . Из полученных оценок видно, что при малых m сложность алгоритма с уменьшением этого параметра будет стремиться к кубической сложности классического метода, однако при m , очень близких к L , алгоритм будет требовать большего числа сложений, чем классический метод, которое будет увеличиваться с ростом m . Отсюда следует, что существует некоторое оптимальное значение m \in [1, L] , которое минимизирует число арифметических операций всего алгоритма. Также при выборе m стоит учитывать роль накладных расходов на реализацию рекурсии при малых размерах матриц.

1.7 Информационный граф

1.8 Ресурс параллелизма алгоритма

1.9 Входные и выходные данные алгоритма

1.10 Свойства алгоритма

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

2.2 Локальность данных и вычислений

2.3 Возможные способы и особенности параллельной реализации алгоритма

2.4 Масштабируемость алгоритма и его реализации

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

2.7 Существующие реализации алгоритма

3 Литература

  1. Тыртышников Е.Е. Матричный анализ и линейная алгебра. М.: 2004-2005.