Прежде чем перейти к статье, хочу вам представить, экономическую онлайн игру Brave Knights, в которой вы можете играть и зарабатывать. Регистируйтесь, играйте и зарабатывайте!
Нас учат, что чтение данных из оперативной памяти — ужасно долгая операция. Приводят аналогии с офисом и удалённым складом, заставляют писать cache-friendly код и внушают смертельный страх перед промахами кэша. Ещё нас учат, что процессоры отлично умеют считать числа, и часто быстрее вычислить результат дважды, чем сохранять его в памяти. Оказывается, это не всегда так.
Эта статья основана на реальном проекте и реальном коде, который был ускорен с помощью кэша почти в полтора раза. Весь код написан на JavaScript.
Задача
Допустим, у нас есть матрица A порядка 2000x2000. Нужно посчитать обратную ей матрицу по простому модулю N. Другими словами, надо найти такую матрицу A-1, что AA-1 mod N = E.
Поскольку вычисления у нас происходят в поле по модулю, итерационные методы обращения нам не подойдут. Будем использовать старый добрый метод Гаусса.
Этот пост посвящён оптимизации метода Гаусса под данный конкретный случай. В реальном проекте вычисление обратной матрицы происходит в отдельном WebWorker, в данном примере обойдёмся главным потоком.
Вспомогательные функции
Для работы программы нам потребуется четыре вспомогательные функции. Первая — вычисление (1 / x) mod N по расширенному алгоритму Евклида:
function invModGcdEx(x, domain)
{
if(x === 1)
{
return 1;
}
else
{
//В случае 0 или делителя нуля возвращается 0, означающий некий "некорректный результат"
if(x === 0 || domain % x === 0)
{
return 0;
}
else
{
//Расширенный алгоритм Евклида, вычисляющий такое число tCurr, что tCurr * x + rCurr * N = 1
//Другими словами, существует такое число rCurr, при котором tCurr * x mod N = 1
let tCurr = 0;
let rCurr = domain;
let tNext = 1;
let rNext = x;
while(rNext !== 0)
{
let quotR = Math.floor(rCurr / rNext);
let tPrev = tCurr;
let rPrev = rCurr;
tCurr = tNext;
rCurr = rNext;
tNext = Math.floor(tPrev - quotR * tCurr);
rNext = Math.floor(rPrev - quotR * rCurr);
}
tCurr = (tCurr + domain) % domain;
return tCurr;
}
}
}Вторая — корректное целочисленное деление по модулю. Наивное вычисление c = a % b во всех языках программирования не будет давать математически верный результат, если a — отрицательное число. Поэтому заведём функцию, которая будет делить правильно:
function wholeMod(x, domain)
{
return ((x % domain) + domain) % domain;
}Последние две функции относятся к операциям над строками матрицы. Первая — вычитание из строки матрицы домноженной на число другой:
function mulSubRow(rowLeft, rowRight, mulValue, domain)
{
for(let i = 0; i < rowLeft.length; i++)
{
rowLeft[i] = wholeMod(rowLeft[i] - mulValue * rowRight[i], domain);
}
}Последняя нужная нам функция — умножение строки матрицы на число:
function mulRow(row, mulValue, domain)
{
for(let i = 0; i < row.length; i++)
{
row[i] = (row[i] * mulValue) % domain;
}
}Обращение матрицы
Начнём с обычной наивной реализации. Создаём единичную матрицу, проходим прямым ходом, потом проходим обратным. На каждом шаге производим одинаковые операции над исходной матрицей и над только что созданной единичной.
function invertMatrix(matrix, domain)
{
let matrixSize = matrix.length;
//Инициализируем обратную матрицу единичной
let invMatrix = [];
for(let i = 0; i < matrixSize; i++)
{
let matrixRow = new Uint8Array(matrixSize);
matrixRow.fill(0);
matrixRow[i] = 1;
invMatrix.push(matrixRow);
}
//Прямой ход: приведение матрицы к ступенчатому виду
for(let i = 0; i < matrixSize; i++)
{
let thisRowFirst = matrix[i][i];
if(thisRowFirst === 0 || (thisRowFirst !== 1 && domain % thisRowFirst === 0)) //Первый элемент строки 0 или делитель нуля, меняем строку местами со следующей строкой, у которой первый элемент не 0
{
for(let j = i + 1; j < matrixSize; j++)
{
let otherRowFirst = matrix[j][i];
if(otherRowFirst !== 0 && (otherRowFirst === 1 || domain % otherRowFirst !== 0)) //Нашли строку с ненулевым первым элементом
{
thisRowFirst = otherRowFirst;
let tmpMatrixRow = matrix[i];
matrix[i] = matrix[j];
matrix[j] = tmpMatrixRow;
let tmpInvMatrixRow = invMatrix[i];
invMatrix[i] = invMatrix[j];
invMatrix[j] = tmpInvMatrixRow;
break;
}
}
}
//Обнуляем первые элементы всех строк после первой, отнимая от них (otherRowFirst / thisRowFirst) * x mod N
let invThisRowFirst = invModGcdEx(thisRowFirst, domain);
for(let j = i + 1; j < matrixSize; j++)
{
let otherRowFirst = matrix[j][i];
let mulValue = invThisRowFirst * otherRowFirst;
if(otherRowFirst !== 0 && (otherRowFirst === 1 || domain % otherRowFirst !== 0))
{
mulSubRow(matrix[j], matrix[i], mulValue, domain);
mulSubRow(invMatrix[j], invMatrix[i], mulValue, domain);
}
}
}
//Обратный ход - обнуление всех элементов выше главной диагонали
let matrixRank = matrixSize;
for(let i = matrixSize - 1; i >= 0; i--)
{
let thisRowLast = matrix[i][i];
let invThisRowLast = invModGcdEx(thisRowLast, domain);
for(let j = i - 1; j >= 0; j--)
{
let otherRowLast = matrix[j][i];
let mulValue = invThisRowLast * otherRowLast;
if(otherRowLast !== 0 && (otherRowLast === 1 || domain % otherRowLast !== 0))
{
mulSubRow(matrix[j], matrix[i], mulValue, domain);
mulSubRow(invMatrix[j], invMatrix[i], mulValue, domain);
}
}
if(thisRowLast !== 0 && domain % thisRowLast !== 0)
{
mulRow(matrix[i], invThisRowLast, domain);
mulRow(invMatrix[i], invThisRowLast, domain);
}
if(matrix[i].every(val => val === 0))
{
matrixRank -= 1;
}
}
return {inverse: invMatrix, rank: matrixRank};
}Проверим скорость на матрице 500 x 500, заполненной случайными значениями из поля Z / 29. После 5 испытаний получаем среднее время выполнения в ~9.4с. Можем ли мы сделать лучше?
Первое, что мы можем заметить — в поле Z / N не больше N обратных элементов. Чтобы избежать многократного вызова алгоритма Евклида, мы можем вычислить все обратные значения один раз и при надобности брать уже готовые. Изменим функцию соответствующим образом:
function invertMatrixCachedInverses(matrix, domain)
{
let matrixSize = matrix.length;
//Инициализируем обратную матрицу единичной
let invMatrix = [];
for(let i = 0; i < matrixSize; i++)
{
let matrixRow = new Uint8Array(matrixSize);
matrixRow.fill(0);
matrixRow[i] = 1;
invMatrix.push(matrixRow);
}
//Вычисляем все обратные элементы заранее
let domainInvs = [];
for(let d = 0; d < domain; d++)
{
domainInvs.push(invModGcdEx(d, domain));
}
//Прямой ход: приведение матрицы к ступенчатому виду
for(let i = 0; i < matrixSize; i++)
{
let thisRowFirst = matrix[i][i];
if(domainInvs[thisRowFirst] === 0) // <--- Первый элемент строки 0 или делитель нуля, меняем строку местами со следующей строкой, у которой первый элемент не 0
{
for(let j = i + 1; j < matrixSize; j++)
{
let otherRowFirst = matrix[j][i];
if(domainInvs[otherRowFirst] !== 0) // <--- Нашли строку с ненулевым первым элементом
{
thisRowFirst = otherRowFirst;
let tmpMatrixRow = matrix[i];
matrix[i] = matrix[j];
matrix[j] = tmpMatrixRow;
let tmpInvMatrixRow = invMatrix[i];
invMatrix[i] = invMatrix[j];
invMatrix[j] = tmpInvMatrixRow;
break;
}
}
}
//Обнуляем первые элементы всех строк после первой, отнимая от них (otherRowFirst / thisRowFirst) * x mod N
let invThisRowFirst = domainInvs[thisRowFirst]; // <---
for(let j = i + 1; j < matrixSize; j++)
{
let otherRowFirst = matrix[j][i];
let mulValue = invThisRowFirst * otherRowFirst;
if(domainInvs[otherRowFirst] !== 0) // <---
{
mulSubRow(matrix[j], matrix[i], mulValue, domain);
mulSubRow(invMatrix[j], invMatrix[i], mulValue, domain);
}
}
}
//Обратный ход - обнуление всех элементов выше главной диагонали
let matrixRank = matrixSize;
for(let i = matrixSize - 1; i >= 0; i--)
{
let thisRowLast = matrix[i][i];
let invThisRowLast = domainInvs[thisRowLast]; // <---
for(let j = i - 1; j >= 0; j--)
{
let otherRowLast = matrix[j][i];
let mulValue = invThisRowLast * otherRowLast;
if(domainInvs[otherRowLast] !== 0) // <---
{
mulSubRow(matrix[j], matrix[i], mulValue, domain);
mulSubRow(invMatrix[j], invMatrix[i], mulValue, domain);
}
}
if(domainInvs[thisRowLast] !== 0) // <---
{
mulRow(matrix[i], invThisRowLast, domain);
mulRow(invMatrix[i], invThisRowLast, domain);
}
if(matrix[i].every(val => val === 0))
{
matrixRank -= 1;
}
}
return {inverse: invMatrix, rank: matrixRank};
}Замерим на тех же условиях и получаем результат в те же ~9.4с. Прироста нет, потому что даже при относительно долгом вычислении алгоритма Евклида он вычисляется всего один раз для каждой строки матрицы и особого вклада во время не приносит. Замерим производительность и посмотрим, что ещё можно улучшить.
72% времени занимает деление по модулю при сложении строк матрицы! Ну что тут сказать, деление по модулю, пусть и немного модифицированное для отрицательных чисел — это элементарная операция и ускорять её некуда. Алгоритм поменять тоже не получится, из чего мы делаем вывод, что дальнейшее улучшение невозможно и статью можно закрывать.
...Или всё же возможно?
Если деление по модулю занимает столько времени, может, все возможные результаты тоже стоит сохранить в кэш? Даже если это не поможет, попытаться всё равно стоит — при текущем времени выполнения функция неюзабельна.
Итак, используется wholeMod()только в функции mulSubRow():
rowLeft[i] = wholeMod(rowLeft[i] - mulValue * rowRight[i], domain);Нам нужно для всех возможных значений x = a - b * c в поле Z / N сохранить результат выражения x mod N. Воспользоваться периодичностью мы не сможем, потому что тогда для вычисления индекса снова придётся использовать деление по модулю. В итоге при 0 <= a, b, c < N получаем N + (N - 1)^2 возможных значений. Много, но деваться некуда.
Из этих значений (N - 1)^2 значений меньше 0. Поскольку отрицательные индексы невозможны, при индексировании значением a - b * c к нему нужно прибавить (N - 1)^2. Тогда функция для сложения строк модифицируется:
function mulSubRowCached(rowLeft, rowRight, mulValue, wholeModCache, cacheIndexOffset)
{
for(let i = 0; i < rowLeft.length; i++)
{
rowLeft[i] = wholeModCache[rowLeft[i] - mulValue * rowRight[i] + cacheIndexOffset];
}
}Заметим, что эта функция накладывает ограничение на mulValue — его значение не может быть больше domain и перед вызовом функции его тоже надо привести в наше поле Z / N. Кроме этого, обычное деление по модулю используется в функции mulRow().
Помимо wholeMod в вычитании строк матриц, используется . Кроме того, появилась вышеуказанная проблема с ограничением mulValue. Во всех этих случаях деление описывается формулой x = (a * b) mod N. Зная, что кэш хранит значения x = (c - a * b) mod N, мы можем вычислить (a * b) mod N, взяв значение кэша при c = 0 и вычтя его из N. Тогда функция для умножения строки на число модифицируется следующим образом:
function mulRowCached(row, mulValue, domain, wholeModCache, cacheIndexOffset)
{
for(let i = 0; i < row.length; i++)
{
row[i] = domain - wholeModCache[cacheIndexOffset - row[i] * mulValue];
}
}И получаем новое обращение матрицы:
function invertMatrix(matrix, domain)
{
let matrixSize = matrix.length;
//Инициализируем обратную матрицу единичной
let invMatrix = [];
for(let i = 0; i < matrixSize; i++)
{
let matrixRow = new Uint8Array(matrixSize);
matrixRow.fill(0);
matrixRow[i] = 1;
invMatrix.push(matrixRow);
}
//Вычисляем все обратные элементы заранее
let domainInvs = [];
for(let d = 0; d < domain; d++)
{
domainInvs.push(invModGcdEx(d, domain));
}
//Вычисляем кэш деления по модулю
const сacheIndexOffset = (domain - 1) * (domain - 1);
let wholeModCache = new Uint8Array((domain - 1) * (domain - 1) + domain);
for(let i = 0; i < wholeModCache.length; i++)
{
let divisor = i - сacheIndexOffset; //[-domainSizeCacheOffset, domainSize - 1]
wholeModCache[i] = wholeMod(divisor, domain); //Whole mod
}
//Прямой ход: приведение матрицы к ступенчатому виду
for(let i = 0; i < matrixSize; i++)
{
let thisRowFirst = matrix[i][i];
if(domainInvs[thisRowFirst] === 0) //Первый элемент строки 0 или делитель нуля, меняем строку местами со следующей строкой, у которой первый элемент не 0
{
for(let j = i + 1; j < matrixSize; j++)
{
let otherRowFirst = matrix[j][i];
if(domainInvs[thisRowFirst] !== 0) //Нашли строку с ненулевым первым элементом
{
thisRowFirst = otherRowFirst;
//Меняем строки местами
let tmpMatrixRow = matrix[i];
matrix[i] = matrix[j];
matrix[j] = tmpMatrixRow;
let tmpInvMatrixRow = invMatrix[i];
invMatrix[i] = invMatrix[j];
invMatrix[j] = tmpInvMatrixRow;
break;
}
}
}
//Обнуляем первые элементы всех строк после первой, отнимая от них (otherRowFirst / thisRowFirst) * x mod N
let invThisRowFirst = domainInvs[thisRowFirst]; // <---
for(let j = i + 1; j < matrixSize; j++)
{
let otherRowFirst = matrix[j][j];
if(domainInvs[otherRowFirst] !== 0)
{
let mulValue = domain - wholeModCache[сacheIndexOffset - otherRowFirst * invThisRowFirst]; // <---
mulSubRowCached(matrix[j], matrix[i], mulValue, wholeModCache, сacheIndexOffset); // <---
mulSubRowCached(invMatrix[j], invMatrix[i], mulValue, wholeModCache, сacheIndexOffset); // <---
}
}
}
//Обратный ход - обнуление всех элементов выше главной диагонали
let matrixRank = matrixSize;
for(let i = matrixSize - 1; i >= 0; i--)
{
let thisRowLast = matrix[i][i];
let invThisRowLast = domainInvs[thisRowLast];
for(let j = i - 1; j >= 0; j--)
{
let otherRowLast = matrix[j][i];
if(domainInvs[otherRowLast] !== 0)
{
let mulValue = domain - wholeModCache[сacheIndexOffset - otherRowLast * invThisRowLast]; // <---
mulSubRowCached(matrix[j], matrix[i], mulValue, wholeModCache, сacheIndexOffset); // <---
mulSubRowCached(invMatrix[j], invMatrix[i], mulValue, wholeModCache, сacheIndexOffset); // <---
}
}
if(domainInvs[thisRowLast] !== 0)
{
mulRowCached(matrix[i], invThisRowLast, domain, wholeModCache, сacheIndexOffset); // <---
mulRowCached(invMatrix[i], invThisRowLast, domain, wholeModCache, сacheIndexOffset); // <---
}
if(matrix[i].every(val => val === 0))
{
matrixRank -= 1;
}
}
return {inverse: invMatrix, rank: matrixRank};
}Замерим производительность. На той же матрице 500x500 по модулю 29 получаем время выполнения в ~5.4с.
Простите, что?
Нет, серьёзно, как это возможно? Кэшируем результат деления. Операции на два такта. В век супермедленной памяти и супербыстрых процессоров. Получаем прирост в 40%. Как?
Да, использование JavaScript создаёт определённый оверхед. Но JIT его нивелирует. Видимо, либо он нивелирует его недостаточно, либо не всё, чему нас учат про cache-friendly код — правда.
И да, размер кэша растёт квадратично. Но если сравнить среднее время в полях по разному модулю, то прирост будет не сильно отличаться:
В реальном проекте, где был применён этот метод, матрицы не рандомные и прирост ещё заметнее.
Заключение
Можно ли ещё больше ускорить вычисление? К сожалению, больше ни одного способа я не знаю. Я думал над распараллеливанием вычислений, но обращение матриц очень плохо параллелизуется. Поэтому пока остаётся так.
Полный код я выложил на Pastebin.
