We are taught that reading data from RAM is a terribly long operation. They give analogies to an office and a remote warehouse, forced to write cache-friendly code and instill a mortal fear of cache misses. We are also taught that processors are great at counting numbers, and it is often faster to compute the result twice than to store it in memory. It turns out that this is not always the case.
This article is based on a real project and real code that has been accelerated by the cache by almost one and a half times. All code is written in JavaScript.
A task
Let's say we have a matrix A of the order of 2000x2000. It is necessary to calculate its inverse matrix modulo N. In other words, it is necessary to find a matrix A -1 such that AA -1 mod N = E.
Since our calculations take place in the field modulo, iterative inversion methods will not work for us. We will use the good old Gauss method.
This post is dedicated to optimizing the Gaussian method for this particular case. In a real project, the inverse matrix is calculated in a separate WebWorker; in this example, we will use the main thread.
Secondary functions
For the program to work, we need four helper functions. The first is the calculation of (1 / x) mod N using the extended Euclidean algorithm:
The code
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 — .
, . , :
, , .
? , . , . .