Нас учат, что чтение данных из оперативной памяти — ужасно долгая операция. Приводят аналогии с офисом и удалённым складом, заставляют писать 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(). Во всех этих случаях деление описывается формулой 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.
