Computing a matrix inverse is one of the most challenging problems in numerical programming. There are roughly a dozen major algorithms (LUP decomposition, QR decomposition, SVD decomposition, etc., etc.) and each algorithm has several variations (LUP: Doolittle, Crout, Gaussian elimination, etc.) The fact that there are so many algorithms for matrix inversion is a direct indication of the exceptional difficulty of the problem.

One interesting algorithm for estimating a matrix inverse (as opposed to computing with high precision) is Newton iteration. Briefly, one of several equivalent versions is:

set X(0)
loop
  X(k+1) = X(k) * (2I - A * X(k))
end-loop
return X(n)

The idea is to start with a guess answer X then loop, updating the current guess until some stopping criterion is reached. The A is the source matrix, X(k) is the previous guess, the X(k+1) is the next guess, I is the identity matrix.

Newton iteration, like most matrix inversion algorithms, can fail in many ways. In particular, the success or failure of Newton iteration completely depends on starting with a good initial guess X(0) otherwise the algorithm quickly blows up.

One way to initialize X(0) is 1/t * AT where AT is the transpose of the source matrix and t is the product of the largest row sum of absolute values and the largest column sum of absolute values. For example, if A =

  1   2   3   4
  8   7  -6   5
  0   2   6   4
  3   1   7   5

the largest row sum of absolute values is 8+7+6+5 = 26 and the largest column sum is 3+6+6+7 = 22, so t = 26 * 22.

I put together a demo using the C# language. The output of the demo is:

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  2.0  6.0  4.0
  3.0  1.0  7.0  5.0

Computing inverse
Done

Inverse:
   -0.416667    0.083333    0.000000    0.250000
   -0.675926    0.157407    0.722222   -0.194444
   -0.472222    0.027778    0.333333    0.083333
    1.046296   -0.120370   -0.611111   -0.027778

Computing A * Ainv =? I:
    1.000000    0.000000    0.000000    0.000000
    0.000000    1.000000    0.000000    0.000000
    0.000000    0.000000    1.000000    0.000000
    0.000000    0.000000    0.000000    1.000000

End demo

Although Newton iteration is conceptually simple, there are several tricky implementation details — the code below explains better than words. The code is long, but it’s not the length of the code (mostly simple helper functions) that is tricky, it’s the details. There is surprisingly little research on the effectiveness of using Newton iteration for computing a matrix inverse. I’ll do some experiments when I get a chance.

Left: Isaac Newton (1643-1727) was brilliant, but by most accounts, a very unpleasant person. Center: One of Newton’s bitter rivals was German mathematician Gottfried Leibniz (1846-1716). Newton and Leibnitz continuously accused each other of stealing their ideas, notably the invention of Calculus. Right: Actress Sarah Jessica Parker (b. 1965) is likely not an intellectual giant, but in the 1980s she had a similar hairstyle to Newton and Leibnitz.

Demo code. Replace “lt” (less than) with Boolean operator symbol (my lame blog editor often chokes on symbols).

using System;
using System.IO;

