Uns wird beigebracht, dass das Lesen von Daten aus dem RAM eine furchtbar lange Operation ist. Sie geben Analogien zu einem BĂŒro und einem Remote-Lagerhaus, die gezwungen sind, cachefreundlichen Code zu schreiben und eine tödliche Angst vor Cache-Fehlern zu wecken. Es wird uns auch beigebracht, dass Prozessoren Zahlen hervorragend zĂ€hlen können und es oft schneller ist, das Ergebnis zweimal zu berechnen, als es im Speicher zu speichern. Es stellt sich heraus, dass dies nicht immer der Fall ist.
Dieser Artikel basiert auf einem realen Projekt und echtem Code, der durch den Cache um fast das Eineinhalbfache beschleunigt wurde. Der gesamte Code ist in JavaScript geschrieben.
Eine Aufgabe
Angenommen, wir haben eine Matrix A in der GröĂenordnung von 2000 x 2000. Es ist notwendig, seine inverse Matrix Modulo N zu berechnen. Mit anderen Worten, es ist notwendig, eine Matrix A -1 zu finden, so dass AA -1 Mod N = E ist.
Da unsere Berechnungen im Feld Modulo stattfinden, funktionieren iterative Inversionsmethoden fĂŒr uns nicht. Wir werden die gute alte GauĂ-Methode verwenden.
Dieser Beitrag widmet sich der Optimierung der GauĂschen Methode fĂŒr diesen speziellen Fall. In einem realen Projekt erfolgt die Berechnung der inversen Matrix in einem separaten WebWorker. In diesem Beispiel werden wir mit dem Hauptthread verwalten.
SekundÀrfunktionen
Damit das Programm funktioniert, benötigen wir vier Hilfsfunktionen. Das erste ist die Berechnung von (1 / x) mod N unter Verwendung des erweiterten euklidischen Algorithmus:
Der 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 â .
, . , :
, , .
? , . , . .