Есть ли алгоритм умножения квадратных матриц на месте?
Наивный алгоритм умножения матриц 4x4 выглядит следующим образом:
void matrix_mul(double out[4][4], double lhs[4][4], double rhs[4][4]) {
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
out[i][j] = 0.0;
for (int k = 0; k < 4; ++k) {
out[i][j] += lhs[i][k] * rhs[k][j];
}
}
}
}
Очевидно, что этот алгоритм дает фиктивные результаты, если out == lhs или out == rhs (здесь == означает ссылочное равенство). Есть ли версия, которая позволяет одному или двум из этих случаев не просто копировать матрицу? Я рад иметь различные функции для каждого случая, если это необходимо.
Я нашел этот документ, но в нем обсуждается алгоритм Штрассен-Винограда, который является излишним для моих маленьких матриц. Ответы на этот вопрос, похоже, указывают на то, что если out == lhs && out == rhs (т.е. Мы пытаемся квадрат матрицы), то это невозможно сделать на месте, но даже там нет убедительных доказательств или доказательств.
Ответы
Ответ 1
Я не в восторге от этого ответа (я публикую его в основном, чтобы заставить замолчать "это, очевидно, не может быть сделано" ), но я скептически отношусь к тому, что можно сделать гораздо лучше, (O (1) дополнительные слова памяти для умножения двух nxn-матриц). Пусть назовем две матрицы мультипликаторами A и B. Предположим, что A и B не сглажены.
Если A было верхнетреугольным, тогда проблема умножения будет выглядеть так.
[a11 a12 a13 a14] [b11 b12 b13 b14]
[ 0 a22 a23 a24] [b21 b22 b23 b24]
[ 0 0 a33 a34] [b31 b32 b33 b34]
[ 0 0 0 a44] [b41 b42 b43 b44]
Мы можем вычислить произведение в B следующим образом. Умножьте первую строку B на a11. Добавьте a12 раз второй ряд B к первому. Добавьте a13 раз третью строку B к первой. Добавьте a14 раз четвертую строку B к первой.
Теперь мы перезаписали первую строку B с правильным продуктом. К счастью, нам это больше не нужно. Умножьте вторую строку B на a22. Добавьте a23 к третьей строке B ко второй. (Вы получаете идею.)
Аналогично, если A было единицей ниже треугольника, то проблема умножения будет выглядеть следующим образом.
[ 1 0 0 0 ] [b11 b12 b13 b14]
[a21 1 0 0 ] [b21 b22 b23 b24]
[a31 a32 1 0 ] [b31 b32 b33 b34]
[a41 a42 a43 1 ] [b41 b42 b43 b44]
Добавьте a43 раз к третьей строке B к четвертой. Добавьте a42 раз, когда вторая строка B - четвертая. Добавьте a41 раз, когда первая строка B будет четвертой. Добавьте a32 раз вторую строку B к третьей. (Вы получаете идею.)
Полный алгоритм состоит в том, чтобы LU-разложить A на месте, умножить UB на B, умножить LB на B и затем LU-undecompose A на месте (я не уверен, что кто-нибудь когда-либо делает это, но это кажется достаточно простым для отмены шагов). Существует около миллиона причин не реализовывать это на практике, два из них состоят в том, что A не может быть LU-разложимым и что A не будет реконструироваться точно в общем случае с арифметикой с плавающей запятой.
Ответ 2
Этот ответ более разумен, чем мой другой, хотя он использует один целый столбец дополнительного хранилища и имеет такое же количество движения данных, как и алгоритм наивного копирования. Чтобы умножить A на B, сохранить продукт в B (опять же, считая, что A и B хранятся отдельно):
For each column of B,
Copy it into the auxiliary storage column
Compute the product of A and the auxiliary storage column into that column of B
Я переключил псевдокод, чтобы сделать копию первой, потому что для больших матриц эффекты кэширования могут привести к более эффективному умножению A на смежный вспомогательный столбец в отличие от несмежных записей в B.
Ответ 3
Этот ответ о матрицах 4x4. Предполагая, что, как вы полагаете, out может ссылаться либо на lhs, либо на rhs, и что A и B имеют ячейки равномерной длины бит, чтобы технически быть в состоянии выполнить умножение на месте, элементы A и B, как знаковые целые числа, обычно не могут быть больше или меньше, чем ± floor (sqrt (2 ^ (cellbitlength - 1) / 4)).
В этом случае мы можем взломать элементы A в B (или наоборот) в виде битового сдвига или комбинации битовых флагов и модульной арифметики и вычислить продукт в прежнюю матрицу. Если A и B были плотно упакованы, за исключением особых случаев или лимитов, мы не могли признать out ссылкой на lhs или rhs.
Использование наивного метода теперь не будет отличаться от описания алгоритма Дэвида, как раз с дополнительным столбцом, хранящимся в или B. В качестве альтернативы мы могли бы реализовать алгоритм Штрассен-Виноград согласно приведенному ниже графику, опять же без хранения вне lhs и rhs. (Формулировка p0,...,p6 и C взята со страницы 166 Джонатана Голана. Линейная алгебра а начинающий выпускник должен знать.)
p0 = (a11 + a12)(b11 + b12), p1 = (a11 + a22)b11, p2 = a11(b12 - b22),
p3 = (a21 - a11)(b11 + b12), p4 = (a11 + a12)b22, p5 = a22(b21 - b11),
p6 = (a12 - a22)(b21 + b22)
┌ ┐
c = │ p0 + p5 - p4 + p6, p2 + p4 │
│ p1 + p5 , p0 - p1 + p2 + p3 │
└ ┘
Расписание:
Каждый p ниже представляет собой квадрант 2x2; "x" означает неназначение; "nc", никаких изменений. Чтобы вычислить каждый p, мы используем неназначенный квадрант 2x2 для наложения (одного или двух) результатов сложения или вычитания блочной матрицы 2x2 с использованием того же метода сдвига бит или модуляции выше; мы затем добавляем их произведение (семь умножений, приводящих к одиночным элементам) непосредственно в целевой блок в любом порядке (обратите внимание, что для 2x2-размера p2 и p4 мы используем юго-западный квадрант rhs, который в этой точке больше не нужен). Например, чтобы записать первый 2x2-размер p6, мы накладываем вычитание блочной матрицы, rhs(a12) - rhs(a22) и добавление блочной матрицы rhs(b21) + rhs(b22) на подматрицу lhs21; затем добавьте каждый из семи одиночных элементов p для этого блочного умножения, (a12 - a22) X (b21 + b22), непосредственно в подматрицу lhs11.
LHS RHS (contains A and B)
(1)
p6 x
x p3
(2)
+p0 x
p0 +p0
(3)
+p5 x
p5 nc
(4)
nc p1
+p1 -p1
(5)
-p4 p4 p4 (B21)
nc nc
(6)
nc +p2 p2 (B21)
nc +p2