namespace MatrixInverseNewton
{
  internal class MatrixInverseNewtonProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin matrix inverse using" +
        " Newton iteration algorithm ");

      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, 2,  6, 4 };
      A[3] = new double[] { 3, 1,  7, 5 };

      Console.WriteLine("\nSource matrix: ");
      MatShow(A, 1, 5);

      Console.WriteLine("\nComputing inverse ");
      double[][] Ainv = MatInverseNewton(A);
      Console.WriteLine("Done ");

      Console.WriteLine("\nInverse: ");
      MatShow(Ainv, 6, 12);

      Console.WriteLine("\nComputing A * Ainv =? I: ");
      double[][] check = MatMult(A, Ainv);
      MatShow(check, 6, 12);

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main

    // ------------------------------------------------------

    static double[][] MatInverseNewton(double[][] A,
      int maxIter = 1000, double epsilon = 1.0e-6)
    {
      // Newton's method
      // X_k+1 = X_k * (2I - A*X_k)

      int n = A.Length; // must be square martix

      double[][] Xprev = MatStart(A);  // Pan algorithm
      double[][] Xnew = MatMake(n, n, 0.0);
      // double[][] Xnew = new double[0]][]; // dummy init OK
      double[][] I = MatIdentity(n);
      double[][] I2 = MatScalarMult(I, 2.0);

      int iter = 0;
      while (iter "lt" maxIter)
      {
        Xnew =
          MatProduct(Xprev, MatSubtract(I2,
          MatProduct(A, Xprev)));
        Xprev = MatCopyOf(Xnew);  // copy by val
        // Xnew = Xprev;  // copy by ref is OK

        if (iter % 10 == 0)
        {
          double[][] check = MatProduct(A, Xnew); // 
          if (MatAreEqual(check, I, epsilon) == true)
            return Xnew;
        }

        ++iter;
      } // while
      Console.WriteLine("WARN: no converge ");
      return Xnew;

      // ----------------------------------------------------
      // nested helpers: MatMake(), MatStart(), MatCopyOf(),
      // MatTranspose(), MatAreEqual(), MatProduct(),
      // MatSubtract(), MatIdentity(), MatScalarMult()
      // ----------------------------------------------------

      static double[][] MatMake(int nRows, int nCols,
        double v)
      {
        double[][] result = new double[nRows][];
        for (int i = 0; i "lt" nRows; ++i)
          result[i] = new double[nCols];
        for (int i = 0; i "lt" nRows; ++i)
          for (int j = 0; j "lt" nCols; ++j)
            result[i][j] = v;
        return result;
      }

      // ----------------------------------------------------

      static double[][] MatStart(double[][] m)
      {
        int n = m.Length;
        double maxRowSum = 0.0;
        double maxColSum = 0.0;
        for (int i = 0; i "lt" n; ++i)
        {
          double rowSum = 0.0;
          for (int j = 0; j "lt" n; ++j)
            rowSum += Math.Abs(m[i][j]);

          if (rowSum "gt" maxRowSum)
            maxRowSum = rowSum;
        }
        for (int j = 0; j "lt" n; ++j)
        {
          double colSum = 0.0;
          for (int i = 0; i "lt" n; ++i)
            colSum += Math.Abs(m[i][j]);
          if (colSum "gt" maxColSum)
            maxColSum = colSum;
        }

        double[][] result = MatTranspose(m);
        double t = 1.0 / (maxRowSum * maxColSum);
        for (int i = 0; i "lt" m.Length; ++i)
          for (int j = 0; j "lt" m.Length; ++j)
            result[i][j] *= t;
        return result;

      }

      // ----------------------------------------------------

      static double[][] MatTranspose(double[][] m)
      {
        int nr = m.Length; int nc = m[0].Length;
        double[][] result = MatMake(nc, nr, 0.0);  // note
        for (int i = 0; i "lt" nr; ++i)
          for (int j = 0; j "lt" nc; ++j)
            result[j][i] = m[i][j];
        return result;
      }

      // ------------------------------------------------------

      static double[][] MatCopyOf(double[][] m)
      {
        int nRows = m.Length; int nCols = m[0].Length;
        double[][] result = MatMake(nRows, nCols, 0.0);
        for (int i = 0; i "lt" nRows; ++i)
          for (int j = 0; j "lt" nCols; ++j)
            result[i][j] = m[i][j];
        return result;
      }

      // ------------------------------------------------------

      static bool MatAreEqual(double[][] matA,
       double[][] matB, double eps)
      {
        int nr = matA.Length; int nc = matB[0].Length;
        for (int i = 0; i "lt" nr; ++i)
          for (int j = 0; j "lt" nc; ++j)
            if (Math.Abs(matA[i][j] - matB[i][j]) "gt" eps)
              return false;
        return true;
      }

      // ----------------------------------------------------

      static double[][] MatProduct(double[][] matA,
      double[][] matB)
      {
        int aRows = matA.Length;
        int aCols = matA[0].Length;
        int bRows = matB.Length;
        int bCols = matB[0].Length;
        if (aCols != bRows)
          throw new Exception("Non-conformable matrices");

        double[][] result = MatMake(aRows, bCols, 0.0);

        for (int i = 0; i "lt" aRows; ++i) // each row of A
          for (int j = 0; j "lt" bCols; ++j) // each col of B
            for (int k = 0; k "lt" aCols; ++k)
              result[i][j] += matA[i][k] * matB[k][j];

        return result;
      }

      // ------------------------------------------------------

      static double[][] MatSubtract(double[][] matA,
        double[][] matB)
      {
        // matA - matB
        int nRows = matA.Length; int nCols = matA[0].Length;
        double[][] result = MatMake(nRows, nCols, 0.0);
        for (int i = 0; i "lt" nRows; ++i)
          for (int j = 0; j "lt" nCols; ++j)
            result[i][j] = matA[i][j] - matB[i][j];
        return result;
      }

      // ----------------------------------------------------

      static double[][] MatIdentity(int n)
      {
        double[][] result = MatMake(n, n, 0.0);
        for (int i = 0; i "lt" n; ++i)
          result[i][i] = 1.0;
        return result;
      }

      // ----------------------------------------------------

      static double[][] MatScalarMult(double[][] m, double u)
      {
        int nRows = m.Length; int nCols = m[0].Length;
        double[][] result = MatMake(nRows, nCols, 0.0);
        for (int i = 0; i "lt" nRows; ++i)
          for (int j = 0; j "lt" nCols; ++j)
            result[i][j] = u * m[i][j];
        return result;
      }

    } // MatInverseNewton()

    // ------------------------------------------------------
    // helpers for Main(): MatShow(), MatMult()
    // ------------------------------------------------------

    static void MatShow(double[][] m, int dec, int wid)
    {
      int nRows = m.Length; int nCols = m[0].Length;
      double small = 1.0 / Math.Pow(10, dec);
      for (int i = 0; i "lt" nRows; ++i)
      {
        for (int j = 0; j "lt" nCols; ++j)
        {
          double v = m[i][j];
          if (Math.Abs(v) "lt" small) v = 0.0;
          Console.Write(v.ToString("F" + dec).
            PadLeft(wid));
        }
        Console.WriteLine("");
      }
    }

    static double[][] MatMult(double[][] matA,
      double[][] matB)
    {
      int aRows = matA.Length;
      int aCols = matA[0].Length;
      int bRows = matB.Length;
      int bCols = matB[0].Length;
      if (aCols != bRows)
        throw new Exception("Non-conformable matrices");

      //double[][] result = MatMake(aRows, bCols, 0.0);
      double[][] result = new double[aRows][];
      for (int i = 0; i "lt" aRows; ++i)
        result[i] = new double[bCols];

      for (int i = 0; i "lt" aRows; ++i) // each row of A
        for (int j = 0; j "lt" bCols; ++j) // each col of B
          for (int k = 0; k "lt" aCols; ++k)
            result[i][j] += matA[i][k] * matB[k][j];

      return result;
    }

    // ------------------------------------------------------

  } // Program class

} // ns