I like to write code almost every day, to entertain myself, and to keep in practice. One morning before work, I realized it had been many, many months since I implemented a function to compute the inverse of a matrix using the Gauss-Jordan elimination to reduced row echelon form.

The Wikipedia article on matrix inverse lists 14 different algorithms, and each algorithm has multiple variations, and each variation has dozens of implementation alternatives. Six 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), 6.) Newton Iteration. The fact that there are so many different techniques is a direct indication of how tricky it is to compute a matrix inverse.

The Gauss-Jordan technique is arguably the simplest way to compute a matrix inverse. But Gauss-Jordan is rarely used in practice because, compared to other techniques, Gauss-Jordan is susceptible to numerical instability, especially when the source matrix is nearly singular (meaning has no inverse). Gauss-Jordan inverse is mostly a teaching technique for students.

Anyway, after a few minutes of coding, I had a C# language demo up and running. The output is:

Matrix inverse using Gauss-Jordan elimination 
  to reduced row echelon form

Source matrix:
  3.0  5.0  9.0
  1.0  1.0  5.0
  2.0  4.0  8.0

Computing inverse
Done

Inverse from Gauss-Jordan:
   1.5000   0.5000  -2.0000
  -0.2500  -0.7500   0.7500
  -0.2500   0.2500   0.2500

Minv * M:
   1.0000  -0.0000  -0.0000
   0.0000   1.0000   0.0000
   0.0000   0.0000   1.0000

M * Minv:
   1.0000  -0.0000   0.0000
  -0.0000   1.0000   0.0000
  -0.0000  -0.0000   1.0000

End demo

The Gauss-Jordan technique starts by creating a matrix that is the source matrix that is augmented by the Identity matrix to the right. For the demo matrix, the augmented matrix is:

 3  5  9   1  0  0
 1  1  5   0  1  0
 2  4  8   0  0  1

Then the algorithm iterates through each row and performs three possible row operations:

* exchange two rows
* multiply a row by a constant
* multiply a row by a constant times another row

The version I implemented doesn’t exchange rows, but some versions do. The goal is to transform the augmented matrix so that the left side of the augmented matrix becomes the identity matrix. For the demo, the result is:

 1  0  0    1.50   0.50  -2.00
 0  1  0   -0.25  -0.75   0.75  
 0  0  1   -0.25   0.25   0.25

When finished, the right side of the transformed augmented matrix is the inverse of the source matrix. Clever.

OK, much fun for me. But as mentioned, the Gauss-Jordan algorithm to compute a matrix inverse is not used very often, so this was just mental exercise for me.

I loved building model airplanes when I was a young man. Two of my favorite models were World War II aircraft that featured inverted gull wing designs.

Right: The U.S. Vought F4U Corsair was a fighter bomber that entered service in mid-war. It was a rare example of a plane that was equally effective as a fighter and as an attack bomber. It was one of the fastest planes of the war with a top speed of 440 mph. The inverted gull wing design allowed a huge 13 foot 4 inch propeller.

Demo code. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols. My blog editor chokes on symbols.

using System;
using System.Collections.Generic;
using System.IO;

namespace MatrixInverseGaussJordan
{
  internal class MatrixInverseGaussJordanProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nMatrix inverse " +
        "using Gauss-Jordan elimination to reduced " +
        "row echelon form ");

      double[][] M = new double[3][];
      M[0] = new double[] { 3, 5, 9 };
      M[1] = new double[] { 1, 1, 5 };
      M[2] = new double[] { 2, 4, 8 };


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

      Console.WriteLine("\nComputing inverse ");
      double[][] Minv = MatInverseGaussJordan(M);
      Console.WriteLine("Done ");

      Console.WriteLine("\nInverse from Gauss-Jordan: ");
      MatShow(Minv, 4, 9);

      double[][] Minv_M = MatMult(Minv, M);
      Console.WriteLine("\nMinv * M: ");
      MatShow(Minv_M, 4, 9);

      double[][] M_Minv = MatMult(M, Minv);
      Console.WriteLine("\nM * Minv: ");
      MatShow(M_Minv, 4, 9);

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

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

    static double[][] MatInverseGaussJordan(double[][] A)
    {
      int n = A.Length;  // assume A is square
      double[][] B = MatMake(n, 2 * n);  // make augmented
      for (int i = 0; i "lt" n; ++i)
      {
        for (int j = 0; j "lt" n; ++j)
        {
          B[i][j] = A[i][j];
        }
        B[i][i + n] = 1.0;
      }

      for (int i = 0; i "lt" n; ++i)
      {
        double pv = B[i][i];
        if (Math.Abs(pv) "lt" 1.0e-12)
          Console.WriteLine("FATAL: matrix is singular ");

        for (int j = 0; j "lt" 2 * n; ++j)
          B[i][j] /= pv;  // make diag element 1.0
        for (int j = 0; j "lt" n; ++j)  // elim other rows
        {
          if (i != j)
          {
            double factor = B[j][i];
            for (int k = 0; k "lt" 2 * n; ++k)
              B[j][k] -= factor * B[i][k];  // tricky
          }
        }
      } // i

      // extract the inverse from augmented
      double[][] result = MatMake(n, n);
      for (int i = 0; i "lt" n; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = B[i][j + n];

      return result;
    }

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

    // helper functions:
    // MatMake(), MatShow(), MatMult(), MatLoad()

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

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

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

    static void MatShow(double[][] M, int dec, int wid)
    {
      for (int i = 0; i "lt" M.Length; ++i)
      {
        for (int j = 0; j "lt" M[0].Length; ++j)
        {
          double v = M[i][j];
          Console.Write(v.ToString("F" + dec).
            PadLeft(wid));
        }
        Console.WriteLine("");
      }
    }

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

    static double[][] MatLoad(string fn, int[] usecols,
      char sep, string comment)
    {
      // load matrix from text file
      List"lt"double[]"gt" result =
        new List"lt"double[]"gt"();
      string line = "";
      FileStream ifs = new FileStream(fn, FileMode.Open);
      StreamReader sr = new StreamReader(ifs);
      while ((line = sr.ReadLine()) != null)
      {
        if (line.StartsWith(comment) == true)
          continue;
        string[] tokens = line.Split(sep);
        List"lt"double"gt" lst = new List"lt"double"gt"();
        for (int j = 0; j "lt" usecols.Length; ++j)
          lst.Add(double.Parse(tokens[usecols[j]]));
        double[] row = lst.ToArray();
        result.Add(row);
      }
      sr.Close(); ifs.Close();
      return result.ToArray();
    }

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

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

      double[][] result = MatMake(aRows, 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] += A[i][k] * B[k][j];

      return result;
    }

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

  } // class Program

} // ns