In machine learning, using the relaxed Moore-Penrose pseudo-inverse function is one way to train a linear model (typically linear regression, or quadratic regression). The MP pseudo-inverse is extrememly complicated — the standard LAPACK implementation used by most libraries is over 10,000 lines of Fortran code.

There are several algorithms that can be used to compute the MP pseudo-inverse: QR Householder decomposition, QR Gram-Schmidt decomposition, QR Givens decomposition, SVD Golub-Reinsch, SVD Jacobi methods, SVD Lanczos, SVD Householder+QR. And as if this isn’t enough, each of these sub-variations has several significantly different implementation designs.

One morning before work, I figured I’d run my current version of relaxed Moore-Penrose pseudo-inverse via QR decomposition using the modified Gram-Schmidt algorithm through some tests. Bottom line: the function passed 1,000 tests, without any problems.

The fully compliant MP pseudo-inverse function must satisify several mathematical conditions. But for machine learning scenarios, it’s enough to satisfy (A * Apinv) * A = A, where A is a matrix (not necessarily square) and Apinv is the pseudo-inverse of A.

My pseudo-inverse function is limited to the case where the source matrix has more rows than columns — almost always the case for machine learning training data.

To test I used this approach:

loop 1000 times
  create random matrix A with between 200 and 10,000 rows
   and between 2 and 20 columns, where each value is
   between -10.0 and +10.0

  compute pinv(A)

  check if (A * Apinv) * A equals A (within 1.0e-8)
end-loop

I’m sure I could come up with some examples where the pseudo-inverse function would fail, for example by using more columns or a bigger range of values in the matrices, but I’m satisfied the relaxed MP pseudo-inverse implementation works well enough for most machine learning scenarios.

The Adventures of Blake and Mortimer is a Belgian comics series created by writer and artist Edgar P. Jacobs (1904-1987). The series is similar to the Tintin series in art style but the stories are more serious. There are currently 31 books in the series. The first 11 were published from 1950 to 1977, then after Jacobs retired, 20 more were created by various authors and artists.

I like the Blake and Mortimer series, but I prefer the Tintin series.

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

// test scratch pseudo-inv via QR Gram-Schmidt
// verify (A * Apinv) * A = A

using System.IO;

namespace MatrixPseudoInverseQR_GramSchmidt_Test
{
  internal class MatrixPseudoInverseQRGSTestProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin test pinv using " +
        " QR Gram-Schmidt ");

      int minRows = 100; int maxRows = 10000;
      int minCols = 2; int maxCols = 20;
      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);

        double[][] Apinv = QRGS.MatPseudoInv(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-8) == true)
        {
          Console.WriteLine("pass");
          ++nPass;
        }
        else
        {
          Console.WriteLine("FAIL");
          Console.WriteLine("nRows = " + A.Length +
            " nCols = " + A[0].Length);
          //Console.ReadLine();
          ++nFail;
        }
        Console.WriteLine("==========");
        // Console.ReadLine();
      } // nTrials

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

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

    // ------------------------------------------------------
    // helpers: 
    // ------------------------------------------------------

    public static double[][] MakeRandomMatrix(int minRows,
      int maxRows, int minCols, int maxCols, Random rnd)
    {
      int nRows = rnd.Next(minRows, maxRows);
      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 QRGS
  {
    public static double[][] MatPseudoInv(double[][] M)
    {
      // relaxed Moore-Penrose pseudo-inverse using QR-GS
      // A = Q*R, pinv(A) = inv(R) * trans(Q)

      int nr = M.Length; int nc = M[0].Length; // aka m, n
      if (nr "lt" nc)
        Console.WriteLine("ERROR: Works only m "gte" n");

      double[][] Q; double[][] R;
      MatDecompQR(M, out Q, out R); // Gram-Schmidt

      double[][] Rinv = MatInvUpperTri(R); // std algo
      double[][] Qinv = MatTranspose(Q);  // is inv(Q)
      double[][] result = MatProduct(Rinv, Qinv);
      return result;
    }

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

    private static void MatDecompQR(double[][] A,
      out double[][] Q, out double[][] R)
    {
      // QR decomposition, modified Gram-Schmidt
      // 'reduced' mode (all machine learning scenarios)
      int m = A.Length; int n = A[0].Length;
      if (m "lt" n)
        throw new Exception("m must be gte n ");

      double[][] QQ = new double[m][]; // working Q mxn
      for (int i = 0; i "lt" m; ++i)
        QQ[i] = new double[n];

      double[][] RR = new double[n][]; // working R nxn
      for (int i = 0; i "lt" n; ++i)
        RR[i] = new double[n];

      for (int k = 0; k "lt" n; ++k) // main loop each col
      {
        double[] v = new double[m];
        for (int i = 0; i "lt" m; ++i)  // col k
          v[i] = A[i][k];
        for (int j = 0; j "lt" k; ++j)  // inner loop
        {
          double[] colj = new double[QQ.Length];
          for (int i = 0; i "lt" colj.Length; ++i)
            colj[i] = QQ[i][j];

          double vecdot = 0.0;
          for (int i = 0; i "lt" colj.Length; ++i)
            vecdot += colj[i] * v[i];
          RR[j][k] = vecdot;

          // v = v - (R[j, k] * Q[:, j])
          for (int i = 0; i "lt" v.Length; ++i)
            v[i] = v[i] - (RR[j][k] * QQ[i][j]);
        } // j

        double normv = 0.0;
        for (int i = 0; i "lt" v.Length; ++i)
          normv += v[i] * v[i];
        normv = Math.Sqrt(normv);
        RR[k][k] = normv;

        // Q[:, k] = v / R[k, k]
        for (int i = 0; i "lt" QQ.Length; ++i)
          QQ[i][k] = v[i] / (RR[k][k] + 1.0e-12);

      } // k

      Q = QQ;
      R = RR;
    } 

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

    private static double[][] MatInvUpperTri(double[][] A)
    {
      // used to invert R from QR
      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[][] 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[][] 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[][] 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 QRGS

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

} // ns