Русская Википедия:Алгоритм умножения матриц
Шаблон:Unsolved Поскольку умножение матриц является центральной операцией во многих численных алгоритмах, много усилий было вложено в повышение эффективности алгоритма умножения матриц. Приложения алгоритма умножения матриц в вычислительных задачах найдены во многих областях, включая Шаблон:Не переведено 5 и распознавания образов, а также во вроде бы не имеющих отношение к матрицам задачах, таких как подсчёт путей через графШаблон:Sfn. Было разработано много различных алгоритмов для умножения матриц на оборудовании различного типа, включая параллельные и распределённые системы, где вычисления распределены на несколько процессоров (и, может быть, по сети).
Прямое использование математического определения умножения матриц даёт алгоритм, который работает за время порядка <math>n^3</math> операций поля для умножения двух <math>n \times n</math> матриц над этим полем (или <math>\Theta(n^3)</math> в нотации «O» большое). Улучшенные асимптотические границы по времени были известны с момента появления алгоритма Штрассена в 1960-х годах, но оптимальное время остаётся неизвестным (то есть, неизвестна сложность задачи). К концу 2020 года алгоритм умножения матриц с лучшей асимптотической сложностью, работающий за время <math>O(n^{2{,}3728596})</math>, был дан Джозефом Алманом и Шаблон:Не переведено 5Шаблон:Sfn, однако этот алгоритм галактического масштаба, то есть только для данных галактического размера, поскольку содержит огромные константы и не может быть реализован на практике.
Итеративный алгоритм
По определению умножения матриц для <math>n \times m</math> матрицы Шаблон:Mvar и по <math>m \times p</math> матрицы Шаблон:Mvar произведением <math>C = AB</math> является <math>n \times p</math> матрица, состоящая из элементов
- <math>c_{ij} = \sum_{k=1}^m a_{ik} b_{kj}</math>.
Отсюда можно построить простой алгоритм путём организации циклов по индексу Шаблон:Mvar от 1 до Шаблон:Mvar и Шаблон:Mvar от 1 до Шаблон:Mvar, и осуществляя вычисления по вышеприведённой формуле с помощью вложенных циклов:
- Вход: матрицы Шаблон:Mvar и Шаблон:Mvar
- Пусть Шаблон:Mvar будет новой матрицей нужного размера
- Для Шаблон:Mvar от 1 до Шаблон:Mvar:
- Для Шаблон:Mvar от 1 до Шаблон:Mvar:
- Положим <math>sum = 0</math>
- Для Шаблон:Mvar от 1 до Шаблон:Mvar:
- Положим <math>sum \leftarrow sum + A_{ik} \times B_{kj}</math>
- Положим <math>C_{ij} \leftarrow sum</math>
- Для Шаблон:Mvar от 1 до Шаблон:Mvar:
- Возвращаем Шаблон:Mvar
Этот алгоритм работает за время <math>\Theta(nmp)</math> (в асимптотических обозначениях) Шаблон:Sfn. Обычно для упрощения анализа алгоритма предполагается, что входными матрицами являются квадратные матрицы размера <math>n \times n</math>, и в этом случае время работы составляет <math>\Theta(n^3)</math>, то есть время зависит кубически от размера матрицШаблон:Sfn.
Поведение кэша
Три цикла при итеративном умножении матриц можно произвольным образом переставлять друг с другом без влияния на правильность алгоритма или асимптотическое время работы. Однако, порядок циклов может влиять на практические характеристики Шаблон:Не переведено 5 и на алгоритм использования кэшаШаблон:Sfn. Какой порядок вычисления лучше, зависит от того, как хранятся входные матрицы — в Шаблон:Не переведено 5, постолбцовом порядке, или в смешении этих порядков.
В частности, в идеальном случае полностью ассоциативного кэша, состоящего из Шаблон:Mvar байт и Шаблон:Mvar байт на строку кэша (то есть, Шаблон:Sfrac строк кэша), вышеприведённый алгоритм является подоптимальным для Шаблон:Mvar и Шаблон:Mvar, хранящихся построчно. Если <math>n > \tfrac{M}{b}</math>, любая итерация внутреннего цикла (одновременный проход по строкам Шаблон:Mvar и столбцам Шаблон:Mvar) приводит к промахам кэша при выборке элемента матрицы Шаблон:Mvar. Это означает, что при работе алгоритма будет <math>\Theta(n^3)</math> промахов кэша в Шаблон:Не переведено 5. К 2010 году для больших матриц скорость доступа к памяти являлась доминирующим фактором, определяющим время работы алгоритма, а не скорость процессора Шаблон:Sfn.
Оптимальным вариантом итерационного алгоритма для матриц Шаблон:Mvar и Шаблон:Mvar при построчном хранении является версия с Шаблон:Не переведено 5, где матрицы в неявном виде разбиты на квадратные части размера <math>\sqrt{M} \times \sqrt{M}</math>Шаблон:SfnШаблон:Sfn:
- Вход: матрицы Шаблон:Mvar и Шаблон:Mvar
- Пусть Шаблон:Mvar будет новой матрицей нужного размера
- Выбираем часть размера <math>T = \Theta(\sqrt{M})</math>
- Для Шаблон:Mvar от 1 до Шаблон:Mvar для части Шаблон:Mvar:
- Для Шаблон:Mvar от 1 до Шаблон:Mvar для части Шаблон:Mvar:
- Для Шаблон:Mvar от 1 до Шаблон:Mvar для части Шаблон:Mvar:
- Умножаем <math>A_{I:I+T, K:K+T}</math> и <math>B_{K:K+T, J:J+T}</math>, помещая результат в <math>C_{I:I+T, J:J+T}</math>, то есть:
- Для Шаблон:Mvar от Шаблон:Mvar до <math>\min(I + T, n)</math>:
- Для Шаблон:Mvar от Шаблон:Mvar до <math>\min(J + T, p)</math>:
- Положим Шаблон:Math
- Для Шаблон:Mvar от Шаблон:Mvar до <math>\min(K + T, m)</math>:
- Положим <math>sum \leftarrow sum + A_{ik} \times B_{kj}</math>
- Положим <math>C_{ij} \leftarrow C_{ij} + sum</math>
- Для Шаблон:Mvar от Шаблон:Mvar до <math>\min(J + T, p)</math>:
- Для Шаблон:Mvar от 1 до Шаблон:Mvar для части Шаблон:Mvar:
- Для Шаблон:Mvar от 1 до Шаблон:Mvar для части Шаблон:Mvar:
- Возвращаем Шаблон:Mvar
В случае идеального кэша алгоритм приводит только к <math>\Theta(\tfrac{n^3}{b \sqrt{M}})</math> промахам. Делитель <math>b\sqrt{M}</math> составляет несколько порядков величины для современных машин, так что вычисления доминируют во времени работы, а не промахи кэшаШаблон:Sfn.
Алгоритм Разделяй-и-властвуй
Альтернативой итерационному алгоритму для умножения матриц является алгоритм разделяй-и-властвуй. Он опирается на разложение на блоки
- <math>C = \begin{pmatrix}
C_{11} & C_{12} \\ C_{21} & C_{22} \\ \end{pmatrix},\, A = \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \\ \end{pmatrix},\, B = \begin{pmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \\ \end{pmatrix}</math>,
что работает для всех квадратных матриц с размерностями, равными степеням двойки, то есть <math>2^n\times 2^n</math> для некоторого Шаблон:Mvar. Произведение матриц тогда равно
- <math>\begin{pmatrix}
C_{11} & C_{12} \\ C_{21} & C_{22} \\ \end{pmatrix} = \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \\ \end{pmatrix} \begin{pmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \\ \end{pmatrix} = \begin{pmatrix} A_{11} B_{11} + A_{12} B_{21} & A_{11} B_{12} + A_{12} B_{22}\\ A_{21} B_{11} + A_{22} B_{21} & A_{21} B_{12} + A_{22} B_{22}\\ \end{pmatrix} </math>
что составляет восемь умножений пар подматриц с последующим шагом сложения. Алгоритм разделяй-и-властвуй вычисляет элементы рекурсивно с помощью скалярного произведения <math>c_{11}=a_{11}b_{11}</math> как в базовом случае.
Сложность этого алгоритма как функция от Шаблон:Mvar задаётся рекурсиейШаблон:Sfn
- <math>T(1) = \Theta(1)</math>;
- <math>T(n) = 8T(n/2) + \Theta(n^2)</math>,
принимающей во внимание восемь рекурсивных вызовов на матрицах размера Шаблон:Math и величину <math>\Theta(n^2)</math> для суммирования четырёх пар полученных матриц. Применяя основную теорему о рекуррентных соотношениях, получим, что эта рекурсия имеет решение <math>\Theta(n^3)</math>, ту же самую сложность, что и для итеративного алгоритмаШаблон:Sfn.
Неквадратные матрицы
Вариант этого алгоритма, который работает для матриц произвольного размера и на практике быстрееШаблон:Sfn, разбивает матрицы на две, а не на четыре подматрицы, что продемонстрировано нижеШаблон:Sfn. Разбиение матрицы теперь означает разделение её на две одинаковые или близкие к одинаковым части, если размер нечётен.
- Вход: матрица Шаблон:Mvar размера <math>n \times m</math> и матрица Шаблон:Mvar размера <math>m \times p</math>.
- Базовый случай: если <math>\max(n, m, p)</math> ниже некоторого порога, используем версию с размоткой итеративного алгоритма.
- Случай рекурсии:
- Если <math>\max(n, m, p) = n</math>, разбиваем матрицу Шаблон:Mvar горизонтально:
- <math>C = \begin{pmatrix} A_1 \\ A_2 \end{pmatrix} {B}
= \begin{pmatrix} A_1 B \\ A_2 B \end{pmatrix}</math>
- В противном случае, если <math>\max(n, m, p) = p</math>, разбиваем матрицу Шаблон:Mvar вертикально:
- <math>C = A \begin{pmatrix} B_1 & B_2 \end{pmatrix}
= \begin{pmatrix} A B_1 & A B_2 \end{pmatrix} </math>
- В противном случае <math>\max(n, m, p) = m</math>. Разбиваем Шаблон:Mvar вертикально, а Шаблон:Mvar горизонтально:
- <math>C = \begin{pmatrix} A_1 & A_2 \end{pmatrix} \begin{pmatrix} B_1 \\ B_2 \end{pmatrix}
= A_1 B_1 + A_2 B_2</math>
Поведение кэша
Количество промахов кэша рекрсивного умножения матриц та же самая, что и у версии aлгоритма Шаблон:Не переведено 5, но, в отличие от этого алгоритма, рекурсивный алгоритм Шаблон:Не переведено 5Шаблон:Sfn — не нужно никакого настоечного параметра для получения оптимального поведения кэша и он работает хорошо в многозадачном окружении, где размер кэша меняется динамически, поскольку другие процессы тоже используют кэшШаблон:Sfn. (Простой итеративный алгоритм нечувствителен к кэшированию также, но на практике много медленнее, если матрицы не подогнаны под алгоритм.)
Число промахов кэша при этом алгоритме на машинах с Шаблон:Mvar линиями идеального кэша, каждая из которых имеет размер в Шаблон:Mvar байт, ограничена величинойШаблон:Sfn
- <math>\Theta \left(m + n + p + \frac{mn + np + mp}{b} + \frac{mnp}{b\sqrt{M}} \right)</math>
Подкубические алгоритмы
Существуют алгоритмы, обеспечивающие лучшее время работы, чем прямолинейные алгоритмы. Первым из таких алгоритмов был открыт алгоритм Штрассена, разработанный Фолькером Штрассеном в 1969 и часто упоминаемый как «быстрое умножение матриц». Алгоритм основан на способе перемножения двух <math>2 \times 2</math> матриц, который требует лишь 7 умножений (вместо обычных 8), но требует выполнения дополнительных операций сложения и вычитания. Применение такого подхода рекурсивно даёт алгоритм с ценой по умножению <math>O( n^{\log_{2}7}) \approx O(n^{2,807})</math>. Алгоритм Штрассена более сложен, а вычислительная устойчивость хуже, чем у наивного алгоритмаШаблон:Sfn, но он быстрее в случае, когда <math>n > 100</math> или где-то около этогоШаблон:Sfn, и алгоритм включён в некоторые библиотеки, такие как BLASШаблон:Sfn. Алгоритм очень полезен для больших матриц над точными областями, такими как конечные поля, где вычислительная устойчивость не играет роли.
Открытым вопросом теоретической информатики является вопрос, насколько можно улучшить алгоритм Штрассена. Показатель умножения матриц <math>\omega</math> — это наименьшее вещественно число, для которого произведение любых двух <math>n\times n</math> матриц над полем может быть вычислено за <math>n^{\omega + o(1)}</math> операций в поле. Текущая лучшая граница <math>\omega</math> равна 2,3728596 (алгоритм Джошуа Алманаи Шаблон:Не переведено 5Шаблон:Sfn. Этот алгоритм, подобно всем другим недавно разработанным алгоритмам в этом направлении исследований, является обобщением алгоритма Копперсмита — Винограда, который представили Дон Копперсмит и Шаблон:Не переведено 5 в 1990 году, и который имеет асимптотическую сложность <math>O(n^{2,376})</math>. Концептуальная идея этих алгоритмов аналогична алгоритму Штрассена — способ разрабатывается для умножения двух <math>k \times k</math> матриц за менее чем <math>k^3</math> умножений и эта техника применяется рекурсивно. Однако, константы, спрятанные в нотации «O большое» так велики, что эти алгоритмы целесообразно применять только для матриц, которые слишком велики, чтобы их можно было обрабатывать на существующих компьютерахШаблон:SfnШаблон:Sfn.
Поскольку любой алгоритм умножения двух <math>n \times n</math> матриц должен обработать все <math>2n^2</math> элемента, имеется асимптотическая нижняя граница числа операций <math>\Omega(n^2)</math>. Рэн Раз доказал нижнюю границу в <math>\Omega(n^2 log(n))</math> ограниченных коэффициентов арифметических цепей над вещественными или комплексными числамиШаблон:Sfn.
Кон с соавторами переложил методы, такие как алгоритмы Штрассена и Копперсмита — Винограда в полностью другой контекст теории групп путём использования троек подмножеств конечных групп, которые удовлетворяют свойством дизъютктности, называемое Шаблон:Не переведено 5 (СТП, Шаблон:Lang-en, TPP). Они показали, что если семейства Шаблон:Не переведено 5 абелевых групп с симметричными группами образуют семейства троек со свойствами, аналогичными СТП, то существуют алгоритмы умножения матриц с фактически квадратичной сложностьюШаблон:SfnШаблон:Sfn. Большинство исследователей верят, что это верно вообщеШаблон:Sfn. Однако Алон, Шпилька и Хрис Уманс недавно показали, что некоторые из гипотез о быстром умножении матриц несовместимы с другой правдоподобной гипотезой, Шаблон:Не переведено 5[1].
Алгоритм Фрейвалдса — это простой алгоритм Монте-Карло, который для заданных матриц <math>A, B</math> и Шаблон:Mvar проверяет, что <math>AB = C</math>, за время <math>\Theta(n^2)</math>.
Параллельные и распределённые алгоритмы
Параллельность с разделением памяти
Алгоритм разделяй-и-властвуй, описанный выше, может быть распараллелен двумя способами для мультипроцессоров с разделяемой памятью. Это основывается на факте, что восемь алгоритмов рекурсивного умножения матриц в
- <math>\begin{pmatrix}
A_{11} B_{11} + A_{12} B_{21} & A_{11} B_{12} + A_{12} B_{22}\\ A_{21} B_{11} + A_{22} B_{21} & A_{21} B_{12} + A_{22} B_{22}\\ \end{pmatrix} </math>
можно осуществлять независимо друг от друга, как и сложение (хотя алгоритм тербует «объединения» умножений перед осуществлением сложения). Воплощая полный параллелизм задачи, получаем алгоритм, который можно выразить в стиле fork–join псевдокодаШаблон:Sfn:
Шаблон:Начало коробки Процедура Шаблон:Math:
- Базовый случай: если <math>n = 1</math>, положим <math>c_{11} \leftarrow a_{11} \times b_{11}</math> (умножаем маленькие блочные матрицы).
- В противном случае распределяем место для новой матрицы Шаблон:Mvar размером <math>n \times n</math>, затем:
- Разбиваем Шаблон:Mvar на <math>A_{11}, A_{12}, A_{21}, A_{22}</math>.
- Разбиваем Шаблон:Mvar на <math>B_{11}, B_{12}, B_{21}, B_{22}</math>.
- Разбиваем Шаблон:Mvar на <math>C_{11}, C_{12}, C_{21}, C_{22}</math>.
- Разбиваем Шаблон:Mvar на <math>T_{11}, T_{12}, T_{21}, T_{22}</math>.
- Распараллеливаем (Fork = вилка, то есть ответвляем процесс):
- Fork <math>multiply(C_{11}, A_{11}, B_{11})</math>.
- Fork <math>multiply(C_{12}, A_{11}, B_{12})</math>.
- Fork <math>multiply(C_{21}, A_{21}, B_{11})</math>.
- Fork <math>multiply(C_{22}, A_{21}, B_{12})</math>.
- Fork <math>multiply(T_{11}, A_{12}, B_{21})</math>.
- Fork <math>multiply(T_{12}, A_{12}, B_{22})</math>.
- Fork <math>multiply(T_{21}, A_{22}, B_{21})</math>.
- Fork <math>multiply(T_{22}, A_{22}, B_{22})</math>.
- Join (Join = объединение, ждём завершения разветвлённых процессов).
- Шаблон:Math.
- Уничтожаем Шаблон:Mvar.
Процедура <math>add(C, T)</math> добавляет Шаблон:Mvar к Шаблон:Mvar поэлементно:
- Базовый случай: если <math>n = 1</math>, полагаем <math>c_{11} \leftarrow c_{11} + t_{11}</math> (или делаем короткий цикл, возможно, с размоткой).
- В противном случае:
- Разбиваем Шаблон:Mvar на Шаблон:Math, Шаблон:Math, Шаблон:Math, Шаблон:Math.
- Разбиваем Шаблон:Mvar на Шаблон:Math, Шаблон:Math, Шаблон:Math, Шаблон:Math.
- Распараллеливаем:
- Fork <math>add(C_{11}, T_{11})</math>.
- Fork <math>add(C_{12}, T_{12})</math>.
- Fork <math>add(C_{21}, T_{21})</math>.
- Fork <math>add(C_{22}, T_{22})</math>.
- Join.
Здесь, fork означает, что вычисления могут осуществляться параллельно к остальной части процедуры, а join означает ожидание завершения всех запущенных в параллельные ветки вычислений. Распараллеливание достигает своей цели лишь передачей указателей.
Этот алгоритм имеет Шаблон:Не переведено 5 в <math>\Theta(\log^2n)</math> шагов, что определяет требуемое время для идеальной машины с неограниченным числом процессоров. Поэтому алгоритм имеет максимальное возможное Шаблон:Не переведено 5 <math>\Theta(n^3/\log^2n)</math> на любом реальном компьютере. Алгоритм не имеет практического значения ввиду неустранимой цены передачи данных во временную матрицу Шаблон:Mvar и из неё, но более практичный вариант, не использующий временных матриц, достигает ускорения <math>\Theta(n^2)</math>Шаблон:Sfn.
Алгоритмы без обмена данных и распределённые алгоритмы
В современных архитектурах с иерархической памятью цена загрузки и выгрузки входной матрицы стремится к определяющей роли при обработке. На одной машине это количество данных, переносимых между RAM и кэшем, в то время как для распределённой памяти машин с несколькими узлами это величина, переносимая между узлами. В любом случае это называется полосой пропускания. Наивный алгоритм, где используются три вложенных цикла, использует полосу <math>\Omega(n^3)</math>.
Шаблон:Не переведено 5, известный также как алгоритм 2D, — это Шаблон:Не переведено 5, который превращает каждую входную матрицу в блочную матрицу, элементами которой являются подматрицы размера <math>\sqrt{\tfrac{M}{3}} \times \sqrt{\tfrac{M}{3}}</math>, где Шаблон:Math является размером быстрой памяти Шаблон:Sfn. Затем используется наивный алгоритм над блоками матриц, вычисляющий произведение подматриц полностью в быстрой памяти. Это сокращает полосу частот канала связи до <math>O(\tfrac{n^3}{\sqrt{M}})</math>, что асимптотически оптимально (для алгоритмов, выполняющих <math>\Omega(n^3)</math> операций)Шаблон:SfnШаблон:Sfn.
В распределённых вычислениях с Шаблон:Mvar процессорами, организованными в двухмерную решётку <math>\sqrt{p} \times \sqrt{p}</math>, каждая из результирующих подматриц может быть назначена процессору и произведение может быть вычислено каждым процессором с передачей <math>O(\tfrac{n^2}{\sqrt{p}})</math> слов, что асимптотически оптимально в предположении, что каждый узел сохраняет минимум <math>O(\tfrac{n^2}{p})</math> элементовШаблон:Sfn. Это может быть улучшено алгоритмом 3D, который распределяет процессоры в трёхмерную кубическую решётку, путём назначения каждого произведения двух входных подматриц одному процессору. Полученные подматрицы затем генерируются путём работы над каждой строкойШаблон:Sfn. Этот алгоритм передаёт <math>O(\tfrac{n^2}{p^{2/3}})</math> слов на процессор, что асимптотически оптимальноШаблон:Sfn. Однако, это требует репликации каждого элемента входной матрицы <math>p^{\tfrac{1}{3}}</math> раз, а потому требует в <math>p^{\tfrac{1}{3}}</math> больше памяти, чем нужно для хранения входных данных. Этот алгоритм можно комбинировать с алгоримом Штрассена для дальнейшего уменьшения времени работыШаблон:Sfn. «2,5D» алгоритмы обеспечивают постоянный обмен между использованием памяти и полосой частот обменаШаблон:Sfn. В современных системах распределённых вычислений, таких как MapReduce, были разработаны специализированные алгоритмы умноженияШаблон:Sfn.
Алгоритмы для ячеистых топологий
Имеется ряд алгоритмов вычисления умножения в ячеистой топологии. Для умножения двух <math>n \times n</math> матриц на стандартной двумерной ячеистой топологии с помощью Шаблон:Не переведено 5 2D можно выполнить умножение за 3n-2 шагов, хотя это число сокращается вдвое для повторных вычисленийШаблон:Sfn. Стандартный массив неэффективен, поскольку данные из двух матриц не приходят одновременно и должны быть дополнены нулевыми значениями.
Результат даже быстрее на двухуровневой решётке с перекрёстными связями, где нужно только 2n-1 шаговШаблон:Sfn. Производительность улучшается далее для повторных вычислений, что приводит к 100% эффективностиШаблон:R. Решётка с перекрёстными связями может рассматриваться как специальный случай неплоской (то есть многослойной) вычислительной структурыШаблон:Sfn.
См. также
- Вычислительная сложность математических операций
- Алгоритм Валианта
- Задача о порядке перемножения матриц
Примечания
Литература
- Шаблон:Статья
- Шаблон:Книга
- Шаблон:Статья
- Шаблон:Книга
- Шаблон:Книга
- Шаблон:Статья
- Шаблон:Книга
- Шаблон:Книга
- Шаблон:Книга
- Шаблон:Книга
- Шаблон:Книга
- Шаблон:Статья Выдержка: «Алгоритм Копперсмита — Виноград непрактичен, поскольку содержит очень большие спрятанные константы в верхней границе числа требуемых умножений.
- Шаблон:Статья
- Шаблон:Книга
- Шаблон:Статья
- Шаблон:Статья
- Шаблон:Статья
- Шаблон:Статья
- Шаблон:Статья
- Шаблон:Статья
- Шаблон:Книга
- Шаблон:Статья
- Шаблон:Статья
Литература для дальнейшего чтения
Шаблон:Численные методы линейной алгебры
- ↑ Alon, Shpilka, Umans, On Sunflowers and Matrix Multiplication Шаблон:Wayback
