Алгоритм штрассена для умножения матриц: Divide&Conquer над алгоритмом Штрассена / Хабр

{2,376})}

Содержание

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

Или же R унитарное кольцо и и B квадратные матрицы размера п , коэффициенты которого являются элементами набора , лежащим в основе R . Для простоты мы рассматриваем случай, когда n является степенью двойки. Мы всегда можем свести себя к этому случаю, возможно добавив столбцы и строки нулей. Алгоритм Штрассен может вычислить матричный продукт C .

Три матрицы A , B и C разделены на матрицы блоками одинакового размера:

Взнак равно[В1,1В1,2В2,1В2,2] , Bзнак равно[B1,1B1,2B2,1B2,2] , ПРОТИВзнак равно[ПРОТИВ1,1ПРОТИВ1,2ПРОТИВ2,1ПРОТИВ2,2]{\ displaystyle \ mathbf {A} = {\ begin {bmatrix} \ mathbf {A} _ {1,1} & \ mathbf {A} _ {1,2} \\\ mathbf {A} _ {2,1 } & \ mathbf {A} _ {2,2} \ end {bmatrix}} {\ mbox {,}} \ mathbf {B} = {\ begin {bmatrix} \ mathbf {B} _ {1,1} & \ mathbf {B} _ {1,2} \\\ mathbf {B} _ {2,1} & \ mathbf {B} _ {2,2} \ end {bmatrix}} {\ mbox {,}} \ mathbf {C} = {\ begin {bmatrix} \ mathbf {C} _ {1,1} & \ mathbf {C} _ {1,2} \\\ mathbf {C} _ {2,1} & \ mathbf {C} _ {2,2} \ end {bmatrix}}}

или же

Вя,j,Bя,j,ПРОТИВя,j∈Fнет/2×Fнет/2.

{n / 2}.}

Тогда у нас есть

ПРОТИВ1,1знак равноВ1,1B1,1+В1,2B2,1{\ displaystyle \ mathbf {C} _ {1,1} = \ mathbf {A} _ {1,1} \ mathbf {B} _ {1,1} + \ mathbf {A} _ {1,2} \ mathbf {B} _ {2,1}}

ПРОТИВ1,2знак равноВ1,1B1,2+В1,2B2,2{\ displaystyle \ mathbf {C} _ {1,2} = \ mathbf {A} _ {1,1} \ mathbf {B} _ {1,2} + \ mathbf {A} _ {1,2} \ mathbf {B} _ {2,2}}

ПРОТИВ2,1знак равноВ2,1B1,1+В2,2B2,1{\ displaystyle \ mathbf {C} _ {2,1} = \ mathbf {A} _ {2,1} \ mathbf {B} _ {1,1} + \ mathbf {A} _ {2,2} \ mathbf {B} _ {2,1}}

ПРОТИВ2,2знак равноВ2,1B1,2+В2,2B2,2{\ displaystyle \ mathbf {C} _ {2,2} = \ mathbf {A} _ {2,1} \ mathbf {B} _ {1,2} + \ mathbf {A} _ {2,2} \ mathbf {B} _ {2,2}}

Этот метод требует 8 умножений матриц для вычисления C _{i, j} , как в классическом произведении.

Сила алгоритма Штрассена заключается в наборе из семи новых матриц M _i, которые будут использоваться для выражения C _{i, j} всего с 7 умножениями вместо 8:

M1знак равно(В1,1+В2,2)(B1,1+B2,2){\ displaystyle \ mathbf {M} _ {1}: = (\ mathbf {A} _ {1,1} + \ mathbf {A} _ {2,2}) (\ mathbf {B} _ {1,1 } + \ mathbf {B} _ {2,2})}

M2знак равно(В2,1+В2,2)B1,1{\ Displaystyle \ mathbf {M} _ {2}: = (\ mathbf {A} _ {2,1} + \ mathbf {A} _ {2,2}) \ mathbf {B} _ {1,1} }

M3знак равноВ1,1(B1,2-B2,2){\ displaystyle \ mathbf {M} _ {3}: = \ mathbf {A} _ {1,1} (\ mathbf {B} _ {1,2} — \ mathbf {B} _ {2,2}) }

M4знак равноВ2,2(B2,1-B1,1){\ displaystyle \ mathbf {M} _ {4}: = \ mathbf {A} _ {2,2} (\ mathbf {B} _ {2,1} — \ mathbf {B} _ {1,1}) }

M5знак равно(В1,1+В1,2)B2,2{\ Displaystyle \ mathbf {M} _ {5}: = (\ mathbf {A} _ {1,1} + \ mathbf {A} _ {1,2}) \ mathbf {B} _ {2,2} }

M6знак равно(В2,1-В1,1)(B1,1+B1,2){\ displaystyle \ mathbf {M} _ {6}: = (\ mathbf {A} _ {2,1} — \ mathbf {A} _ {1,1}) (\ mathbf {B} _ {1,1 } + \ mathbf {B} _ {1,2})}

M7знак равно(В1,2-В2,2)(B2,1+B2,2){\ Displaystyle \ mathbf {M} _ {7}: = (\ mathbf {A} _ {1,2} — \ mathbf {A} _ {2,2}) (\ mathbf {B} _ {2,1 } + \ mathbf {B} _ {2,2})}

Тогда C _{i, j} выражаются как

ПРОТИВ1,1знак равноM1+M4-M5+M7{\ displaystyle \ mathbf {C} _ {1,1} = \ mathbf {M} _ {1} + \ mathbf {M} _ {4} — \ mathbf {M} _ {5} + \ mathbf {M} _ {7}}

ПРОТИВ1,2знак равноM3+M5{\ Displaystyle \ mathbf {C} _ {1,2} = \ mathbf {M} _ {3} + \ mathbf {M} _ {5}}

ПРОТИВ2,1знак равноM2+M4{\ Displaystyle \ mathbf {C} _ {2,1} = \ mathbf {M} _ {2} + \ mathbf {M} _ {4}}

ПРОТИВ2,2знак равноM1-M2+M3+M6{\ displaystyle \ mathbf {C} _ {2,2} = \ mathbf {M} _ {1} — \ mathbf {M} _ {2} + \ mathbf {M} _ {3} + \ mathbf {M} _ {6}}

