Dozens of machine learning algorithms require computing the inverse of a matrix. Computing a matrix inverse is conceptually easy, but implementation is one of the most difficult tasks in numerical programming.
The Wikipedia article on matrix inverse lists 14 different algorithms, and each algorithm has multiple variations, and each variation has dozens of implementation alternatives. Five examples include 1.) LUP decomposition algorithm (Crout version and others), 2.) QR decomposition algorithm (Householder version and others), 3.) SVD decomposition algorithm (Jacobi version and others), 4.) Cholesky decomposition (Banachiewicz version and others), 5.) Cayley-Hamilton algorithm (Faddeev–LeVerrier version and others). The fact that there are so many different techniques is a direct indication of how tricky it is to compute a matrix inverse.
One of my favorite techniques to compute a matrix inverse is called Newton iteration. It is simple, easy to implement, easy to modify, and easy to diagnose. One morning before work, I decided to brush off my C# implementation of Newton iteration and refactor it to JavaScript.
In pseudo-code, Newton iteration matrix inverse is:
let A equal the source matrix set X_k to an initial starting matrix loop max_iteration times set X_k+1 = X_k * (2I - (A * X_k)) set X_k = X_k+1 end-loop return X_k
Inside the processing loop, the * indicates matrix multiplication, and the – indicates matrix subtraction. The 2I indicates the scalar value 2 times the identity matrix, which is a matrix with 2.0 values on the diagonal, and 0.0 values elsewhere.
The X_k matrix slowly gets closer and closer to the inverse of the source matrix A. In most scenarios, you want the algorithm to periodically (perhaps once every 10 iterations) check to see if X_k is close enough (to within a small value epsilon) to the goal inverse, and early-exit if true.
Although it’s not obvious from the pseudo-code, the key to the Newton iteration is setting a good starting matrix X_k. An X_k with random values will usually not converge to the inverse of the source matrix.
My demo program uses the Pan-Reif technique to set a starting matrix. In words, the starting X_k matrix is (1/t) * A_tr where t is the product of the largest row sum of absolute values and the largest column sum of absolute values, and A_tr is the transpose of the source matrix A.
Suppose, as in the demo, source matrix A =
1.0 2.0 3.0 1.0 5.0 0.0 -5.0 4.0 1.0 4.0 6.0 1.0 0.0 -2.0 2.0 1.0 -4.0 5.0 3.0 2.0 0.0 2.0 4.0 0.0 -1.0
The largest row sum of absolute values is 15.0 from the fourth row. The largest column sum of absolute values is 16.0 from the third column. Therefore t = 15.0 * 16.0 = 240.0 and 1/t = 0.0042. The transpose of a matrix swaps rows and columns, so A_tr =
1.0 0.0 6.0 1.0 0.0 2.0 -5.0 1.0 -4.0 2.0 3.0 4.0 0.0 5.0 4.0 1.0 1.0 -2.0 3.0 0.0 5.0 4.0 2.0 2.0 -1.0
And therefore the starting X_k = (1/t) * A_tr =
0.0042 0.0000 0.0250 0.0042 0.0000 0.0083 -0.0208 0.0042 -0.0167 0.0083 0.0125 0.0167 0.0000 0.0208 0.0167 0.0042 0.0042 -0.0083 0.0125 0.0000 0.0208 0.0167 0.0083 0.0083 -0.0042
The main challenge when implementing Newton iteration matrix inverse is coding up a lot of small but relatively simple helper functions.
The output of my demo program is:
Begin matrix inverse using Newton iteration algorithm with JavaScript Source matrix: 1.0 2.0 3.0 1.0 5.0 0.0 -5.0 4.0 1.0 4.0 6.0 1.0 0.0 -2.0 2.0 1.0 -4.0 5.0 3.0 2.0 0.0 2.0 4.0 0.0 -1.0 Setting maxIter = 1000 Setting eps = 1.0e-8 Setting verbose = true Computing inverse Done. Required 10 iterations. Inverse: -0.0316 -0.1190 0.1472 0.1483 -0.0428 0.1227 -0.1264 -0.0186 -0.0112 0.0483 -0.0242 0.0855 0.0067 -0.0160 0.2026 0.1152 -0.3309 -0.0781 0.3532 -0.1970 0.1487 0.0892 -0.0104 -0.0862 -0.0929 Checking A * Ainv = I 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 Check PASS End demo
The eps (epsilon) value controls how close to the correct inverse the result is before short-circuit exiting the processing iteration loop. The demo program checks if the computed inverse times the source matrix equals the identity matrix using a matAreEqual() function.
A good way to spend some time one morning before work.
Unlike most of my colleagues, I really enjoy coding with JavaScript. JavaScript is a language where you have to be very, very careful because there’s a lot that could go wrong. Here are two signs where the “What could possibly go wrong?” question is implied.
Demo code. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols (my blog editor consistently chokes on symbols).
// matrix_inverse_newton_iteration.js
// node.js
let FS = require('fs'); // for loadTxt()
// ----------------------------------------------------------
function matInverseNewton(A, maxIter, eps, verbose)
{
// Newton's method
// X_k+1 = X_k * (2I - A*X_k)
let n = A.length; // must be square martix
let Xcurr = matStart(A); // Pan-Reif algorithm
let Xnew = matMake(n, n, 0.0);
let I = matIdentity(n);
let I2 = matScalarMult(2.0, I);
let iter = 0;
while (iter "lt" maxIter) {
Xnew =
matMult(Xcurr, matSubtract(I2, matMult(A, Xcurr)));
// Xcurr = matCopyOf(Xnew); // by value
Xcurr = Xnew; // by ref
if (iter % 10 == 0) {
let check = matMult(A, Xcurr); //
if (matAreEqual(check, I, eps) == true) {
if (verbose == true) {
console.log("Done. Required " +
iter.toString() + " iterations. ");
}
return Xcurr;
}
}
++iter;
} // while
console.log("WARN: no converge ");
return Xnew;
}
// ----------------------------------------------------------
function matStart(A)
{
// Pan-Reif initial matrix for Newton iteration inverse
let n = A.length;
let maxRowSum = 0.0;
let maxColSum = 0.0;
for (let i = 0; i "lt" n; ++i) {
let rowSum = 0.0;
for (let j = 0; j "lt" n; ++j) {
rowSum += Math.abs(A[i][j]);
}
if (rowSum "gt" maxRowSum) {
maxRowSum = rowSum;
}
}
for (let j = 0; j "lt" n; ++j) {
let colSum = 0.0;
for (let i = 0; i "lt" n; ++i) {
colSum += Math.abs(A[i][j]);
}
if (colSum "gt" maxColSum) {
maxColSum = colSum;
}
}
let result = matTranspose(A);
let t = 1.0 / (maxRowSum * maxColSum);
for (let i = 0; i "lt" n; ++i) {
for (let j = 0; j "lt" n; ++j) {
result[i][j] *= t;
}
}
return result;
}
// ----------------------------------------------------------
function loadTxt(fn, delimit, usecols, comment) {
// efficient but mildly complicated
let all = FS.readFileSync(fn, "utf8"); // giant string
all = all.trim(); // strip final crlf in file
let lines = all.split("\n"); // array of lines
// count number non-comment lines
let nRows = 0;
for (let i = 0; i "lt" lines.length; ++i) {
if (!lines[i].startsWith(comment))
++nRows;
}
let nCols = usecols.length;
let result = matMake(nRows, nCols, 0.0);
let r = 0; // into lines
let i = 0; // into result[][]
while (r "lt" lines.length) {
if (lines[r].startsWith(comment)) {
++r; // next row
}
else {
let tokens = lines[r].split(delimit);
for (let j = 0; j "lt" nCols; ++j) {
result[i][j] = parseFloat(tokens[usecols[j]]);
}
++r;
++i;
}
}
return result;
}
// ----------------------------------------------------------
function vecMake(n, val)
{
let result = [];
for (let i = 0; i "lt" n; ++i) {
result[i] = val;
}
return result;
}
// ----------------------------------------------------------
function matMake(nRows, nCols, val)
{
let result = [];
for (let i = 0; i "lt" nRows; ++i) {
result[i] = [];
for (let j = 0; j "lt" nCols; ++j) {
result[i][j] = val;
}
}
return result;
}
// ----------------------------------------------------------
function vecShow(vec, dec, wid, nl)
{
for (let i = 0; i "lt" vec.length; ++i) {
let x = vec[i];
if (Math.abs(x) "lt" 1.0/Math.pow(10,dec)) {
x = 0.0; // avoid "-0.00"
}
let xx = x.toFixed(dec);
let s = xx.toString().padStart(wid, ' ');
process.stdout.write(s);
process.stdout.write(" ");
}
if (nl == true)
process.stdout.write("\n");
}
// ----------------------------------------------------------
function matShow(A, dec, wid)
{
let nRows = A.length;
let nCols = A[0].length;
for (let i = 0; i "lt" nRows; ++i) {
for (let j = 0; j "lt" nCols; ++j) {
let x = A[i][j];
if (Math.abs(x) "lt" 1.0/Math.pow(10,dec)) {
x = 0.0; // avoid "-0.00"
}
if (x "gte" 0.0) {
process.stdout.write(" "); // + or - space
}
process.stdout.write(x.toFixed(dec));
process.stdout.write(" ");
}
process.stdout.write("\n");
}
}
// ----------------------------------------------------------
function matMult(A, B)
{
let aRows = A.length;
let aCols = A[0].length;
let bRows = B.length;
let bCols = B[0].length;
if (aCols != bRows)
console.log("** FATAL ** Non-conformable matrices");
let result = matMake(aRows, bCols, 0.0);
for (let i = 0; i "lt" aRows; ++i) { // each row of A
for (let j = 0; j "lt" bCols; ++j) { // each col of B
for (let k = 0; k "lt" aCols; ++k) {
result[i][j] += A[i][k] * B[k][j];
}
}
}
return result;
}
// ----------------------------------------------------------
function matAdd(A, B)
{
let nRows = A.length;
let nCols = A[0].length;
let result = matMake(nRows, nCols, 0.0);
for (let i = 0; i "lt" nRows; ++i) {
for (let j = 0; j "lt" nCols; ++j) {
result[i][j] = A[i][j] + B[i][j];
}
}
return result;
}
// ----------------------------------------------------------
function matSubtract(A, B)
{
let nRows = A.length;
let nCols = A[0].length;
let result = matMake(nRows, nCols, 0.0);
for (let i = 0; i "lt" nRows; ++i) {
for (let j = 0; j "lt" nCols; ++j) {
result[i][j] = A[i][j] - B[i][j];
}
}
return result;
}
// ----------------------------------------------------------
function matTrace(A)
{
let result = 0.0;
for (let i = 0; i "lt" A.length; ++i) {
result += A[i][i];
}
return result;
}
// ----------------------------------------------------------
function matCopyOf(A)
{
let nRows = A.length;
let nCols = A[0].length;
let result = matMake(nRows, nCols, 0.0);
for (let i = 0; i "lt" nRows; ++i) {
for (let j = 0; j "lt" nCols; ++j) {
result[i][j] = A[i][j];
}
}
return result;
}
// ----------------------------------------------------------
function matScalarMult(u, A)
{
let nRows = A.length;
let nCols = A[0].length;
let result = matMake(nRows, nCols, 0.0);
for (let i = 0; i "lt" nRows; ++i) {
for (let j = 0; j "lt" nCols; ++j) {
result[i][j] = u * A[i][j];
}
}
return result;
}
// ----------------------------------------------------------
function matTranspose(A)
{
let nRows = A.length;
let nCols = A[0].length;
let result = matMake(nCols, nRows, 0.0); // note reverse
for (let i = 0; i "lt" nRows; ++i) {
for (let j = 0; j "lt" nCols; ++j) {
result[j][i] = A[i][j];
}
}
return result;
}
// ----------------------------------------------------------
function matIdentity(n)
{
let result = matMake(n, n, 0.0);
for (let i = 0; i "lt" n; ++i)
result[i][i] = 1.0;
return result;
}
// ----------------------------------------------------------
function matAreEqual(A, B, eps)
{
let nr = A.length;
let nc = B[0].length;
for (let i = 0; i "lt" nr; ++i)
for (let j = 0; j "lt" nc; ++j)
if (Math.abs(A[i][j] - B[i][j]) "gt" eps)
return false;
return true;
}
// ----------------------------------------------------------
function main()
{
console.log("\nBegin matrix inverse using " +
"Newton iteration algorithm with JavaScript ");
let A = [];
A[0] = [ 1.0, 2.0, 3.0, 1.0, 5.0 ];
A[1] = [ 0.0, -5.0, 4.0, 1.0, 4.0 ];
A[2] = [ 6.0, 1.0, 0.0, -2.0, 2.0 ];
A[3] = [ 1.0, -4.0, 5.0, 3.0, 2.0 ];
A[4] = [ 0.0, 2.0, 4.0, 0.0, -1.0 ];
// A = loadTxt(".\\dummy_data.txt",
// ",", [0,1,2,3,4], "#");
console.log("\nSource matrix: ");
matShow(A, 1, 5);
let maxIter = 1000;
let eps = 1.0e-8;
let verbose = true;
console.log("\nSetting maxIter = " +
maxIter.toString());
console.log("Setting eps = " +
eps.toExponential(1).toString());
console.log("Setting verbose = " +
verbose.toString());
console.log("\nComputing inverse ");
let Ainv = matInverseNewton(A, maxIter, eps, verbose);
console.log("\nInverse: ");
matShow(Ainv, 4, 9);
console.log("\nChecking A * Ainv = I ");
let check = matMult(A, Ainv);
matShow(check, 4, 9);
let I = matIdentity(A.length);
if (matAreEqual(check, I, eps) == true)
console.log("\nCheck PASS ");
else
console.log("\nCheck FAIL ");
console.log("\nEnd demo ");
}
main()
.NET Test Automation Recipes
Software Testing
SciPy Programming Succinctly
Keras Succinctly
R Programming
2026 Visual Studio Live
2025 Summer MLADS Conference
2025 DevIntersection Conference
2025 Machine Learning Week
2025 Ai4 Conference
2025 G2E Conference
2025 iSC West Conference
You must be logged in to post a comment.