A matrix inverse applies only to square matrices. A pseudo-inverse applies to any shape matrix. In machine learning, the matrix pseudo-inverse has many uses, for example, to directly compute the weights and bias of a linear regression model (instead of an iterative SGD approach).

There are several different types of pseudo-inverse. The most common is the Moore-Penrose pseudo-inverse. There are several ways to compute a Moore-Penrose pseudo-inverse. The two most common that I come across are using SVD (singular value decomposition) and QR decomposition. There are several ways to compute both SVD and QR — for example, Jacobi and Golub-Kahan-Reisch for SVD, and Householder and Gram-Schmidt for QR. The point here is that there’s a bewildering number of ways to compute a pseudo-inverse.

One foggy Pacific Northwest morning, I decided to kill some time (I have a sleep disorder and can only sleep about one hour per day, at most, no matter how hard I try) by implementing Moore-Penrose pseudo-inverse using QR decomposition using the Householder algorithm using the C# language.

Here’s the output of my demo:

Begin Moore-Penrose pseudo-inverse using
  QR (Householder) with C#

Source matrix A:
   4.0000   7.0000   1.0000
   6.0000   0.0000   3.0000
   8.0000   1.0000   9.0000
   2.0000   5.0000   6.0000
   1.0000   5.0000   4.0000

Computing pseudo-inverse

The pseudo-inverse:
   0.0882   0.1016   0.0299  -0.0721  -0.0574
   0.0937  -0.0202  -0.0455   0.0323   0.0455
  -0.1041  -0.0478   0.0609   0.0825   0.0511

Checking A * Apinv * A
   4.0000   7.0000   1.0000
   6.0000  -0.0000   3.0000
   8.0000   1.0000   9.0000
   2.0000   5.0000   6.0000
   1.0000   5.0000   4.0000

End demo

I validated my demo using the Python NumPy library linald.pinv() function. The output:

Begin pseudo-inverse using NumPy

Source A:
[[4.0000 7.0000 1.0000]
 [6.0000 0.0000 3.0000]
 [8.0000 1.0000 9.0000]
 [2.0000 5.0000 6.0000]
 [1.0000 5.0000 4.0000]]

pinv(A):
[[ 0.0882  0.1016  0.0299 -0.0721 -0.0574]
 [ 0.0937 -0.0202 -0.0455  0.0323  0.0455]
 [-0.1041 -0.0478  0.0609  0.0825  0.0511]]

End demo

OK, that was fun. (If you’re reading this blog post you understand — the rest of the world doesn’t understand what fun is for guys like us.)

The “pseudo” in “pseudo-inverse” means “sort of”. I’ve read that the game of chess is “pseudo-warfare”. I don’t agree with that description — chess is just a game.

Right: I played Samuel Reshevshy (champion 1936, 1938, 1940, 1941, 1942, 1946, 1969) at the Whittier Chess Club. I played the black side of a French Defense, MacCutcheon variation and held him to a draw. Reshevsky was an unpleasant little man (about 5 feet 3 inches tall) and he was never invited back to Whittier.

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

using System;
using System.IO;
using System.Collections.Generic;