Процесс повторяется рекурсивно до тех пор, пока матрицы A и B не станут «маленькими». (Де) Фолькер Штрассен, « Исключение Гаусса не оптимально » , Numerische Mathematik , vol. 13,1969 г., стр. 354-356.

↑ (in) Томас Х. Кормен , Чарльз Э. Лейзерсон , Рональд Л. Ривест и Клиффорд Стейн , Введение в алгоритмы , MIT Press ,2009 г., 3 ^е изд. [ деталь издания ], примечания гл. 4, стр.

112 .

Библиография

(ru) Томас Х. Кормен , Чарльз Э. Лейзерсон , Рональд Л. Ривест и Клиффорд Штайн , Введение в алгоритмы , MIT Press и McGraw-Hill ,2001 г., 2- ^е изд. [ деталь издания ], гл. 31, раздел 31.2 (« Алгоритм Штрассена для умножения матриц »),

с. 739-745 .

Внешние ссылки

(en) Эрик В. Вайстейн , « Формулы Штрассена » , на сайте MathWorld.
(ru) Анализ алгоритма Штрассена с использованием Maple
Жерар Сукахет, алгоритм В. Штрассена для быстрого умножения матриц

Умножение

Почтальон
- Множимое
- Множитель
Продукт
Крест умножения
Таблица умножения

Характеристики

Распределительность
Ассоциативность
Коммутативность

Примеры

Скалярное произведение
Векторный продукт
Умножение на скаляр
Матричный продукт

Алгоритмы умножения

Умножение целых чисел	Древний Египет русский Древний Китай Скольжением Ревностью Метод Трахтенберга Алгоритм умножения Бута Карацуба Тоом-Кук Шёнхаге-Штрассен Фюрер
Умножение матриц	Адамар Кронекер Стразы Медник-Виноград Цепной алгоритм умножения матриц

Алгоритмы верификации

Умножение целых чисел	Доказательство девяти
Умножение матриц	Алгоритм проверки Фрейвальдса

Алгоритм умножения матриц

Нерешенная проблема в информатике:

Какой алгоритм умножения матриц самый быстрый?

(больше нерешенных проблем в информатике)

Потому что матричное умножение такая центральная операция во многих численные алгоритмы, много работы было вложено в создание алгоритмы матричного умножения эффективный. Применение матричного умножения в вычислительных задачах можно найти во многих областях, включая научные вычисления и распознавание образов и в, казалось бы, не связанных между собой задачах, таких как подсчет путей через график.^[1] Для умножения матриц на различных типах оборудования было разработано множество различных алгоритмов, включая параллельно и распределен системы, в которых вычислительная работа распределяется между несколькими процессорами (возможно, по сети).

Непосредственное применение математического определения умножения матриц дает алгоритм, который занимает время в порядке п³ умножить два п × п матрицы (Θ (п³) в нотация большой O ). Лучшие асимптотические оценки времени, необходимого для умножения матриц, были известны со времен работы Штрассена в 1960-х годах, но до сих пор неизвестно, каково оптимальное время (т.е. сложность проблемы является).

Содержание

1 Итерационный алгоритм
- 1.1 Поведение кеша
2 Алгоритм разделяй и властвуй
- 2.1 Неквадратные матрицы
- 2.2 Поведение кеша
3 Субкубические алгоритмы
4 Параллельные и распределенные алгоритмы
- 4.1 Параллелизм с общей памятью
- 4.2 Распределенные алгоритмы предотвращения общения
- 4.3 Алгоритмы для сеток
5 Смотрите также
6 использованная литература
7 дальнейшее чтение

Итерационный алгоритм

В определение умножения матриц это если C = AB для п × м матрица А и м × п матрица B, тогда C является п × п матрица с записями

cяj=∑k=1маяkбkj{displaystyle c_ {ij} = sum _ {k = 1} ^ {m} a_ {ik} b_ {kj}}.

Исходя из этого, можно построить простой алгоритм, который перебирает индексы я с 1 по п и j с 1 по п, вычисляя указанное выше, используя вложенный цикл:

Вход: матрицы А и B
Позволять C быть новой матрицей подходящего размера
Для я от 1 до п:
- Для j от 1 до п:
  - Позволять сумма = 0
  - Для k от 1 до м:
    - Набор сумма ← сумма + А_ik × B_кДж
  - Набор C_ij ← сумма
Вернуть C

Этот алгоритм занимает время Θ (nmp) (в асимптотическая запись ).^[1] Обычное упрощение с целью анализ алгоритмов предполагает, что все входные данные представляют собой квадратные матрицы размера п × п, в этом случае время работы Θ (п³), т.е. кубический по размеру измерения.^[2]

Поведение кеша

Иллюстрация порядка строк и столбцов

Три цикла в итеративном умножении матриц можно произвольно менять местами без влияния на правильность или асимптотическое время выполнения. Однако порядок может существенно повлиять на практические характеристики из-за шаблоны доступа к памяти и тайник использование алгоритма;^[1]какой порядок лучше, также зависит от того, хранятся ли матрицы в рядовой порядок, порядок столбцов или их сочетание.

В частности, в идеализированном случае полностью ассоциативный кеш состоящий из M байты и б байтов на строку кэша (т.е. M/б строки кэша), приведенный выше алгоритм не является оптимальным для А и B хранятся в строчном порядке. Когда п > M/б, каждая итерация внутреннего цикла (одновременное прохождение строки А и столбец B) вызывает промах в кэше при доступе к элементу B. Это означает, что алгоритм требует Θ (п³) кеш промахивается в худшем случае. По состоянию на 2010 г., скорость памяти по сравнению со скоростью процессоров такова, что промахи в кэше, а не фактические вычисления, доминируют во времени работы для значительных матриц.^[3]

Оптимальный вариант итерационного алгоритма для А и B в строчном макете — это выложенный плиткой версия, где матрица неявно разделена на квадратные плитки размером √M от √M:^[3]^[4]

