Inverting a matrix is one of the most difficult problems in numerical programming. A special case of inverting a matrix M is when M is square, symmetric, and positive definite. In machine learning, one scenario where this happens is when M is a covariance matrix, which is an intermediate result during Gaussian process regression.

To invert a special M, you can use one of several standard algorithms, such as inversion with Crout’s decomposition. But for a square, symmetric, positive definite (PD) matrix, you can use Cholesky decomposition, which is a bit easier and slightly more efficient.

I implemented a C# function MatInverseFromCholesky(double[][] L) that accepts the result of Cholesky decomposition and returns the inverse of the source matrix. The calling code looks like:

double[][] source = // set up PD matrix values
double[][] lower = MatCholesky(source);  // decompose
double[][] inverse = MatInverseFromCholesky(lower); 

It would be possible to combine the MatCholesky() and MatInverseFromCholesky() functions into a single MatInverseViaCholesky() function.

The MatInverseFromCholesky() function implementation below uses a direct approach, which is most efficient. An alternative approach uses helper functions. The idea is based on the fact that computing the inverse of a lower triangular matrix is relatively easy, and computing the inverse of an upper triangular matrix is relatively easy, and the inverse(A * B) = inverse(B) * inverse(A):

A = L * trans(L) 
inv(A) = inv(L * trans(L))
inv(A) = inv(trans(L) * inv(L)
inv(A) = inverse of an upper * inverse of a lower 

In code:

// L is Cholesky decomp of a special PD matrix A
double[][] Lt = MatTranspose(L);
double[][] Lti = MatInverseUpperTri(Lt);
double[][] Li = MatInverseLowerTri(L);
double[][] result = MatProduct(Lti, Li);
return result;

A fun exercise.

“The Matrix” (1999) is considered by many people, including me, to be one of the best sci-fi movies of all time. According to one of my friends who is an artist, a fairly common idea for graphic artists is creating alternative advertising posters. Left: The original poster from 1999 emphasized the actors, all of whom were reasonably well-known but not big stars at the time. Center: Here’s an alternative design that has a retro minimalist feel. Nice. Right: Here’s an alternative design that has a sort of mysterious feel. Nice again.

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

using System;
using System.IO;

namespace MatrixLibraryCholesky
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# matrix library demo ");
      Console.WriteLine("Inverse Cholesky decomposition ");

      double[][] mat = MatCreate(4, 4);
      mat[0] = new double[] { 1.0, 0.2, 0.3, 0.4 };
      mat[1] = new double[] { 0.2, 1.0, 0.5, 0.6 };
      mat[2] = new double[] { 0.3, 0.5, 1.0, 0.7 };
      mat[3] = new double[] { 0.4, 0.6, 0.7, 1.0 };

      double[][] lower = MatCholesky(mat); // decompose

      Console.WriteLine("\nSource matrix: ");
      MatShow(mat, 1, 5);
      Console.WriteLine("\nCholesky lower: ");
      MatShow(lower, 4, 8);

      double[][] inverse = MatInverseFromCholesky(lower);
      Console.WriteLine("\nCholesky inverse: ");
      MatShow(inverse, 4, 8);

      double[][] product1 = MatProduct(mat, inverse);
      Console.WriteLine("\nsource * inverse: ");
      MatShow(product1, 4, 8);

      double[][] product2 = MatProduct(inverse, mat);
      Console.WriteLine("\ninverse * source: ");
      MatShow(product2, 4, 8);

      Console.WriteLine("\nEnd ");
      Console.ReadLine();

    }

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

      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) // could use bRows
              result[i][j] += matA[i][k] * matB[k][j];

        return result;
      }

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

      static double[][] MatCholesky(double[][] m)
      {
        // Cholesky decomposition
        // m is square, symmetric, positive definite
        int n = m.Length;
        double[][] result = MatCreate(n, n);  // all 0.0
        for (int i = 0; i "lt" n; ++i)
        {
          for (int j = 0; j "lte" i; ++j)
          {
            double sum = 0.0;
            for (int k = 0; k "lt" j; ++k)
              sum += result[i][k] * result[j][k];
            if (i == j)
            {
              double tmp = m[i][i] - sum;
              if (tmp "lt" 0.0)
                throw new Exception("MatCholesky fatal error ");
              result[i][j] = Math.Sqrt(tmp);
            }
            else
            {
              if (result[j][j] == 0.0)
                throw new Exception("MatCholesky fatal error ");
              result[i][j] = (1.0 / result[j][j] * (m[i][j] - sum));
            }
          } // j
        } // i
        return result;
      }

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

      public static double[][]
        MatInverseFromCholesky(double[][] L)
      {
        // inverse of a special PD source matrix 
        // that has been decomposed to L
        int n = L.Length;
        double[][] result = MatIdentity(n);

        for (int k = 0; k "lt" n; ++k)
        {
          for (int j = 0; j "lt" n; j++)
          {
            for (int i = 0; i "lt" k; i++)
            {
              result[k][j] -= result[i][j] * L[k][i];
            }
            result[k][j] /= L[k][k];
          }
        }

        for (int k = n - 1; k "gte" 0; --k)
        {
          for (int j = 0; j "lt" n; j++)
          {
            for (int i = k + 1; i "lt" n; i++)
            {
              result[k][j] -= result[i][j] * L[i][k];
            }
            result[k][j] /= L[k][k];
          }
        }
        return result;
      }

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

      static double[][] MatIdentity(int n)
      {
        double[][] result = MatCreate(n, n);
        for (int i = 0; i "lt" n; ++i)
          result[i][i] = 1.0;
        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;  // avoid "-0.00"
            Console.Write(v.ToString("F" + dec).PadLeft(wid));
          }
          Console.WriteLine("");
        }
      }

    } // Program

  } // ns