Участник:DenisAnuprienko/Метод Штрассена

Основные авторы описания: Д.В.Ануприенко.

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

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

Метод Штрассена предназначен для умножения матриц. Здесь будет рассмотрен вариант метода, который можно применять к квадратным матрицам размера [math]\displaystyle{ N = 2^n }[/math]. В таком случае две матрицы можно умножить быстрее, чем за [math]\displaystyle{ O(N^3) }[/math].

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

Пусть имеются две матрицы [math]\displaystyle{ A, B \in \mathbb{R}^{N\times N} }[/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] При обычном умножении матриц пришлось бы совершить 8 умножений подматриц порядка [math]\displaystyle{ N/2 }[/math]. В методе Штрассена предлагается обойтись всего 7 умножениями. Находятся 7 вспомогательных подматриц [math]\displaystyle{ M_1, ..., M_7 }[/math] по следующим формулам:
[math]\displaystyle{ M_1 = (A_{11} + A_{22})(B_{11} + B_{22}) }[/math]
[math]\displaystyle{ M_2 = (A_{21} + A_{22})B_{11} }[/math]
[math]\displaystyle{ M_3 = A_{11}(B_{12} - B_{22}) }[/math]
[math]\displaystyle{ M_4 = A_{22}(B_{21} - B_{11}) }[/math]
[math]\displaystyle{ M_5 = (A_{11} + A_{12})B_{22} }[/math]
[math]\displaystyle{ M_6 = (A_{21} - A_{22})(B_{11} + B_{12}) }[/math]
[math]\displaystyle{ M_7 = (A_{12} - A_{22})(B_{21} + B_{22}) }[/math]

После этого матрица [math]\displaystyle{ C }[/math], являющаяся произведением [math]\displaystyle{ A }[/math] и [math]\displaystyle{ B }[/math], находится по формулам
[math]\displaystyle{ C_{11} = M_1 + M_4 - M_5 + M_7 }[/math]
[math]\displaystyle{ C_{12} = M_3 + M_5 }[/math]
[math]\displaystyle{ C_{21} = M_2 + M_4 }[/math]
[math]\displaystyle{ C_{22} = M_1 - M_2 + M_3 + M_6 }[/math]

Если и умножения подматриц, необходимые для нахождения [math]\displaystyle{ M_i }[/math], проводить по такой же схеме, получается рекурсивный алгоритм. Всего в нем понадобится выполнить [math]\displaystyle{ O(7^{log_2N}) = O(N^{log_27}) \approx O(N^{2.81}) }[/math] умножений. На практике рекурсию можно не разворачивать до конца, а использовать обычное умножение уже на матрицах размера 512.

[math]\displaystyle{ }[/math]

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

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

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

1. Если размер матриц меньше некоторого числа [math]\displaystyle{ N_{min} }[/math], умножить их обычным способом.
2. Иначе
   1. Сформировать множители для матрицы [math]\displaystyle{ M_1 }[/math]
   2. Применить метод Штрассена для этих множителей
   3. Сформировать множители для матрицы [math]\displaystyle{ M_2 }[/math]
   4. Применить метод Штрассена для этих множителей
   5. ...
   6. Сформировать множители для матрицы [math]\displaystyle{ M_7 }[/math]
   7. Применить метод Штрассена для этих множителей
   8. Сформировать результат из матриц [math]\displaystyle{ M_1, ..., M_7 }[/math].

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

Метод Штрассена позволяет сократить число умножений, поэтому оценивается именно оно. В методе Штрассена число умножений составляет [math]\displaystyle{ O(7^{log_2N}) = O(N^{log_27}) \approx O(N^{2.81}) }[/math].

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

Рассмотрим блок-схему для последовательной реализации рекурсивной функции Strassen для числа [math]\displaystyle{ N_{min} }[/math], равного 512.

Здесь 7 рекурсивных вызовов функции Strassen можно выполнять параллельно:

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

Умножения, которые необходимы для нахождения матриц [math]\displaystyle{ M_1, ..., M_7 }[/math], можно провести параллельно. Их можно предоставить 7 процессам (1 хозяин, который раздает задания 6 рабочим и получает результаты, а также работает вместе с рабочими) или 8 процессам (1 хозяин, который раздает задания 7 рабочим и получает результаты).

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

Нет никаких предположений насчет структуры матриц. Считается, что это обычные плотные матрицы. Они хранятся в виде одномерного массива, что позволяет легко выделять из них подматрицы.

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

• Алгоритм устойчив
• Алгоритм детерминирован

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

Исходный код

Исходный код

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

Матрицы хранятся как одномерные массивы. Благодаря этому в составе матрицы легко выделить подматрицу, зная размеры подматрицы и родительской матрицы, а также указатель на начало подматрицы.
В этом состоит выгодное отличие от некоторых существующих реализаций, где при использовании подматриц выделяется отдельная память.

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

Все действия с двумя главными матрицами проводит только процесс-хозяин. Каждый процесс также создает и освобождает для себя вспомогательные матрицы [math]\displaystyle{ M_i }[/math], а также две матрицы для записи множителей для нахождения [math]\displaystyle{ M_i }[/math].

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

Устройство метода Штрассена накладывает ограничения на количество процессов для распараллеливания.

• Как уже было отмечено, распараллеливать можно 7 умножений, из которых получаются матрицы [math]\displaystyle{ M_1, ..., M_7 }[/math], с помощью 7 или 8 процессов.
• Распараллеливание этих 7 умножений с помощью меньшего, чем 7, или большего, чем 8, числа процессов не рассматривается, так как в первом случае количество пересылок будет тем же, а время работы - большим, а во втором случае получается более 1 процесса на 1 умножение, что является дополнительным усложнением с негарантированной пользой.
• В таком случае, распараллеливание [math]\displaystyle{ n }[/math] уровней рекурсии требует как минимум [math]\displaystyle{ 7^n }[/math] процессов. Количество процессов меняется тогда и только тогда, когда меняется число распараллеливаемых уровней рекурсии.
• Здесь будет рассмотрен вариант метода Штрассена, где распараллеливается 1 уровень рекурсии с помощью 8 процессов: 1 хозяина и 7 рабочих. Эта версия работает несколько быстрее, чем версия с 7 процессами.

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

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

Сравнение работы последовательной и параллельной реализаций

Эксперименты проводились на суперкомпьютере "Ломоносов". Использовались gcc и mpicc 4.4.7 с флагом компиляции -O3, а также OpenMPI 1.8.4. Параллельная версия задействовала 8 процессов.

Результаты запусков на "Ломоносове"

Размер матрицы	Последовательная реализация	Параллельная реализация с использованием MPI	Ускорение
512	0.18 c	0.035 c	5.14
1024	1.99 c	0.40 c	4.98
2048	16.1 c	3.45 c	4.67
4096	230.1 c	71.0 c	3.23

Литература

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