Вход: матрицы А и B
Позволять C быть новой матрицей подходящего размера
Выберите размер плитки Т = Θ (√M)
Для я от 1 до п в шагах от Т:
- Для J от 1 до п в шагах от Т:
  - Для K от 1 до м в шагах от Т:
    - Умножить А_{я:я+Т, K:K+Т} и B_{K:K+Т, J:J+Т} в C_{я:я+Т, J:J+Т}, это:
    - Для я от я к мин (я + Т, п):
      - Для j от J к мин (J + Т, п):
        Позволять сумма = 0
        Для k от K к мин (K + Т, м):
        Набор сумма ← сумма + А_ik × B_кДж
        Набор C_ij ← C_ij + сумма
Вернуть C

В идеализированной модели кэша этот алгоритм требует только Θ (п³/б √M) промахи в кэше; делитель б √M составляет несколько порядков на современных машинах, так что фактические вычисления преобладают над временем выполнения, а не промахами в кэше. ^[3]

Алгоритм разделяй и властвуй

Альтернативой итерационному алгоритму является разделяй и властвуй алгоритм для матричного умножения. Это зависит от блочное разбиение

C=(C11C12C21C22),А=(А11А12А21А22),B=(B11B12B21B22){displaystyle C = {egin {pmatrix} C_ {11} & C_ {12} C_ {21} & C_ {22} end {pmatrix}} ,, A = {egin {pmatrix} A_ {11} & A_ {12} A_ {21} & A_ {22} end {pmatrix}} ,, B = {egin {pmatrix} B_ {11} & B_ {12} B_ {21} & B_ {22} end {pmatrix}}},

который работает для всех квадратных матриц, размерность которых равна степени двойки, т.е. 2^п × 2^п для некоторых п. Матричный продукт теперь

(C11C12C21C22)=(А11А12А21А22)(B11B12B21B22)=(А11B11+А12B21А11B12+А12B22А21B11+А22B21А21B12+А22B22){displaystyle {egin {pmatrix} C_ {11} & C_ {12} C_ {21} & C_ {22} end {pmatrix}} = {egin {pmatrix} A_ {11} & A_ {12} A_ {21} & A_ {22} end {pmatrix}} {egin {pmatrix} B_ {11} & B_ {12} B_ {21} & B_ {22} end {pmatrix}} = {egin {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}}}

который состоит из восьми умножений пар подматриц с последующим этапом сложения. {2})},

учет восьми рекурсивных вызовов матриц размера п/2 и Θ (п²) для поэлементного суммирования четырех пар полученных матриц. Применение основная теорема для повторений «разделяй и властвуй» показывает эту рекурсию, чтобы получить решение Θ (п³), так же, как итерационный алгоритм.^[2]

Неквадратные матрицы

Вариант этого алгоритма, который работает для матриц произвольной формы и быстрее на практике^[3] разбивает матрицы на две вместо четырех подматриц следующим образом.^[5]Разделение матрицы теперь означает разделение ее на две части равного размера или как можно более близких к равным размерам в случае нечетных размеров.

Входы: матрицы А размера п × м, B размера м × п.
Базовый случай: если Максимум(п, м, п) ниже некоторого порога, используйте развернутый версия итерационного алгоритма.
Рекурсивные случаи:

Если Максимум(п, м, п) = п, Трещина А горизонтально:

C=(А1А2)B=(А1BА2B){displaystyle C = {egin {pmatrix} A_ {1} A_ {2} end {pmatrix}} {B} = {egin {pmatrix} A_ {1} B A_ {2} Bend {pmatrix}}}

Иначе, если Максимум(п, м, п) = п, Трещина B вертикально:

C=А(B1B2)=(АB1АB2){displaystyle C = A {egin {pmatrix} B_ {1} & B_ {2} end {pmatrix}} = {egin {pmatrix} AB_ {1} & AB_ {2} end {pmatrix}}}

В противном случае, Максимум(п, м, п) = м. Трещина А вертикально и B горизонтально:

C=(А1А2)(B1B2)=А1B1+А2B2{displaystyle C = {egin {pmatrix} A_ {1} & A_ {2} end {pmatrix}} {egin {pmatrix} B_ {1} B_ {2} end {pmatrix}} = A_ {1} B_ {1} + A_ {2} B_ {2}}

Поведение кеша

Частота промахов в кэше рекурсивного умножения матриц такая же, как у выложенный плиткой итерационная версия, но в отличие от этого алгоритма рекурсивный алгоритм не обращающий внимания на тайник:^[5] нет параметра настройки, необходимого для получения оптимальной производительности кеша, и он хорошо себя ведет в мультипрограммирование среда, в которой размеры кэша фактически динамичны из-за того, что другие процессы занимают кеш-пространство.^[3](Простой итерационный алгоритм также не учитывает кеш, но на практике он намного медленнее, если макет матрицы не адаптирован к алгоритму.)

Количество промахов в кэше, вызванных этим алгоритмом на машине с M строки идеального кэша, каждая размером б байтов, ограничено^[5]^:13

Θ(м+п+п+мп+пп+мпб+мппбM){displaystyle Theta left (m + n + p + {frac {mn + np + mp} {b}} + {frac {mnp} {b {sqrt {M}}}} ight)}

Субкубические алгоритмы

Низший ω такое, что умножение матриц, как известно, О(п^ω), построенный против времени. {2.807})}. Алгоритм Штрассена более сложен, и числовая стабильность уменьшается по сравнению с наивным алгоритмом,^[6]но это быстрее в случаях, когда п > 100 или так^[1] и присутствует в нескольких библиотеках, таких как BLAS.^[7] Это очень полезно для больших матриц в точных областях, таких как конечные поля, где числовая стабильность не является проблемой.

Электрический ток О(п^k) алгоритм с наименьшим известным показателем k является обобщением Алгоритм Копперсмита – Винограда имеющая асимптотическую сложность О(п^2.3728639), автор Франсуа Ле Галль.^[8] Алгоритм Ле Галля и алгоритм Копперсмита – Винограда, на котором он основан, похожи на алгоритм Штрассена: изобретен способ умножения двух k × k-матрицы с числом меньше k³ умножения, и этот метод применяется рекурсивно. Однако постоянный коэффициент, скрытый Обозначение Big O настолько велик, что эти алгоритмы имеют смысл только для матриц, которые слишком велики для обработки на современных компьютерах. ^[9]^[10]

