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

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ \tilde{A} = \begin{bmatrix} A & 0 \\ 0 & 0 \end{bmatrix}, \: \: \tilde{B} = \begin{bmatrix} B & 0 \\ 0 & 0 \end{bmatrix}, }[/math]

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

[math]\displaystyle{ \tilde{C} = \tilde{A}\tilde{B} = \begin{bmatrix} AB & 0 \\ 0 & 0 \end{bmatrix} = \begin{bmatrix} C & 0 \\ 0 & 0 \end{bmatrix}. }[/math]

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

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

[math]\displaystyle{ 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}. }[/math]

Их произведение [math]\displaystyle{ C = AB }[/math], согласно алгоритму Штрассена, вычисляется следующим образом:

[math]\displaystyle{ \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*} }[/math]

[math]\displaystyle{ \alpha_{i} = S_{i} \cdot T_{i}, \: i = \overline{1,7}, }[/math]

[math]\displaystyle{ C_{11} = \alpha_{1} + \alpha_{4} - \alpha_{5} + \alpha_{7}, }[/math]

[math]\displaystyle{ C_{12} = \alpha_{3} + \alpha_{5}, }[/math]

[math]\displaystyle{ C_{21} = \alpha_{2} + \alpha_{4}, }[/math]

[math]\displaystyle{ C_{22} = \alpha_{1} + \alpha_{3} - \alpha_{2} + \alpha_{6}, }[/math]

[math]\displaystyle{ C = \begin{bmatrix} C_{11} & C_{12} \\ C_{21} & C_{22} \end{bmatrix}. }[/math]

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

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

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

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

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

Алгоритм Штрассена состоит из [math]\displaystyle{ m+1 }[/math] рекурсивных шагов и может быть представлен в виде двух этапов.

Первый этап соответствует первым [math]\displaystyle{ m }[/math] рекурсивным шагам, на каждом из которых можно выделить следующие три последовательные макрооперации: 1. вычисление матриц [math]\displaystyle{ S_{i}, T_{i}, \: i = \overline{1, 7} }[/math]; 2. рекурсивный вызов для вычисления произведений [math]\displaystyle{ \alpha_{i} = S_{i} \cdot T_{i}, \: i = \overline{1, 7} }[/math], таким образом при переходе на следующий рекурсивный шаг задача размера [math]\displaystyle{ k }[/math] сводится к семи задачам размера [math]\displaystyle{ \frac{k}{2} }[/math]; 3. вычисление матрицы [math]\displaystyle{ C }[/math] из матриц [math]\displaystyle{ \alpha_{i} }[/math].

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

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

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

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

        call Strassen(a, b, c, n, m)

	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

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

1. Вычисление семи пар матриц [math]\displaystyle{ S_{i}, T_{i} }[/math]. Так как они определены независимо друг от друга, порядок их вычисления может быть произвольным, в том числе параллельным.

2. Выполнение семи рекурсивных вызовов функции Strassen -- по одному для каждой тройки [math]\displaystyle{ \alpha_{i} = S_{i} \cdot T_{i} }[/math]. Порядок этих вызовов также может быть произвольным и параллельным.

3. Вычисление матрицы [math]\displaystyle{ С }[/math] из полученных после рекурсивных вызовов матриц [math]\displaystyle{ \alpha_{i} }[/math].

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

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

        call Strassen(a, b, c, n, m)

	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

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

1. Вычисление [math]\displaystyle{ i }[/math]-ой пары матриц [math]\displaystyle{ S_{i}, T_{i} }[/math].

2. Выполнение рекурсивного вызова функции Strassen для вычисления [math]\displaystyle{ \alpha_{i} = S_{i} \cdot T_{i} }[/math].

3. Изменение значений блоков матрицы [math]\displaystyle{ С }[/math] в соответствии с формулами алгоритма в зависимости от номера полученной после рекурсивного вызова матрицы [math]\displaystyle{ \alpha_{i} }[/math].

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