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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[math] 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] C = AB [/math], согласно алгоритму Штрассена, вычисляется следующим образом:

[math] \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] \alpha_{i} = S_{i} \cdot T_{i}, \: i = \overline{1,7}, [/math]

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

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

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

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

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

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

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

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

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

[math] r_{ij} = p_{ij} + q_{ij}, \: i,j = \overline{1, k} [/math].

В таком виде могут быть представлены все вычисления на каждом из первых [math] m [/math] рекурсивных шагов, а именно вычисление матриц [math] S_{i}, \: T_{i}, \: i = \overline{1,7} [/math] и блоков матрицы [math] C [/math] из матриц [math] \alpha_{i}, \: i = \overline{1,7} [/math].

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

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