Поскольку любой алгоритм умножения двух п × п-матрицы должны обрабатывать все 2п² элементов, существует асимптотическая нижняя оценка Ω (п²) операции. Раз доказал нижнюю оценку Ω (п² журнал(п)) для арифметических схем с ограниченными коэффициентами над действительными или комплексными числами.^[11]

Кон и другие. использовать такие методы, как алгоритмы Штрассена и Копперсмита – Винограда, в совершенно ином теоретико-групповой контекста, используя тройки подмножеств конечных групп, которые удовлетворяют свойству дизъюнктности, называемому свойство тройного продукта (TPP). Они показывают, что если семьи венки из Абелевы группы с симметричными группами реализуют семейства троек подмножеств с одновременной версией TPP, тогда существуют алгоритмы умножения матриц с существенно квадратичной сложностью.^[12]^[13] Большинство исследователей считают, что это действительно так. ^[10] Однако Алон, Шпилка и Крис Уманс недавно показали, что некоторые из этих гипотез, предполагающих быстрое матричное умножение, несовместимы с другой правдоподобной гипотезой, догадка подсолнечника.^[14]

Алгоритм Фрейвальдса это простой Алгоритм Монте-Карло что, учитывая матрицы А, B и C, проверяет в Θ (п²) время, если AB = C.

Параллельные и распределенные алгоритмы

Параллелизм с общей памятью

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

(А11B11+А12B21А11B12+А12B22А21B11+А22B21А21B12+А22B22){displaystyle {egin {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}}}

могут выполняться независимо друг от друга, как и четыре суммирования (хотя алгоритм должен «соединить» умножения перед выполнением суммирования). Используя полный параллелизм задачи, можно получить алгоритм, который может быть выражен в виде вилка – стиль соединения псевдокод:^[15]

Процедура умножить (C, А, B):

Базовый случай: если п = 1, набор c₁₁ ← а₁₁ × б₁₁ (или умножьте небольшую блочную матрицу).
В противном случае выделите место для новой матрицы Т формы п × п, тогда:
- Раздел А в А₁₁, А₁₂, А₂₁, А₂₂.
- Раздел B в B₁₁, B₁₂, B₂₁, B₂₂.
- Раздел C в C₁₁, C₁₂, C₂₁, C₂₂.
- Раздел Т в Т₁₁, Т₁₂, Т₂₁, Т₂₂.
- Параллельное исполнение:
  - Вилка умножить (C₁₁, А₁₁, B₁₁).
  - Вилка умножить (C₁₂, А₁₁, B₁₂).
  - Вилка умножить (C₂₁, А₂₁, B₁₁).
  - Вилка умножить (C₂₂, А₂₁, B₁₂).
  - Вилка умножить (Т₁₁, А₁₂, B₂₁).
  - Вилка умножить (Т₁₂, А₁₂, B₂₂).
  - Вилка умножить (Т₂₁, А₂₂, B₂₁).
  - Вилка умножить (Т₂₂, А₂₂, B₂₂).
- Присоединиться (дождитесь завершения параллельных вилок).
- Добавить(C, Т).
- Освободить Т.

Процедура Добавить(C, Т) добавляет Т в C, поэлементно:

Базовый случай: если п = 1, набор c₁₁ ← c₁₁ + т₁₁ (или сделайте короткий цикл, возможно, развернутый).
В противном случае:
- Раздел C в C₁₁, C₁₂, C₂₁, C₂₂.
- Раздел Т в Т₁₁, Т₁₂, Т₂₁, Т₂₂.
- В параллели:
  - Вилка Добавить(C₁₁, Т₁₁).
  - Вилка Добавить(C₁₂, Т₁₂).
  - Вилка Добавить(C₂₁, Т₂₁).
  - Вилка Добавить(C₂₂, Т₂₂).
- Присоединиться.

Вот, вилка — ключевое слово, которое сигнализирует, что вычисление может выполняться параллельно с остальной частью вызова функции, а присоединиться ожидает завершения всех ранее «разветвленных» вычислений. раздел достигает своей цели только манипуляциями с указателями.

Этот алгоритм имеет длина критического пути из Θ (журнал²п) шагов, то есть на идеальной машине с бесконечным количеством процессоров требуется столько времени; следовательно, он имеет максимально возможное ускорение из Θ (п³/журнал²п) на любом реальном компьютере. Алгоритм непрактичен из-за затрат на связь, связанных с перемещением данных во временную матрицу и из нее. Т, но в более практичном варианте достигается Θ (п²) ускорение без использования временной матрицы.^[15]

Блочное умножение матриц. В 2D-алгоритме каждый процессор отвечает за одну подматрицу C. В 3D-алгоритме каждая пара подматриц из А и B то, что умножается, назначается одному процессору.

Распределенные алгоритмы предотвращения общения

В современных архитектурах с иерархической памятью стоимость загрузки и хранения элементов матрицы ввода имеет тенденцию преобладать над затратами на арифметику. На одной машине это объем данных, передаваемых между ОЗУ и кешем, тогда как на многоузловой машине с распределенной памятью это объем, передаваемый между узлами; в любом случае это называется пропускная способность связи. Наивный алгоритм, использующий три вложенных цикла, использует Ω (п³) пропускная способность связи.

Алгоритм Кэннона, также известный как 2D алгоритм, это алгоритм предотвращения общения который разбивает каждую входную матрицу на блочную матрицу, элементы которой являются подматрицами размера √M/3 от √M/3, где M это размер быстрой памяти.^[16] Затем наивный алгоритм используется над блочными матрицами, вычисляя продукты подматриц полностью в быстрой памяти. Это снижает пропускную способность связи до О(п³/√M), что является асимптотически оптимальным (для алгоритмов, выполняющих Ω (п³) вычисление).^[17]^[18]

В распределенной среде с п процессоры расположены в √п от √п 2D-сетка, одна подматрица результата может быть назначена каждому процессору, и произведение может быть вычислено с каждым процессором, передающим О(п²/√п) слов, что является асимптотически оптимальным при условии, что каждый узел хранит минимум О(п²/п) элементы. ^[18] Это можно улучшить с помощью 3D алгоритм, который объединяет процессоры в трехмерную кубическую сетку, назначая каждое произведение двух входных подматриц одному процессору. Затем подматрицы результатов генерируются путем сокращения по каждой строке.^[19] Этот алгоритм передает О(п²/п^2/3) слов на процессор, что является асимптотически оптимальным.^[18] Однако это требует репликации каждого элемента входной матрицы. п^1/3 раз, и поэтому требуется фактор п^1/3 больше памяти, чем необходимо для хранения входных данных. Этот алгоритм можно комбинировать со Strassen для дальнейшего сокращения времени выполнения.^[19] Алгоритмы «2.5D» обеспечивают постоянный компромисс между использованием памяти и пропускной способностью связи.^[20] В современных распределенных вычислительных средах, таких как Уменьшение карты, разработаны специализированные алгоритмы умножения. ^[21]

