I wrote an article titled “Matrix Inverse Using Newton Iteration with C#” in the May 2025 edition of Microsoft Visual Studio Magazine. See https://visualstudiomagazine.com/articles/2025/05/15/matrix-inverse-using-newton-iteration-with-c.aspx.

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. The fact that there are so many different techniques is a direct indication of how tricky it is to compute a matrix inverse.

My article presents a complete end-to-end demonstration of computing a matrix inverse using an algorithm called Newton iteration. Compared to other algorithms, Newton iteration is simple and easy to customize. The main disadvantage of Newton iteration is slower performance.

In ordinary arithmetic, 0.25 is the inverse of 4 because 0.25 * 4 = 1. In matrix algebra, the inverse of some matrix A is a matrix Ainv if A * Ainv = Ainv * A = I, where * indicates matrix multiplication, and I is the identity matrix (1.0 values on the upper-left to lower-right diagonal, 0.0 values elsewhere). Matrix inverse applies only to square matrices (equal number of rows and columns). Not all square matrices have an inverse.

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 periodically (perhaps once every 10 iterations) to check to see if X_k is close enough 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.

The demo program uses the Pan-Reif technique to set a starting matrix. In words, starting X_k 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.

The time required to invert a matrix using Newton iteration is dominated by calls to the MatMult() function which requires approximately (dim * dim * dim) operations. Therefore, time increases quickly and at some point around size 1000-by-1000 Newton iteration becomes impractical. Other matrix inverse algorithms can handle larger matrices, but for huge matrix sizes, special hardware and software are often needed.

I posted this blog about matrix inverse while I was visiting the city of Amsterdam in the Netherlands (sometimes aka Holland).

I visited the Max Euwe Museum in Amsterdam. It was wonderful, especially the huge chess library.