On nous apprend que la lecture de données depuis la RAM est une opération terriblement longue. Ils donnent des analogies avec un bureau et un entrepôt distant, obligés d'écrire du code compatible avec le cache et de susciter une peur mortelle des échecs de cache. On nous apprend également que les processeurs sont excellents pour compter les nombres et qu'il est souvent plus rapide de calculer le résultat deux fois que de le stocker en mémoire. Il s'avère que ce n'est pas toujours le cas.
Cet article est basé sur un vrai projet et du vrai code qui a été accéléré par le cache de près d'une fois et demie. Tout le code est écrit en JavaScript.
Une tâche
Disons que nous avons une matrice A de l'ordre de 2000x2000. Il faut calculer sa matrice inverse modulo N. En d'autres termes, il faut trouver une matrice A -1 telle que AA -1 mod N = E.
Puisque nos calculs ont lieu dans le champ modulo, les méthodes d'inversion itératives ne fonctionneront pas pour nous. Nous utiliserons la bonne vieille méthode de Gauss.
Cet article est dédié à l'optimisation de la méthode gaussienne pour ce cas particulier. Dans un projet réel, le calcul de la matrice inverse a lieu dans un WebWorker séparé, dans cet exemple nous allons gérer avec le thread principal.
Fonctions secondaires
Pour que le programme fonctionne, nous avons besoin de quatre fonctions d'assistance. Le premier est le calcul de (1 / x) mod N en utilisant l'algorithme euclidien étendu:
Le 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 — .
, . , :
, , .
? , . , . .