Алгоритмы для сеток

Умножение матриц выполнено за 2n-1 шагов для двух матриц размера n × n на перекрестной сетке.

Существует множество алгоритмов умножения на сетки. Для умножения двух п×п на стандартной двумерной сетке с использованием 2D Алгоритм Кэннона, можно завершить умножение на 3п-2 шага, хотя при повторных вычислениях это число уменьшается вдвое.^[22] Стандартный массив неэффективен, потому что данные из двух матриц не поступают одновременно, и он должен быть дополнен нулями.

Результат будет еще быстрее на двухслойной сетке с перекрестной связью, где всего 2пТребуется -1 шаг.^[23] Производительность улучшается еще больше для повторных вычислений, что приводит к 100% эффективности.^[24] Сетчатый массив с перекрестными связями можно рассматривать как частный случай неплоской (то есть многослойной) структуры обработки.^[25]

Смотрите также

Вычислительная сложность математических операций
Алгоритм CYK, алгоритм §Valiant
Умножение матричной цепочки
Умножение разреженной матрицы на вектор

использованная литература

^ ^а ^б ^c ^d Скиена, Стивен (2008). Как, S (1988). «Многослойные массивные вычисления». Информационные науки. 45 (3): 347–365. CiteSeerX 10.1.1.90.4753. Дои:10.1016/0020-0255(88)90010-2.

дальнейшее чтение

Буттари, Альфредо; Лангу, Жюльен; Курзак, Якуб; Донгарра, Джек (2009). «Класс параллельных мозаичных алгоритмов линейной алгебры для многоядерных архитектур». Параллельные вычисления. 35: 38–53. arXiv:0709.1272. Дои:10.1016 / j.parco.2008.10.002. S2CID 955.
Гото, Казушигэ; ван де Гейн, Роберт А. (2008). «Анатомия высокопроизводительного матричного умножения». Транзакции ACM на математическом ПО. 34 (3): 1–25. CiteSeerX 10.1.1.140.3583. Дои:10.1145/1356052.1356053. S2CID 9359223.
Van Zee, Field G .; ван де Гейн, Роберт А. (2015). «BLIS: структура для быстрого создания экземпляра функциональности BLAS». Транзакции ACM на математическом ПО. 41 (3): 1–33. Дои:10.1145/2764454. S2CID 1242360.
Как оптимизировать GEMM

Алгоритм Штрассена

Алгоритм Штрассена.

Содержание

1	Вступление и постановка задачи	2
2	Алгоритм Штрассена	2
3	Умножение матриц по определению, встроенная функция
	matmul и оптимизация	4
4	Реализация алгоритма	8

1Вступление и постановка задачи

Одна из главных задач данной работы — реализовать алгоритм Штрассена. Это рекурсивный алгоритм для умножения матриц, который позволяет значительно ускорить процесс умножения. Надо заметить, что данный алгоритм эффективен только для матриц достаточно большого размера. Определить размер матриц, с которого имеет смысл начинать умножать по алгоритму Штрассена, можно экспериментальным путем. Кроме того, я уже упоминала, что алгоритм является рекурсивным. В нерекурсивной ветви алгоритма используется «простое»умножение матриц, т.е. умножение по определению. Сделать нерекурсивную ветвь наиболее эффективной и организовать выход из рекурсии в нужный момент — еще одна задача моей работы.

Итак, наметим цели работы:

Реализовать алгоритм Штрассена;

Экспериментально найти размер матрицы, для которого следует переходить на умножение по алгоритму;

Организовать выход из рекурсии в нужный момент;

Сделать не рекурсивную ветвь наиболее эффективной.

2Алгоритм Штрассена

Самая затратная по времени операция — умножение. Т.е. чтобы ускорить процесс, нужно сократить их число. Штрассен придумал алгоритм [1, 2, 3], который позволяет избежать одного умножения. Замечательно, что формулы не требуют коммутативности сложения, т.е. они приминимы и для матриц, что позволяет применять алгоритм рекурсивно.

Алгоритм работает с квадратными матрицами размера 2n. Заметим, что для умножения матриц такого размера требуется 23n умножений. Если же матрицы не квадратные, то можно расширить их до нужного размера, заполнив недостающие строки и столбцы нулями. Допустим, есть две матрицы A и B размера 2n, которые требуется перемножить: C = AB. Разобьем каждую матрицу на четыре подматрицы:

A2;1	A2;2	B2;1	B2;2	=	C2;1	C2;2
A1;1	A1;2	B1;1	B1;2		C1;1	C1;2

Тогда матрицы Ai;j и Bi;j имеют размер 2n 1. А элементы матрицы можно вычислить следующим образом:

C1;1 = A1;1B1;1 + A1;2B2;1

C1;2 = A1;1B1;2 + A1;2B2;2

C2;1 = A2;1B1;1 + A2;2B2;1

C2;2 = A2;1B1;2 + A2;2B2;2

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

Штрассен ввел семь новых элементов, которые вычисляются с помощью одного умножения, а матрица C получается из них сложением или вычитанием. Определим эти семь элементов:

M1 = (A1;1 + A2;2)(B1;1 + B2;2)

M2 = (A2;1 + A2;2)B1;1

M3 = A1;1(B1;2 B2;2)

M4 = A2;2(B2;1 B1;1)

M5 = (A1;1 + A1;2)B2;2

M6 = (A2;1 A1;1)(B1;1 + B1;2)

M7 = (A1;2 A2;2)(B2;1 + B2;2)

Элементы матрицы C вычисляются по формулам:

C1;1 = M1 + M4 M5 + M7

C1;2 = M3 + M5

C2;1 = M2 + M4

C2;2 = M1 M2 + M3 + M6

