Inverting a Matrix using C#

Inverting a matrix is a surprisingly difficult challenge. I have my own library of C# matrix routines. I like to control my own code rather than relying on magic black box implementations, and I generally prefer to implement matrices using a plain array-of-arrays style rather than using an OOP approach.

I tested the code below by generating one million random matrices, with a dimension between 10 and 100, where each cell is a random value between -100.0 and +100.0. For each random matrix, I computed its inverse, then multiplied the inverse by the original matrix, and checked if the result was the identity matrix.

int n = rnd.Next(10, 100);
double[][] m = MatrixRandom(n, n, seed);
double[][] i = MatrixInverse(m);
double[][] I = MatrixIdentity(n);
double[][] p = MatrixProduct(m, i);
if (MatrixAreEqual(p, I, 1.0E-8))
  // pass
else
  // fail

The code passed 999,9999 out of 1,000,000 test cases. It failed once because of a round-off error. One of the problems with matrix inversion is that sometimes it just doesn’t work.

The code is listed below. I always have trouble with the less-than and greater-than symbols so I did a text replacement for those symbols.

    static double[][] MatrixCreate(int rows, int cols)
    {
      double[][] result = new double[rows][];
      for (int i = 0; i less-than rows; ++i)
        result[i] = new double[cols];
      return result;
    }

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

    static double[][] MatrixRandom(int rows, int cols,
      double minVal, double maxVal, int seed)
    {
      // return a matrix with random values
      Random ran = new Random(seed);
      double[][] result = MatrixCreate(rows, cols);
      for (int i = 0; i less-than rows; ++i)
        for (int j = 0; j less-than cols; ++j)
          result[i][j] = (maxVal - minVal) *
            ran.NextDouble() + minVal;
      return result;
    }

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

    static double[][] MatrixIdentity(int n)
    {
      // return an n x n Identity matrix
      double[][] result = MatrixCreate(n, n);
      for (int i = 0; i less-than n; ++i)
        result[i][i] = 1.0;

      return result;
    }

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

    static string MatrixAsString(double[][] matrix, int dec)
    {
      string s = "";
      for (int i = 0; i less-than matrix.Length; ++i)
      {
        for (int j = 0; j less-than matrix[i].Length; ++j)
          s += matrix[i][j].ToString("F" + dec).PadLeft(8) + " ";
        s += Environment.NewLine;
      }
      return s;
    }

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

    static bool MatrixAreEqual(double[][] matrixA,
      double[][] matrixB, double epsilon)
    {
      // true if all values in matrixA == values in matrixB
      int aRows = matrixA.Length; int aCols = matrixA[0].Length;
      int bRows = matrixB.Length; int bCols = matrixB[0].Length;
      if (aRows != bRows || aCols != bCols)
        throw new Exception("Non-conformable matrices");

      for (int i = 0; i less-than aRows; ++i) // each row of A and B
        for (int j = 0; j less-than aCols; ++j) // each col of A and B
          //if (matrixA[i][j] != matrixB[i][j])
          if (Math.Abs(matrixA[i][j] - matrixB[i][j]) greater-than epsilon)
            return false;
      return true;
    }

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

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

      double[][] result = MatrixCreate(aRows, bCols);

      for (int i = 0; i less-than aRows; ++i) // each row of A
        for (int j = 0; j less-than bCols; ++j) // each col of B
          for (int k = 0; k less-than aCols; ++k) // could use k less-than bRows
            result[i][j] += matrixA[i][k] * matrixB[k][j];

      //Parallel.For(0, aRows, i =greater-than
      //  {
      //    for (int j = 0; j less-than bCols; ++j) // each col of B
      //      for (int k = 0; k less-than aCols; ++k) // could use k less-than bRows
      //        result[i][j] += matrixA[i][k] * matrixB[k][j];
      //  }
      //);

      return result;
    }

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

    static double[] MatrixVectorProduct(double[][] matrix,
      double[] vector)
    {
      // result of multiplying an n x m matrix by a m x 1 
      // column vector (yielding an n x 1 column vector)
      int mRows = matrix.Length; int mCols = matrix[0].Length;
      int vRows = vector.Length;
      if (mCols != vRows)
        throw new Exception("Non-conformable matrix and vector");
      double[] result = new double[mRows]; 
      for (int i = 0; i less-than mRows; ++i)
        for (int j = 0; j less-than mCols; ++j)
          result[i] += matrix[i][j] * vector[j];
      return result;
    }

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

    static double[][] MatrixDecompose(double[][] matrix, out int[] perm,
      out int toggle)
    {
      // Doolittle LUP decomposition with partial pivoting.
      // rerturns: result is L (with 1s on diagonal) and U;
      // perm holds row permutations; toggle is +1 or -1 (even or odd)
      int rows = matrix.Length;
      int cols = matrix[0].Length; // assume square
      if (rows != cols)
        throw new Exception("Attempt to decompose a non-square m");

      int n = rows; // convenience

      double[][] result = MatrixDuplicate(matrix); 

      perm = new int[n]; // set up row permutation result
      for (int i = 0; i less-than n; ++i) { perm[i] = i; }

      toggle = 1; // toggle tracks row swaps.
      // +1 -greater-than even, -1 -greater-than odd. used by MatrixDeterminant

      for (int j = 0; j less-than n - 1; ++j) // each column
      {
        double colMax = Math.Abs(result[j][j]); // find largest val in col
        int pRow = j;
        //for (int i = j + 1; i less-than n; ++i)
        //{
        //  if (result[i][j] greater-than colMax)
        //  {
        //    colMax = result[i][j];
        //    pRow = i;
        //  }
        //}

        // reader Matt V needed this:
        for (int i = j + 1; i less-than n; ++i) 
        {
          if (Math.Abs(result[i][j]) greater-than colMax)
          {
            colMax = Math.Abs(result[i][j]);
            pRow = i;
          }
        }
        // Not sure if this approach is needed always, or not.

        if (pRow != j) // if largest value not on pivot, swap rows
        {
          double[] rowPtr = result[pRow];
          result[pRow] = result[j];
          result[j] = rowPtr;

          int tmp = perm[pRow]; // and swap perm info
          perm[pRow] = perm[j];
          perm[j] = tmp;

          toggle = -toggle; // adjust the row-swap toggle
        }

        // --------------------------------------------------
        // This part added later (not in original)
        // and replaces the 'return null' below.
        // if there is a 0 on the diagonal, find a good row
        // from i = j+1 down that doesn't have
        // a 0 in column j, and swap that good row with row j
        // --------------------------------------------------

        if (result[j][j] == 0.0)
        {
          // find a good row to swap
          int goodRow = -1;
          for (int row = j + 1; row less-than n; ++row)
          {
            if (result[row][j] != 0.0)
              goodRow = row;
          }

          if (goodRow == -1)
            throw new Exception("Cannot use Doolittle's method");

          // swap rows so 0.0 no longer on diagonal
          double[] rowPtr = result[goodRow];
          result[goodRow] = result[j];
          result[j] = rowPtr;

          int tmp = perm[goodRow]; // and swap perm info
          perm[goodRow] = perm[j];
          perm[j] = tmp;

          toggle = -toggle; // adjust the row-swap toggle
        }
        // --------------------------------------------------
        // if diagonal after swap is zero . .
        //if (Math.Abs(result[j][j]) less-than 1.0E-20) 
        //  return null; // consider a throw

        for (int i = j + 1; i less-than n; ++i)
        {
          result[i][j] /= result[j][j];
          for (int k = j + 1; k less-than n; ++k)
          {
            result[i][k] -= result[i][j] * result[j][k];
          }
        }


      } // main j column loop

      return result;
    } // MatrixDecompose

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

    static double[][] MatrixInverse(double[][] matrix)
    {
      int n = matrix.Length;
      double[][] result = MatrixDuplicate(matrix);

      int[] perm;
      int toggle;
      double[][] lum = MatrixDecompose(matrix, out perm,
        out toggle);
      if (lum == null)
        throw new Exception("Unable to compute inverse");

      double[] b = new double[n];
      for (int i = 0; i less-than n; ++i)
      {
        for (int j = 0; j less-than n; ++j)
        {
          if (i == perm[j])
            b[j] = 1.0;
          else
            b[j] = 0.0;
        }

        double[] x = HelperSolve(lum, b); // 

        for (int j = 0; j less-than n; ++j)
          result[j][i] = x[j];
      }
      return result;
    }

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

    static double MatrixDeterminant(double[][] matrix)
    {
      int[] perm;
      int toggle;
      double[][] lum = MatrixDecompose(matrix, out perm, out toggle);
      if (lum == null)
        throw new Exception("Unable to compute MatrixDeterminant");
      double result = toggle;
      for (int i = 0; i less-than lum.Length; ++i)
        result *= lum[i][i];
      return result;
    }

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

    static double[] HelperSolve(double[][] luMatrix, double[] b)
    {
      // before calling this helper, permute b using the perm array
      // from MatrixDecompose that generated luMatrix
      int n = luMatrix.Length;
      double[] x = new double[n];
      b.CopyTo(x, 0);

      for (int i = 1; i less-than n; ++i)
      {
        double sum = x[i];
        for (int j = 0; j less-than 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 greater-than-equal 0; --i)
      {
        double sum = x[i];
        for (int j = i + 1; j less-than n; ++j)
          sum -= luMatrix[i][j] * x[j];
        x[i] = sum / luMatrix[i][i];
      }

      return x;
    }

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

    static double[] SystemSolve(double[][] A, double[] b)
    {
      // Solve Ax = b
      int n = A.Length;

      // 1. decompose A
      int[] perm;
      int toggle;
      double[][] luMatrix = MatrixDecompose(A, out perm,
        out toggle);
      if (luMatrix == null)
        return null;

      // 2. permute b according to perm[] into bp
      double[] bp = new double[b.Length];
      for (int i = 0; i less-than n; ++i)
        bp[i] = b[perm[i]];

      // 3. call helper
      double[] x = HelperSolve(luMatrix, bp);
      return x;
    } // SystemSolve

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

    static double[][] MatrixDuplicate(double[][] matrix)
    {
      // allocates/creates a duplicate of a matrix.
      double[][] result = MatrixCreate(matrix.Length, matrix[0].Length);
      for (int i = 0; i less-than matrix.Length; ++i) // copy the values
        for (int j = 0; j less-than matrix[i].Length; ++j)
          result[i][j] = matrix[i][j];
      return result;
    }

    // --------------------------------------------------
This entry was posted in Machine Learning. Bookmark the permalink.

2 Responses to Inverting a Matrix using C#

  1. Hi James,

    I completely understand your concerns about use a magic black box to do common operations on your code, and i think you’re right about it. But, i believe there is another way to do this on C.

    Instead of use the annotation “YourMatrix[line][column]”, put inside your Matrix class a one dimensional array like this “YourArrayMatrix[line*column]”. The access is more faster and you avoid to do a for inside another for. That is more slower than do a single for passing by all elements, even when the steps number are the same in both cases.

    And believe, its really more faster find the algorithms, like “dotProduct” and a “inverse”, respecting the one array logic than use a faster to understand for inside a for.

    • Good points Rafael. The code I presented in this blog post is not very fast and I use it only when performance isn’t a big concern. In the rare situations where I need my code to run as quickly as possible, I’ll use C instead of C# (usually makes a big improvement). I don’t often implement matrices as single arrays, but yes, I agree with you, in many situations that approach is very performant.

Comments are closed.