The Data Science Lab
Matrix Inverse Using Newton Iteration with C#
Dr. James McCaffrey from Microsoft Research presents a complete end-to-end demonstration of computing a matrix inverse using the Newton iteration algorithm. Compared to other algorithms, Newton iteration is simple and easy to customize, but the technique is relatively slow.
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. Four examples include 1.) LUP decomposition algorithm (Crout version), 2.) QR decomposition algorithm (Householder version), 3.) SVD decomposition algorithm (Jacobi version), 4.) Cholesky decomposition (Banachiewicz version). The fact that there are so many different techniques is a direct indication of how tricky it is to compute a matrix inverse.
This 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 very simple and easy to customize, which is often necessary in practice. 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.
[Click on image for larger view.] Figure 1: Computing a Matrix Inverse Using Newton Iteration
A good way to see where this article is headed is to take a look at the screenshot of a demo program in Figure 1. The demo begins by setting up a 4-by-4 matrix:
Begin matrix inverse using Newton iteration algorithm
Source matrix:
1.0 -2.0 3.0 4.0
8.0 7.0 -6.0 5.0
0.0 -5.0 1.0 9.0
3.0 1.0 -7.0 5.0
Next, the demo computes and displays the inverse of the source matrix:
Setting tolerance = 0.00000001
Setting maxIter = 1000
Computing inverse
Done
Inverse:
1.30000000 -0.30000000 -0.80000000 0.70000000
-1.07941176 0.40294118 0.65294118 -0.71470588
-0.02352941 0.08235294 0.08235294 -0.21176471
-0.59705882 0.21470588 0.46470588 -0.37352941
When working with matrix inverses, it's necessary to set a tolerance value which is roughly the number of decimal places of accuracy. The demo uses 1.0e-8 which is 0.00000001 (8 decimal places), which is a typical tolerance. The Newton iteration algorithm is iterative and requires a maximum number of iterations. The demo sets the maxIter parameter to 1,000. The maxIter parameter must be determined by trial and error.
The demo program concludes by verifying that the computed inverse is correct:
Checking A * Ainv == I
1.00000000 0.00000000 0.00000000 0.00000000
0.00000000 1.00000000 0.00000000 0.00000000
0.00000000 0.00000000 1.00000000 0.00000000
0.00000000 0.00000000 0.00000000 1.00000000
Verified A * Ainv == I
End demo
Behind the scenes, the demo program calls a function MatAreEqual() that compares A * Ainv with the identity matrix to check that the two matrices are the same to within the 1.0e-8 tolerance.
This article assumes you have intermediate or better programming skill but doesn't assume you know anything about the Newton iteration algorithm for computing the inverse of a matrix. The demo is implemented using C# but you should be able to refactor the demo code to another C-family language if you wish. Most normal error checking has been removed to keep the main ideas as clear as possible.
The source code for the demo program is a bit too long to be presented in its entirety in this article. The complete code is available in the accompanying file download, and is also available online.
Understanding the Newton Iteration Algorithm
The Newton iteration algorithm is deceptively simple looking. In pseudo-code:
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.
Suppose, as in the demo, source matrix A =
1.0 -2.0 3.0 4.0
8.0 7.0 -6.0 5.0
0.0 -5.0 1.0 9.0
3.0 1.0 -7.0 5.0
The largest row sum of absolute values is 26.0 from the second row. The largest column sum of absolute values is 23.0 from the last column. Therefore t = 26.0 * 23.0 = 598.0 and 1/t = 0.00167224. The transpose of a matrix swaps rows and columns, so A_tr =
1.0 8.0 0.0 3.0
-2.0 7.0 -5.0 1.0
3.0 -6.0 1.0 -7.0
4.0 5.0 9.0 5.0
And therefore the starting X_k = (1/t) * A_tr =
0.00167224 0.01337793 0.00000000 0.00501672
-0.00334448 0.01170569 -0.00836120 0.00167224
0.00501672 -0.01003344 0.00167224 -0.01170569
0.00668896 0.00836120 0.01505017 0.00836120
Notice that if the source matrix A has all 0.0 values, t = 0.0 and 1/t does not exist.
The Demo Program
I used Visual Studio 2022 (Community Free Edition) for the demo program. I created a new C# console application and checked the "Place solution and project in the same directory" option. I specified .NET version 8.0. I named the project MatrixInverseNewtonIteration. I checked the "Do not use top-level statements" option to avoid the weird program entry point shortcut syntax.
The demo has no significant .NET dependencies and any relatively recent version of Visual Studio with .NET (Core) or the older .NET Framework will work fine. You can also use the Visual Studio Code program if you like.
After the template code loaded into the editor, I right-clicked on file Program.cs in the Solution Explorer window and renamed the file to the more descriptive MatrixInverseNewtonProgram.cs. I allowed Visual Studio to automatically rename class Program.
The overall program structure is presented in Listing 1. All the control logic is in the Main() method in the Program class. All of the matrix inverse functionality is contained in static functions/methods in a Utils class.
Listing 1: Overall Program Structure
using System;
using System.IO;
using System.Collections.Generic;
namespace MatrixInverseNewtonIteration
{
internal class MatrixInverseNewtonProgram
{
static void Main(string[] args)
{
Console.WriteLine("Begin matrix inverse using" +
" Newton iteration algorithm ");
// 1. set up and display source matrix
// 2. compute and display inverse matrix
// 3. verify inverse is correct
Console.WriteLine("End demo ");
Console.ReadLine();
}
} // class Program
public class Utils
{
public static double[][] MatMake(int nRows,
int nCols) { . . }
public static void MatShow(double[][] A,
int dec, int wid) { . . }
public static double[][] MatLoad(string fn,
int[] usecols, char sep, string comment) { . . }
public static double[][] MatInverseNewton(double[][] A,
int maxIter = 1000, double eps = 1.0e-8) { . . }
public static double[][] MatStart(double[][] A) { . . }
public static double[][] MatMult(double[][] A,
double[][] B) { . . }
public static double[][] MatSubtract(double[][] A,
double[][] B) { . . }
public static double[][] MatCopyOf(double[][] A) { . . }
public static double[][] MatScalarMult(double[][] A,
double u) { . . }
public static double[][] MatTranspose(double[][] A) { . . }
public static double[][] MatIdentity(int n) { . . }
public static bool MatAreEqual(double[][] A,
double[][] B, double eps) { . . }
}
} // ns
The demo program begins by setting up a source matrix using these statements:
double[][] A = new double[4][];
A[0] = new double[] { 1, -2, 3, 4 };
A[1] = new double[] { 8, 7, -6, 5 };
A[2] = new double[] { 0, -5, 1, 9 };
A[3] = new double[] { 3, 1, -7, 5 };
Console.WriteLine("Source matrix: ");
Utils.MatShow(A, 1, 6); // 1 decimal, 6 cols wide
The demo defines a matrix as an array of arrays of type double. It's possible to use a program-defined matrix class, but in my opinion this approach creates unnecessary complexity. The most recent versions of .NET allow you to use a shortcut syntax to initialize an array, for example, A[0] = [1, -2, 3, 4] but using this syntax creates a version dependency.
You access each cell of an array-of-arrays style using syntax like A[0][3] = 1.23 where the first index is the row and the second index is the column. Unlike most programming languages, C# also has a 2D array type that is accessed like A[0,3] = 1.23 but I do not recommend using this data type mostly because all published matrix algorithms assume an array-of-arrays style matrix.
In a non-demo scenario, you might want to load data into a matrix from a text file. The demo code has a MatLoad() function that isn't used but could be called like:
string fn = "dummy_data.txt";
double[][] A = Utils.MatLoad(fn, new int[] { 0, 1, 2, 3 },
',', "#");
This syntax means load the data in file dummy_data.txt, using columns [0] to [3], where the data is delimited by the comma character, and comments to be ignored are lines that begin with the # character.
Computing the Matrix Inverse
The demo program computes and displays the inverse of the source matrix using these statements:
double tol = 1.0e-8;
int maxIer = 1000;
Console.WriteLine("Computing inverse ");
double[][] Ainv = Utils.MatInverseNewton(A, 1000, tol);
Console.WriteLine("Done ");
Console.WriteLine("Inverse: ");
Utils.MatShow(Ainv, 8, 12);
For matrices that hold type double, precision to 6 or 8 decimal places is common. The maximum number of iterations varies a lot depending on the size and complexity of the source matrix, but 1,000 usually works well in practice as a starting point.
The code for the Utils.MatInverseNewton() function is shown in Listing 2. The code presented in this article is explicitly intended for you to modify to meet your needs. For example, you might want to verify that the source matrix A has the same number of rows and columns.
Listing 2: The MatInverseNewton Function
public static double[][] MatInverseNewton(double[][] A,
int maxIter = 1000, double eps = 1.0e-8)
{
// X_k+1 = X_k * (2I - A*X_k)
int n = A.Length; // must be square martix
double[][] Xcurr = MatStart(A); // Pan-Reif algorithm
double[][] Xnew = MatMake(n, n);
double[][] I = MatIdentity(n);
double[][] I2 = MatScalarMult(I, 2.0);
int iter = 0;
while (iter < maxIter)
{
Xnew =
MatMult(Xcurr, MatSubtract(I2,
MatMult(A, Xcurr)));
// Xcurr = MatCopyOf(Xnew); // by value
Xcurr = Xnew; // by ref
if (iter % 10 == 0)
{
double[][] check = MatMult(A, Xcurr);
if (MatAreEqual(check, I, eps) == true)
return Xcurr;
}
++iter;
} // while
Console.WriteLine("WARN: no converge ");
return Xnew;
}
The code for the MatInverseNewton() function is short because most of the work is performed by helper functions MatStart(), MatMake(), MatIdentity(), MatScalrMult(), MatMult(), MatSubtract(), and MatAreEqual().
The demo implementation checks every 10 iterations to see if the current solution matrix is close enough, to within a small epsilon value, to the goal inverse matrix. Checking is a relatively expensive operation so in most situations you don't want to check after every iteration.
Not all matrices have an inverse. For example, a matrix that has a pair of identical rows, or a pair of identical columns, does not have a matrix. In such cases the demo implementation will not fail, but instead the demo inverse function will return a result matrix where each cell is double.NaN (not a number). You can check for this before passing the source matrix to the MatInverseNewton() function, or inside the function implementation, or you can check the return result.
If the solution matrix does not converge to a result, the demo prints a warning message. An alternative design is to throw an Exception. Another modification you can make is to programmatically increase the maximum number of iterations.
Checking the Matrix Inverse
The demo concludes by checking that the computed matrix inverse is correct:
Console.WriteLine("Checking A * Ainv == I ");
double[][] check = Utils.MatMult(A, Ainv);
Utils.MatShow(check, 8, 12);
double[][] I = Utils.MatIdentity(4); // 4x4
if (Utils.MatAreEqual(check, I, tol) == true)
Console.WriteLine("Verified A * Ainv == I ");
else
Console.WriteLine("Fail A * Ainv == I ");
Recall that for a source matrix A, the definition of a matrix inverse Ainv is one where A * Ainv = Ainv * A = I. The demo program checks just A * Ainv. In a non-demo scenario, it'd be more thorough to also check Ainv * A.
Wrapping Up
Compared to other matrix inverse techniques, the Newton iteration algorithm is easy to modify and is highly interpretable if you insert diagnostic print statements. The main disadvantage of Newton iteration is speed performance. Here are some timing results for inverting dim-by-dim matrices filled with random values between -1.0 and +1.0, using a tolerance of 1.0e-8, run on a standard laptop:
dim time
--------------------
100 0.3 seconds
200 3.1 seconds
300 13.8 seconds
400 33.6 seconds
500 85.1 seconds
600 134.1 seconds ( ~2 minutes)
700 238.8 seconds ( ~4 minutes)
800 359.6 seconds ( ~6 minutes)
900 755.2 seconds (~12 minutes)
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.