Теперь на каждом этапе рекурсии требуется выполнять только семь умножений вместо восьми. Предполагается, что итерационный процесс продолжается 2n раз, до тех пор, пока матрицы Ci;j не выродятся в числа. Заметим, что для маленьких матриц алгоритм теряет свою эффективность. Также этот алгоритм не имеет большой устойчивости [3], т. е. для вещественных матриц, в случае, если нужна очень высокая точность, не

стоит умножать по Штрассену. Этот существенный минус легко объяснить большим количесвом сложений и вычитаний. Но для целочисленных матриц, конечно, проблемы не возникает.

Рекурсивный алгоритм Штрассена умножает две матрицы за время (n2:81) [1], тогда как простое умножение справляется с этой задачей за (n3). Для матриц большого размера(начиная примерно с 2000) эта разница существенна.

3Умножение матриц по определению, встроенная функция matmul и оптимизация

В этом разделе пойдет речь об умножении матриц без помощи алгоритма Штрассена. Прежде чем приступать непосредственно к алгоритму, нужно выяснить, нельзя ли ускорить процесс умножения с помощью других возможностей фортрана. Также в фортране есть специальная встроенная функция matmul, которая возвращает результат умножения двух матриц. Каким образом эта функция производит умножение, неизвестно. Попробуем написать код для умножения матриц по опредлению. Пусть A и B — умножаемые матрицы, а C — их произведение.

do i = 1, N

do j = 1, N C (i , j) = 0

do k = 1, N

C(i , j) = C(i , j) + A(i , k )* B(k , j)

end do

Назовем такой вариант def mul. Посмотрим, как быстро умножает матрицы matmul и как медленно происходит умножение по определению. В первой строке таблицы указан размер умножаемых матриц, во второй — время, за которое умножает matmul, а втретьей — время умножения по определению. Время указано в секундах, т.о. матрицы размера 2048 умножаются в первом случае — за 13 секунд, а во втором — за 15 минут!

size	128	256	512	1024	2048
matmul	4.00000019E-03	2.00009998E-02	0.21601403	1.7041068	13.596855
def mul	2. 80019976E-02	0.28801799	6.7804241	98.874176	900.42432

Откуда возникает такая разница? В фортране обращение к элементам массива происходит по столбцам [4]. Т.е., к примеру, чтобы обратиться к первому элементу пятого столбца, компилятор сначала просчитает все впереди стоящие элементы по столбцам. Значит, имеет смысл попробовать переписать умножение так, чтобы индекс столбца менялся быстрее,

чем индекс строки. Для этого надо поменять внешний и внутренний циклы местами, а средний оставить на месте:

do k = 1, N

do j = 1, N

do i = 1, N

C(i , j) = C(i , j) + A(i , k )* B(k , j)

end do

Теперь к предыдущей таблице добавим результат получившейся программы, назовем ее def mul tr.

size	128	256	512	1024	2048
matul	4. 00000019E-03	2.00009998E-02	0.21601403	1.7041068	13.596855
def mul tr	4.00019996E-02	0.26401600	2.1721358	17.385088	138.56866
def mul	2.80019976E-02	0.28801799	6.7804241	98.874176	900.42432

Нетрудно догадаться, что транпонировав одну из матриц, можно перемножать только столбцы, что еще сократит время умножения. Для транспонирования можно воспользоваться встроенной функцией transpose. Далее можно избавившись от одного цикла. Это позволяет сделать возможным описание сечений (в данном случае — столбцов) и использование встроенной функции сложения массивов sum. Тогда умножение двух матриц приобретёт вид:

A = transpose (A) do i = 1, N

do j = 1, N

C(j , i) = sum (A(: , j )* B(: , i ))

end do

Посмотрим, на что способен такой вариант(его название — two loops mul):

size	128	256	512	1024	2048
matmul	4. 00000019E-03	2.00009998E-02	0.21601403	1.7041068	13.596855
two loops mul	8.00000131E-03	5.60030043E-02	0.49202991	3.6642284	29.393829
def mul tr	4.00019996E-02	0.26401600	2.1721358	17.385088	138.56866
def mul	2.80019976E-02	0.28801799	6.7804241	98.874176	900.42432

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

Чтобы усовершенствовать наш метод, нужно обратить внимание на компилятор. Что вообще такое компилятор? Это программа, которая переводит текст с языка программирования высокого уровня(в данном случае на Фортране) на язык машинного уровня. Компилятор позволяет осуществить связь между тем, что пишет программист, и процессором компьютера. Я пользуюсь компилятором GF ortran из GNU Compiler Collection, который позволяет оптимизировать компиляцию[5]. Существует несколько уровней оптимизации: O1; O2; O3; Ofast. Каждый имеет свои особенности, которые я не буду подробно описывать. Для умножения матриц и алгоритма Штрассена лучше всего подходит уровень Ofast, т.к. он включает в себя наибольшее количество опций и является самым быстрым. Также можно использовать опцииmarch = native и mtune = native [6], которые генерируют код для указанной архитектуры процессора. На место native компилятор сам ставит значение, соответствующее данному процессору.

До этого момента я использовала стандартную компиляцию:

gfortran -o test_mul test_mul . f90

А теперь при компиляции добавлю вышеописанные опции. Тогда компиляция будет выглядеть так:

gfortran -o test_mul — Ofast — mtune = native — march = native test_mul . f90

Посмотрим, как изменятся результаты:

size	128	256	512	1024	2048
two loops mul	0. 0000000	8.00000131E-03	9.20059681E-02	0.69604301	5.7083435
def mul tr	0.0000000	1.20010078E-02	0.16401005	1.2800789	9.8126221
matmul	4.00000019E-03	2.40010004E-02	0.21201299	1.6761041	13.372837
def mul	1.20009994E-02	8.80059898E-02	1.9721229	22.437403	220.77379

