Kami diajari bahwa membaca data dari RAM adalah operasi yang sangat lama. Mereka memberikan analogi ke kantor dan gudang terpencil, dipaksa untuk menulis kode ramah-cache dan menanamkan ketakutan yang mematikan akan kehilangan cache. Kita juga diajari bahwa prosesor hebat dalam menghitung angka, dan seringkali lebih cepat untuk menghitung hasilnya dua kali daripada menyimpannya dalam memori. Ternyata tidak selalu demikian.
Artikel ini didasarkan pada proyek nyata dan kode nyata yang telah dipercepat oleh cache hampir satu setengah kali. Semua kode ditulis dalam JavaScript.
Sebuah tugas
Misalkan kita memiliki matriks A dengan orde 2000x2000. Untuk itu diperlukan perhitungan matriks invers modulo N. Dengan kata lain, matriks A -1 perlu dicari sedemikian rupa sehingga AA -1 mod N = E.
Karena kalkulasi kita dilakukan di field modulo, metode inversi iteratif tidak akan bekerja untuk kita. Kami akan menggunakan metode Gauss lama yang baik.
Posting ini didedikasikan untuk mengoptimalkan metode Gaussian untuk kasus khusus ini. Dalam proyek nyata, perhitungan matriks terbalik dilakukan di WebWorker terpisah, dalam contoh ini kita akan mengelola dengan utas utama.
Fungsi sekunder
Agar program dapat bekerja, kita membutuhkan empat fungsi pembantu. Yang pertama adalah perhitungan (1 / x) mod N menggunakan algoritma extended Euclidean:
Kode
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 — .
, . , :
, , .
? , . , . .