I have a personal library of C# matrix and vector functions. One of the key functions is a matrix inverse function — one of the most complex problems in numerical programming.

One of my matrix inverse functions, which I’ve used for a long time, is made of four related functions: MatInverse(), MatDecompose(), MatCreate(), Reduce(). The top-level MatInverse() function calls MatCreate(), MatDecompose(), Reduce(). This version uses LUP (“lower, upper, permutation”) decomposition via Crout’s algorithm. Other ways to compute an inverse include QR decomposition (but you can’t get the determinant), and SVD (“singular value”) which is more complex than LUP.

One morning before work, I decided to refactor the MatInverse() function so that it is completely self-contained, by placing the three helpers as nested functions inside the top-level function. I named the self-contained function MatInverseLU() because the MatDecompose() function uses Crout’s algorithm, as opposed to Doolittle’s or some other algorithm.

The refactored version of matrix inverse worked as expected.

A downside to the self-contained design is that a MatDeterminant() function also calls MatDecompose() and so you’d have to replicate that LU decomposition code inside the MatDeterminant() function.

One of the problems with music albums in the 1960s was that a typical album had a few good songs and a few throw-away songs. Here are three albums that are self-contained in the sense that every song is good (in my opinion of course — music is incredibly subjective and depends on your memory of what you were doing when you heard the song). Clips: “You Baby”, “A Day in the Life”, “Sit Down I Think I Love You”.

Demo code. Replace “lt”, “gt”, “lte”, “gte” with Boolean operator symbols.

using System;
using System.IO;

namespace MatrixLibrary
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# matrix library demo ");
      Console.WriteLine("Matrix inverse self-contained LU ");

      double[][] A = MatCreate(4, 4);
      A[0] = new double[] { 1, 2, 3, 4 };
      A[1] = new double[] { 6, 5, 6, 7 };
      A[2] = new double[] { 8, 9, 8, 9 };
      A[3] = new double[] { 2, 4, 3, 7 };

      Console.WriteLine("\nSource matrix A = ");
      MatShow(A, 2, 8);

      double[][] Ai = MatInverseLU(A);
      Console.WriteLine("\nInverse from LU decomp: ");
      MatShow(Ai, 4, 12);

      double[][] AiA = MatProduct(Ai, A);
      Console.WriteLine("\nVerifying Ai * A = I ");
      MatShow(AiA, 4, 12);

      double[][] AAi = MatProduct(A, Ai);
      Console.WriteLine("\nVerifying A * Ai = I ");
      MatShow(AAi, 4, 12);

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

    static double[][] MatInverseLU(double[][] m)
    {
      // inverse from Crout's decomp
      // helpers fully contained
      int n = m.Length;
      double[][] result = MatCreateHelp(n, n); // make copy
      for (int i = 0; i "lt" n; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = m[i][j];

      double[][] lum; // combined lower and upper
      int[] perm;  // out parameter
      MatDecomposeHelp(m, out lum, out perm);

      double[] b = new double[n];
      for (int i = 0; i "lt" n; ++i)
      {
        for (int j = 0; j "lt" n; ++j)
          if (i == perm[j])
            b[j] = 1.0;
          else
            b[j] = 0.0;

        double[] x = ReduceHelp(lum, b); // 
        for (int j = 0; j "lt" n; ++j)
          result[j][i] = x[j];
      }
      return result;

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

      static double[][] MatCreateHelp(int rows, int cols)
      {
        double[][] result = new double[rows][];
        for (int i = 0; i "lt" rows; ++i)
          result[i] = new double[cols];
        return result;
      }

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

      static int MatDecomposeHelp(double[][] m,
        out double[][] lum, out int[] perm)
      {
        // Crout's LU decomposition for matrix 
        // determinant and inverse
        // stores combined lower & upper in lum[][]
        // stores row permuations into perm[]
        // returns +1 or -1 according to even or
        // odd number of row permutations
        // lower gets dummy 1.0s on diagonal (0.0s above)
        // upper gets lum values on diagonal (0.0s below)

        int toggle = +1;
        int n = m.Length;

        // make a copy of m[][] into result lu[][]
        lum = MatCreateHelp(n, n);
        for (int i = 0; i "lt" n; ++i)
          for (int j = 0; j "lt" n; ++j)
            lum[i][j] = m[i][j];

        // make perm[]
        perm = new int[n];
        for (int i = 0; i "lt" n; ++i)
          perm[i] = i;

        for (int j = 0; j "lt" n - 1; ++j) // note n-1 
        {
          double max = Math.Abs(lum[j][j]);
          int piv = j;

          for (int i = j + 1; i "lt" n; ++i) // pivot index
          {
            double xij = Math.Abs(lum[i][j]);
            if (xij "gt" max)
            {
              max = xij;
              piv = i;
            }
          } // i

          if (piv != j)
          {
            double[] tmp = lum[piv]; // swap rows j, piv
            lum[piv] = lum[j];
            lum[j] = tmp;

            int t = perm[piv]; // swap perm elements
            perm[piv] = perm[j];
            perm[j] = t;

            toggle = -toggle;
          }

          double xjj = lum[j][j];
          if (xjj != 0.0)
          {
            for (int i = j + 1; i "lt" n; ++i)
            {
              double xij = lum[i][j] / xjj;
              lum[i][j] = xij;
              for (int k = j + 1; k "lt" n; ++k)
                lum[i][k] -= xij * lum[j][k];
            }
          }

        } // j

        return toggle;  // for determinant
      } // MatDecomposeHelp

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

      static double[] ReduceHelp(double[][] luMatrix, double[] b) 
      {
        int n = luMatrix.Length;
        double[] x = new double[n];
        //b.CopyTo(x, 0);
        for (int i = 0; i "lt" n; ++i)
          x[i] = b[i];

        for (int i = 1; i "lt" n; ++i)
        {
          double sum = x[i];
          for (int j = 0; j "lt" i; ++j)
            sum -= luMatrix[i][j] * x[j];
          x[i] = sum;
        }

        x[n - 1] /= luMatrix[n - 1][n - 1];
        for (int i = n - 2; i "gt"= 0; --i)
        {
          double sum = x[i];
          for (int j = i + 1; j "lt" n; ++j)
            sum -= luMatrix[i][j] * x[j];
          x[i] = sum / luMatrix[i][i];
        }

        return x;
      } // ReduceHelp

    } // MatInverseLU

    // ======================================================

    static double[][] MatCreate(int rows, int cols)
    {
      double[][] result = new double[rows][];
      for (int i = 0; i "lt" rows; ++i)
        result[i] = new double[cols];
      return result;
    }

    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 = MatCreate(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] += matA[i][k] * matB[k][j];

      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];
          if (Math.Abs(v) "lt" 1.0e-5) v = 0.0; // hack
          Console.Write(v.ToString("F" + dec).PadLeft(wid));
        }
        Console.WriteLine("");
      }
    }

  } // Program
} // ns