Matmul занял на этот раз третье место, и время его работы почти не изменилось. Это объясняется тем, что эта встроенная функция уже скомпилирована, а вот написанные мной программы — нет, и на них опции компиляции влияют очень сильно. Теперь первое место в таблице занимает умножение с двумя циклами, второе место досталось умножению по столбцам, а умножение по определению занимает почётное последнее место. Последние результаты показывают, что преуспеть можно и без алгоритма Штрассена, а лишь используя знания о том, как устроен язык программирования. Опции компилятора Ofast; mtune = native; march = native сыграли решающую роль в оптимизации, что говорит о том, насколько важно собирать программу под конкретный процессор. В дальнейшем я всегда буду использовать эти три опции. Самое быстрое умножение с двумя циклами лучше всего подойдет для нерекурсивной ветви алгоритма Штрассена, так что часть задачи уже решена.

4Реализация алгоритма

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

все остальные матрицы нужно пополнить до нужного размера нулевыми строками и столбцами. Как это будет реализовано непосредственно на языке программирования, не имеет значения. После совершения этой операции мы можем иметь дело с матрицами любого размера. Написать рекурсивную подпрограмму, выполняющую алгоритм Штрассена, не очень сложно: для этого надо просто воспроизвести все то, что описано в соответствующем разделе, на языке программирования. Пусть есть две матрицы A и B, заполненные случайными числами, которые требуется перемножить, и пусть эти матрицы уже имеет размер 2n, а матрица C — результат умножения. В ходе алгоритма потребуется семь новых элементов M1 M7, которые будут вычисляться рекурсивно. Пока оставим вопрос о том, начиная с какого номера лучше выходить из рекурсии, будем умножать по алгоритму, пока матрица не станет размера 8 8.

recursive function strassen (A , B ,			N)	result (C)
real C(N , N), AT (N ,		N)
integer , intent ( in ) :: N
integer P , i , j ,	k
real , intent ( in )	:: A(N , N), B(N ,		N)

real , dimension (: ,:) , allocatable :: A11 , A12					,	A21 ,	A22
real ,	dimension (: ,:) ,	allocatable ::	B11 ,	B12	,	B21 ,	B22
real ,	dimension (: ,:) ,	allocatable ::	M1 ,	M2 ,	M3 , M4 ,		M5	, M6 , M7

if (N == 8) then
AT = transpose (A)
do i = 1, N
do j	= 1, N
	C(j , i) = sum ( AT (: , j )* B(: , i ))
end	do
end do
return
end if

P=N /2

allocate ( A11 (P ,P), A12 (P ,P), A21 (P ,P), A22 (P ,P )) allocate ( B11 (P ,P), B12 (P ,P), B21 (P ,P), B22 (P ,P ))

A11	=	A	(1:P , 1: P)
A12	=	A	(1:P , P +1: N)
A21	=	A	(P +1:N , 1: P)
A22	=	A	(P +1:N , P +1: N)
B11	=	B	(1:P , 1: P)
B12	=	B	(1:P , P +1: N)
B21	=	B	(P +1:N , 1: P)
B22	=	B	(P +1:N , P +1: N)

allocate ( M1 (P ,P), M2 (P ,P), M3 (P ,P), M4 (P ,P )) allocate ( M5 (P ,P), M6 (P ,P), M7 (P ,P ))


M1	=	strassen ( A11 + A22 , B11 + B22 , P)
M2	=	strassen ( A21 + A22 , B11 , P)
M3	=	strassen (A11 , B12			— B22 , P)
M4	=	strassen (A22 , B21			— B11 , P)
M5	=	strassen ( A11	+	A12 ,		B22 ,	P)
M6	=	strassen ( A21	—	A11 ,		B11	+ B12 , P)
M7	=	strassen ( A12	—	A22 ,		B21	+ B22 , P)

deallocate (A11 , A12 , A21 , A22 , B11 , B12 , B21 , B22 )

C	(1:P ,	1: P) = M1		+	M4	— M5 + M7
C	(1:P ,	P +1: N)	=	M3	+	M5
C	(P +1:N , 1: P)		=	M2	+	M4

C (P +1:N , P +1: N) = M1 — M2	+ M3 + M6
deallocate (M1 , M2 , M3 , M4 ,	M5 , M6 , M7 )
end function

Разберем тут же еще один вариант. Этот вариант не предполагает объявления семи матриц M1 M7, и тогда четыре блока матрицы C можно найти следующим образом:

(2.81)).

Вскоре после объявления результата Штрассена С. Виноград нашел способ умножить матрицы 2 X 2, используя 7 умножений и 15 добавления/удаления. Это наблюдение улучшает константу, но делает не меняет асимптотического времени вычислений. Впоследствии Виноград показал, что 7 умножений и 15 сложений/вычитаний оптимальны для матрицы 2 X 2 умножение. Его подход основан на следующих уравнениях.

[ А11 А12 ] [ В11 В12 ] [ С11 С12 ]
Пусть А = [ ] и B = [ ] и С = [ ] = А*В.
[ А21 А22] [ В21 В22 ] [ С21 С22 ]

S1 = A21 + A22, T1 = B12 – B11
S2 = S1 – A11, T2 = B22 – T1
S3 = A11 – A21, T3 = B22 – B12
S4 = A12 – S2, T4 = B21 – T2

Р1 = А11*В11
Р2 = А12*В21
Р3 = С1*Т1
Р4 = С2*Т2
Р5 = С3*Т3
Р6= С4*В22
Р7 = А22*Т4

9к. Запишите и решите рекуррентное соотношение для числа операций, используемых алгоритмом Штрассена и вариантом Винограда.

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

Вычислить отношение количества операций, используемых вариантом Винограда алгоритма Штрассена к стандартному алгоритму. Постройте отношение как k становится достаточно большим, так что асимптотическое поведение становится очевидным. Постройте соотношение. Когда алгоритм Штрассена становится быстрее?

Каково предельное соотношение числа операций, используемых Штрассеном? оригинальный алгоритм по сравнению с вариантом Винограда?

Точка пересечения, отмеченная в предыдущем вопросе, убедила многих люди, что алгоритм Штрассена (вариант Винограда) непрактичен; однако, если используется гибридный алгоритм, где обычный алгоритм используется, когда это быстрее, чем применение алгоритма Штрассена, тогда крест над точкой намного меньше. Найдите оптимальный гибридный алгоритм и определите точка, в которой он становится быстрее обычного алгоритма.

k — 1,9к.

На каждом рекурсивном шаге алгоритма Штрассена, если n нечетно, добавьте дополнительный строка и столбец нулей, так что A и B вложены в (n+1) X (n+1) матрицы А’ и Б’. Если n четное, встраивание не требуется.

