Bottom line: Right pseudo-inverse is virtually never used in machine learning scenarios, but I put together a demo anyway.

One morning before work, I was waiting for the rain to stop so I could give my two dogs their early walk. I figured I’d put together a demo of computing the right pseudo-inverse of a matrix using the C# language. I have implemented left-pseudo inverse many times — it’s a useful technique for training linear models — but I haven’t used right-pseudo-inverse in recent memory.

Briefly, left-pseudo-inverse is helpful in machine learning scenarios where a matrix of training data has more rows than columns. The right-pseudo inverse is helpful in scientific scenarios where a matrix of some sort of data has more columns than rows.

If you have a square matrix A, the inverse of A, Ainv, is a matrix so that A * Ainv = Ainv * A = I, where * is matrix multiplication and I is the identity matrix (1.0s on the upper-left to lower-right diagonal, 0.0s elsewhere). Not all square matrices have an inverse.

For non-square matrices, it’s useful to compute a pseudo-inverse. In machine learning scenarios, the matrix to invert is usually training data, or a design matrix (training data with a leading column of 1.0s added). In these scenarios, the number of rows (number of training items) is almost always greater than the number of columns (predictor variables).

There are two main categories of pseudo-inverse algorithms used in machine learning training scenarios: 1.) relaxed Moore-Penrose pseudo-inverse (using either QR or SVD inverse, each with several variations), 2.) left pseudo-inverse (almost always using Cholesky inverse). Relaxed Moore-Penrose pseudo-inverse is extremely complex. Left pseudo-inverse is simpler but works only when nRows is greater than nCols — which is the case for machine learning scenarios.

If you have a matrix A with more columns than rows (sometimes occurs in scientific programming but rarely in machine learning scenarios), the right pseudo-inverse via normal equations is computed as right-pinv(A) = At * inv(A * At). Here inv() is any regular matrix inverse function, but because A * At is square symmetric positive definite (assuming you condition the diagonal by adding a small constant, typically about .0e-8), you can apply the relatively simple Cholesky inverse function.

In code:

public class Cholesky
{
  public static double[][] MatRightPseudoInv(double[][] A)
  {
    // right-pinv(A) = At * inv(A * At)
    double[][] At = MatTranspose(A);
    double[][] AAt = MatProduct(A, At);
    for (int i = 0; i < A.Length; ++i) // condition AAt
      A[i][i] += 1.0e-8;
    double[][] AAtinv = MatInvCholesky(AAt);
    double[][] result = MatProduct(At, AAtinv);
    return result;
  }
  . . . 

I wrote a demo program that generates random matrices, computes the right pseudo-inverse, and checks the result. The right pseudo-inverse worked reasonably well (precision to within 1.0e-6 but no better than that).

The main problem with both left-pseudo inverse and right pseudo-inverse is the computation of At * A or A * At. If there are many small values, or many large values in At, then A * At can generate arithmetic overflow or underflow.

The 1930s, 40s, and 50s were the golden age of science fiction magazines. Big bugs made frequent cover appearances.

Left: Amazing Stories, May 1941.
Center: Thrilling Wonder Stories, December 1939.
Right: Fantastic Adventures, February 1940.

Demo program. Replace "lt" (less than), "gt", "lte", "gte" with Boolean operator symbols -- my blog editor chokes on symbols.

using System.IO;

namespace MatrixRightPseudoInverse
{
  internal class MatrixRightPseudoInverseProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin right pseudo-inverse" +
        " via normal equations (Cholesky) ");
      Console.WriteLine("Used when nCols gt nRows (rare) ");

      int minRows = 2; int maxRows = 10;
      int minCols = 10; int maxCols = 1000;
      Random rnd = new Random(0);
      int nTrials = 1000;
      int nPass = 0; int nFail = 0;
      for (int i = 0; i "lt" nTrials; ++i)
      {
        Console.WriteLine("\n==========");
        Console.WriteLine("trial " + i);
        double[][] A =
          MakeRandomMatrix(minRows, maxRows,
          minCols, maxCols, rnd);

        if (A.Length "gt" A[0].Length)
          continue;  // right pseudo-inv doesn't apply
        double[][] Apinv = Cholesky.MatRightPseudoInv(A);

        // check (A * Apinv) * A = A
        double[][] AApinv = MatProduct(A, Apinv);
        double[][] C = MatProduct(AApinv, A);

        //Console.WriteLine("\nA = ");
        //MatShow(A, 4, 9);
        //Console.WriteLine("\nApinv = ");
        //MatShow(Apinv, 4, 9);
        //Console.WriteLine("\n(A * Apinv) * A = ");
        //MatShow(C, 4, 9);

        if (MatAreEqual(C, A, 1.0e-6) == true)
        {
          Console.WriteLine("pass");
          ++nPass;
        }
        else
        {
          Console.WriteLine("FAIL");
          Console.WriteLine("nRows = " + A.Length +
            " nCols = " + A[0].Length);
          Console.ReadLine();
          ++nFail;
        }
        Console.WriteLine("==========");
     } // nTrials