namespace MatrixPseudoInverseQR
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin Moore-Penrose" +
        " pseudo-inverse using QR (Householder)" +
        " with C# ");

      // to load matrix from file:
      // double[][] A = MatLoad("data.txt",
      //  new int[] { 0, 1, 2 }, ',', "#");

      double[][] A = new double[5][];  // 5x3
      A[0] = new double[] { 4.0, 7.0, 1.0 };
      A[1] = new double[] { 6.0, 0.0, 3.0 };
      A[2] = new double[] { 8.0, 1.0, 9.0 };
      A[3] = new double[] { 2.0, 5.0, 6.0 };
      A[4] = new double[] { 1.0, 5.0, 4.0 };

      Console.WriteLine("\nSource matrix A: ");
      MatShow(A, 4, 9);

      Console.WriteLine("\nComputing pseudo-inverse ");
      double[][] Apinv = MatPseudoInv(A);
      Console.WriteLine("\nThe pseudo-inverse: ");
      MatShow(Apinv, 4, 9);

      Console.WriteLine("\nChecking A * Apinv * A ");
      double[][] C = MatProduct(A, MatProduct(Apinv, A));
      MatShow(C, 4, 9);

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

    // ------------------------------------------------------
    // helpers for Main: MatShow, MatLoad, MatProduct
    // ------------------------------------------------------

    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];
          Console.Write(v.ToString("F" + dec).
            PadLeft(wid));
        }
        Console.WriteLine("");
      }
    }

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

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

      return result;
    }

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

    static double[][] MatLoad(string fn, int[] usecols,
      char sep, string comment)
    {
      List"lt"double[]"gt" result = 
        new List"lt"double[]"gt"();
      string line = "";
      FileStream ifs = new FileStream(fn, FileMode.Open);
      StreamReader sr = new StreamReader(ifs);
      while ((line = sr.ReadLine()) != null)
      {
        if (line.StartsWith(comment) == true)
          continue;
        string[] tokens = line.Split(sep);
        List"lt"double"gt" lst = new List"lt"double"gt"();
        for (int j = 0; j "lt" usecols.Length; ++j)
          lst.Add(double.Parse(tokens[usecols[j]]));
        double[] row = lst.ToArray();
        result.Add(row);
      }
      sr.Close(); ifs.Close();
      return result.ToArray();
    }

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

    static double[][] MatPseudoInv(double[][] M)
    {
      // Moore-Penrose pseudoinverse using QR decomposition
      // 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;
      MatDecomposeQR(M, out Q, out R, true); // reduced

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

      // ----------------------------------------------------
      // nested helpers: MatDecomposeQR (which has many
      // nested sub-helpers), MatInverseUpperTri,
      // MatTranspose, MatProduct
      // ----------------------------------------------------

      static void MatDecomposeQR(double[][] M,
         out double[][] Q, out double[][] R, bool reduced)
      {
        // QR decomposition, Householder algorithm.
        // see rosettacode.org/wiki/QR_decomposition
        int m = M.Length;
        int n = M[0].Length;

        if (m "lt" n)
          throw new Exception("Cannot do rows less than cols");

        double[][] QQ = MatIdentity(m); // working Q
        double[][] RR = MatCopy(M); // working R

        int end;
        if (m == n) end = n - 1;
        else end = n;

        for (int i = 0; i "lt" end; ++i)
        {
          double[][] H = MatIdentity(m);
          double[] a = new double[m - i]; // corr
          int k = 0;
          for (int ii = i; ii "lt" m; ++ii) // corr
            a[k++] = RR[ii][i];

          double normA = VecNorm(a);
          if (a[0] "lt" 0.0 && normA "gt" 0.0) // corr
            normA = -normA;
          else if (a[0] "gt" 0.0 && normA "lt" 0.0)
            normA = -normA;

          double[] v = new double[a.Length];
          for (int j = 0; j "lt" v.Length; ++j)
            v[j] = a[j] / (a[0] + normA);
          v[0] = 1.0;

          // Householder algorithm
          double[][] h = MatIdentity(a.Length);
          double vvDot = VecDot(v, v);
          double[][] A = VecToMat(v, v.Length, 1);
          double[][] B = VecToMat(v, 1, v.Length);
          double[][] AB = MatMult(A, B);

          for (int ii = 0; ii "lt" h.Length; ++ii)
            for (int jj = 0; jj "lt" h[0].Length; ++jj)
              h[ii][jj] -= (2.0 / vvDot) * AB[ii][jj];

          // copy h[][] into lower right corner of H[][]
          int d = m - h.Length; // corr
          for (int ii = 0; ii "lt" h.Length; ++ii)
            for (int jj = 0; jj "lt" h[0].Length; ++jj)
              H[ii + d][jj + d] = h[ii][jj];

          QQ = MatMult(QQ, H);
          RR = MatMult(H, RR);
        } // i

        if (reduced == false)
        {
          Q = QQ; // working results into the out params
          R = RR;
          return;
        }
        //else if (reduced == true)
        {
          int qRows = QQ.Length; int qCols = QQ[0].Length;
          int rRows = RR.Length; int rCols = RR[0].Length;
          // assumes m "gte" n !!

          // square-up R
          int dim = Math.Min(rRows, rCols);
          double[][] Rsquared = MatMake(dim, dim);
          for (int i = 0; i "lt" dim; ++i)
            for (int j = 0; j "lt" dim; ++j)
              Rsquared[i][j] = RR[i][j];

          // Q needs same number columns as R
          // so that inv(R) * trans(Q) works
          double[][] Qtrimmed = MatMake(qRows, dim);
          for (int i = 0; i "lt" qRows; ++i)
            for (int j = 0; j "lt" dim; ++j)
              Qtrimmed[i][j] = QQ[i][j];

          Q = Qtrimmed;
          R = Rsquared;
          return;
        }

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

        // nested helpers for MatDecomposeQR:
        // MatMake, MatIdentity, MatCopy, MatMult
        // VecNorm, VecDot, VecToMat

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

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

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

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

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

        static double[][] MatCopy(double[][] M)
        {
          int nr = M.Length; int nc = M[0].Length;
          double[][] result = MatMake(nr, nc);
          for (int i = 0; i "lt" nr; ++i)
            for (int j = 0; j "lt" nc; ++j)
              result[i][j] = M[i][j];
          return result;
        }

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

        static double[][] MatMult(double[][] matA,
          double[][] matB)
        {
          // aka MatProduct. diff name to avoid collision
          // with the one in primary MatPseudoInv()
          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 = MatMake(aRows, bCols);

          for (int i = 0; i "lt" aRows; ++i) // each row of A
            for (int j = 0; j "lt" bCols; ++j)
              for (int k = 0; k "lt" aCols; ++k)
                result[i][j] += matA[i][k] * matB[k][j];

          return result;
        }

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

        static double VecNorm(double[] vec)
        {
          int n = vec.Length;
          double sum = 0.0;
          for (int i = 0; i "lt" n; ++i)
            sum += vec[i] * vec[i];
          return Math.Sqrt(sum);
        }

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

        static double[][] VecToMat(double[] vec,
            int nRows, int nCols)
        {
          double[][] result = MatMake(nRows, nCols);
          int k = 0;
          for (int i = 0; i "lt" nRows; ++i)
            for (int j = 0; j "lt" nCols; ++j)
              result[i][j] = vec[k++];
          return result;
        }

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

        static double VecDot(double[] v1, double[] v2)
        {
          double result = 0.0;
          int n = v1.Length;
          for (int i = 0; i "lt" n; ++i)
            result += v1[i] * v2[i];
          return result;
        }
      } // MatDecomposeQR() primary helper

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

      static double[][] MatInverseUpperTri(double[][] U)
      {
        int n = U.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] * U[i][k];
            }
            result[j][k] /= U[k][k];
          }
        }
        return result;
      }

      static double[][] MatTranspose(double[][] m)
      {
        int nr = m.Length;
        int nc = m[0].Length;
        double[][] result = MatMake(nc, nr);  // note
        for (int i = 0; i "lt" nr; ++i)
          for (int j = 0; j "lt" nc; ++j)
            result[j][i] = m[i][j];
        return result;
      }

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

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

      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 = MatMake(aRows, bCols);

        for (int i = 0; i "lt" aRows; ++i) // each row of A
          for (int j = 0; j "lt" bCols; ++j)
            for (int k = 0; k "lt" aCols; ++k)
              result[i][j] += matA[i][k] * matB[k][j];

        return result;
      }

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

    } // MatPseudoInv()

  } // Program

} // ns