На каждом рекурсивном шаге алгоритма Штрассена, если n нечетное представление A и B как матрицы 2 X 2, у которых верхний левый блок A11 и B11 равен (n-1) X (n-1) матрицы. Выполните умножение A11*B11 с помощью рекурсивный вызов алгоритма Штрассена. Остальные умножения равны выполненный используя обычное умножение матриц.

Алгоритм Штрассена для умножения матриц: моделирование, анализ и реализация Стивен

Идентификатор корпуса: 8734043

  title={Алгоритм Штрассена для умножения матриц: моделирование, анализ и реализация, Стивен},
  автор = {- Гус и Ледерман, Элейн и М. и Джейкобсон Дж},
  год = {1996}
}

-. Huss, LedermanElaine, JacobsonJ
Опубликовано в 1996 г.
Информатика

В этой статье мы сообщаем о разработке эффективной и портативной реализации алгоритма умножения матриц Штрассена для матриц произвольного размера. Наша реализация предназначена для использования вместо DGEMM, процедуры умножения матриц BLAS уровня 3. Наш код разработан таким образом, чтобы обеспечить эффективную производительность для всех размеров и форм матриц, а дополнительная память, необходимая для временных переменных, была сведена к минимуму. Замена DGEMM нашей процедурой должна обеспечить значимость…

Реализация алгоритма Штрассена для умножения матриц

Реализация предназначена для использования вместо DGEMM, процедуры умножения матриц BLAS уровня 3, и подтверждает, что алгоритм Штрассена применим для матриц реалистичного размера.

2356-5608 Эффективные рекурсивные реализации операций линейной алгебры

Р. Махфуди, Сэми Ачур, З. Махджуб
Информатика
2014

0082 В этой статье разработаны рекурсивные алгоритмы D&C для решения четырех основных задач линейной алгебры, а именно: умножение матриц (MM), решение системы треугольных матриц, LU-факторизация, решение системы плотных матриц и алгоритм Штрассена.
Adaptive Winograd’s matrix multiplications
- P. D’Alberto, A. Nicolau
- Computer Science
  TOMS
- 2009
This work presents a novel, hybrid, and adaptive recursive Strassen-Winograd’s matrix multiplication ( MM), который использует автоматически настраиваемое программное обеспечение линейной алгебры (ATLAS) или GotoBLAS и представляет экспериментальные доказательства, основанные на установленных методологиях, найденных в литературе, что алгоритм для семейства матриц так же точен, как и классические алгоритмы.
Основа для практического параллельного быстрого умножения матриц
Показано, что новые алгоритмы быстрого умножения матриц могут значительно превзойти реализации поставщиками классического алгоритма и быстрого алгоритма Штрассена для задач небольшого размера и формы, и что лучший выбор быстрого алгоритма зависит от только от размера матриц, но и от формы.
Матричное умножение по методу Штрассена, подходящее для возведения в квадрат и более мощных вычислений
- Марко Бодрато
- Информатика, математика
  ISSAC
- 2010
Новый многовариантный вариант метода Штрассена для матричного умножения, который может сохранить, по крайней мере, не хуже, чем уже известные простые операции с матрицами. либо для цепных изделий, либо для возведения в квадрат, и может оказаться оптимальным для этих задач.
Фреймворк для высокопроизводительного умножения матриц на основе иерархических абстракций, алгоритмов и оптимизированных низкоуровневых ядер
В этой статье сделано несколько усовершенствований алгоритмов, включая осциллирующий итерационный алгоритм для умножения матриц и переменный критерий отсечения рекурсии для алгоритма Штрассена, а также необходимость стандартизации интерфейсов ядра линейной алгебры, отличных от BLAS, для написания портативных высокопроизводительных код выставлен.
Планирование алгоритма умножения матриц Штрассена-Винограда с эффективным использованием памяти
- Брис Бойер, Ж. Дюма, Клеман Перне, В. Чжоу
- Информатика
  ISSAC ’09
- 2009
Мы предлагаем несколько новых графиков для алгоритма умножения матриц Штрассена-Винограда. Умножение матриц на графических процессорах
В этой работе представлены эффективные реализации алгоритма Штрассена, а также варианта Винограда относительно кода sgemm в CUBLAS 3.0 с одинарной точностью и целочисленным графическим процессором при умножении матриц 16384×16384.
Смешанные параллельные реализации шага верхнего уровня алгоритмов Штрассена и Винограда для умножения матриц
В этой статье представлены параллельные реализации верхнего уровня алгоритмов Штрассена и Винограда для умножения матриц, которые используют смешанный параллелизм, т. е. одновременное использование данных и …
Плотная линейная алгебра над простыми полями размером в слово: пакеты FFLAS и FFPACK
- Дж. Дюма, Паскаль Джорджи, Клеман Перне
- Информатика, математика
  TOMS
- 2008
В этой статье рассматриваются высокопроизводительные реализации основных подпрограмм линейной алгебры над простыми полями размером слова: особенно умножение матриц, с целью предоставить точную альтернативу числовой библиотеке BLAS. .
ПОКАЗАНЫ 1-10 ИЗ 27 ССЫЛОК
СОРТИРОВАТЬ ПОРелевантности Наиболее влиятельные документыНедавность
Использование быстрого матричного умножения на уровне 3 BLAS
- Н. Хайэм
- Информатика
  TOMS
- 1990
Алгоритм для операций BLAS3, который асимптотически быстрее, чем обычный метод умножения для матриц, теперь распознается методом Фастсена, основанным на Страссене быть практически полезным методом, когда размеры матрицы превышают примерно 100.
Сверхскоростное умножение матриц на Cray-2
- D. Bailey
- Business
- 1988
В этой статье описывается, как умножение матриц может быть выполнено еще быстрее, в два раза выше, чем указано выше, за счет использования алгоритма умножения матриц Штрассена для уменьшения количества выполняемых операций с плавающей запятой и использования локальная память на Cray-2, чтобы избежать потери производительности из-за конфликта банков памяти.
Быстрое умножение прямоугольных матриц и QR-разложение
- П. Найт
- Математика
- 1995
Гауссовая ликвидация не является оптимальной
- V. Strassen
- Математика
- 1969