      Console.WriteLine("\nNumber pass = " + nPass);
      Console.WriteLine("Number fail = " + nFail);

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

    // ------------------------------------------------------
    // helpers: MakeRandomMatrix(), MatProduct(), 
    // MatShow(), MatAreEqual()
    // ------------------------------------------------------

    public static double[][] MakeRandomMatrix(int minRows,
      int maxRows, int minCols, int maxCols, Random rnd)
    {
      int nRows = rnd.Next(minRows, maxRows);  // [2,maxN)
      int nCols = rnd.Next(minCols, maxCols);
      double[][] result = new double[nRows][];
      for (int i = 0; i "lt" nRows; ++i)
        result[i] = new double[nCols];

      double lo = -10.0; double hi = 10.0;
      for (int r = 0; r "lt" nRows; ++r)
        for (int c = 0; c "lt" nCols; ++c)
          result[r][c] = (hi - lo) * rnd.NextDouble() + lo;
      return result;
    }

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

    public static double[][] MatProduct(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 = 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] += A[i][k] * B[k][j];

      return result;
    }

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

    public 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("");
      }
    }

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

    public static bool MatAreEqual(double[][] A, double[][] B,
      double epsilon)
    {
      int nRows = A.Length; int nCols = A[0].Length;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          if (Math.Abs(A[i][j] - B[i][j]) "gt" epsilon)
            return false;
      return true;
    }

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

  } // Program

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

  public class Cholesky
  {
    public static double[][] MatRightPseudoInv(double[][] A)
    {
      // right-pinv(A) = At * inv(A * At)
      double[][] At = MatTranspose(A);
      double[][] AAt = MatProduct(A, At);
      // condition AAt
      for (int i = 0; i "lt" A.Length; ++i)
        A[i][i] += 1.0e-8;
      double[][] AAtinv = MatInvCholesky(AAt);
      double[][] result = MatProduct(At, AAtinv);
      return result;
    }
    public static double[][] MatLeftPseudoInv(double[][] A)
    {
      // left-pinv(A) = inv(At * A) * At
      double[][] At = MatTranspose(A);
      double[][] AtA = MatProduct(At, A);
      // condition AtA
      for (int i = 0; i "lt" A.Length; ++i)
        A[i][i] += 1.0e-8;
      double[][] AtAinv = MatInvCholesky(AtA);
      double[][] result = MatProduct(AtAinv, At);
      return result;
    }

    private static double[][] MatInvCholesky(double[][] A)
    {
      // A must be square symmetric positive definite
      // A = L * Lt, inv(A) = inv(Lt) * inv(L)
      double[][] L = MatDecompCholesky(A); // lower tri
      double[][] Lt = MatTranspose(L); // upper tri
      double[][] invL = MatInvLowerTri(L);
      double[][] invLt = MatInvUpperTri(Lt);
      double[][] result = MatProduct(invLt, invL);
      return result;
    }

    private static double[][] MatDecompCholesky(double[][] A)
    {
      // A must be square, symmetric, positive definite
      int m = A.Length; int n = A[0].Length;  // m == n
      double[][] L = new double[n][];
      for (int i = 0; i "lt" n; ++i)
        L[i] = new double[n];

      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 += L[i][k] * L[j][k];
          if (i == j)
          {
            double tmp = A[i][i] - sum;
            if (tmp "lt" 0.0)
              throw new
                Exception("decomp Cholesky fatal");
            L[i][j] = Math.Sqrt(tmp);
          }
          else
          {
            if (L[j][j] == 0.0)
              throw new
                Exception("decomp Cholesky fatal ");
            L[i][j] = (A[i][j] - sum) / L[j][j];
          }
        } // j
      } // i

      return L;
    }

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

    private static double[][] MatInvLowerTri(double[][] A)
    {
      int n = A.Length;  // must be square matrix
      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] * A[k][i];
          }
          result[k][j] /= (A[k][k] + 1.0e-8);
        }
      }
      return result;
    }

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

    private static double[][] MatInvUpperTri(double[][] A)
    {
      int n = A.Length;  // must be square matrix
      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[j][k] -= result[j][i] * A[i][k];
          }
          result[j][k] /= (A[k][k] + 1.0e-8); // avoid 0
        }
      }
      return result;
    }

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

    private 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;
    }

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

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

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

    private static double[][] MatTranspose(double[][] M)
    {
      int nRows = M.Length; int nCols = M[0].Length;
      double[][] result = MatMake(nCols, nRows);  // note
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[j][i] = M[i][j];
      return result;
    }

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

    private static double[][] MatProduct(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 Cholesky

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

} // ns