I recently (yesterday) made major revisions to my personal C# matrix QR decomposition using the Householder algorithm. In my work environment, I use matrix QR decomposition for matrix pseudo-inverse of a matrix of training data. So, my next step after a new QR decomposition function, presented in yesterday’s blog post, was to implement a new pseudo-inverse function.

Suppose matrix X holds training data (typically a design matrix, where a leading column of 1.0 values has been added to deal with the model bias). Then, to find the pseudo-inverse of X, pinv(X):

1. X = Q * R
2. pinv(X) = inv(Q * R)
3. pinv(X) = inv(Q) * inv(R)
4. pinv(X) = tr(Q) * inv(R)

In step 1, the matrix of training data X is decomposed using QR-Householder decomposition. The QR decomposition applies to any shape matrix, but regular inverse applies only to square matrices. Because training data almost always has more rows than columns, it’s possible to use a relaxed form that assumes this condition.

In step 2, the matrix inverse is applied to both sides of the equation.

In step 3, the matrix property inv(A*B) = inv(B) * inv(A) is applied.

In step 4, the inverse of Q is the transpose of Q, which is very easy to calculate. This is a special property of the Q result of QR decomposition, and is why QR decomposition is used. The R matrix is upper-triangular, and that inverse is a special case which is easy to compute.

Expressed in code:

public static double[][] MatPinv(double[][] M)
{
  double[][] Q; double[][] R;
  MatDecompQR(M, out Q, out R);  // Householder
  double[][] Ri = MatInvUpperTri(R);
  double[][] Qi = MatTranspose(Q);
  double[][] result = MatProduct(Ri, Qi);
  return result;
}

Most of the work is done in the MatDecompQR() function.

Here’s output of a demo run:

Relaxed  MP pseudo-inverse using QR decomposition
 Householder algorithm

Source (tall) matrix A:
   1.0   2.0   3.0   4.0   5.0
   0.0  -3.0   5.0  -7.0   9.0
   2.0   0.0  -2.0   0.0  -2.0
   4.0  -1.0   5.0   6.0   1.0
   3.0   6.0   8.0   2.0   2.0
   5.0  -2.0   4.0  -4.0   3.0

Computing pinv
Done

pseudo-inverse =
   0.0911  -0.0561   0.1650  -0.0250  -0.0144   0.1445
   0.0921  -0.0417   0.0755  -0.1379   0.0939   0.0053
  -0.1618   0.0430  -0.1421   0.1041   0.0790  -0.0414
   0.0664  -0.0250  -0.0145   0.0753  -0.0372  -0.0454
   0.1808   0.0376   0.0624  -0.0541  -0.0470   0.0104

End demo

I partially validated my relaxed MP pseudo-inverse by feeding the same source matrix to the NumPy np.linalg.pinv() function to make sure I got the same result.

OK. Good fun. My next step will be to use this new relaxed MP pseudo-inverse via QR-Householder functionality to train a machine learning regression problem (probably quadratic regression or maybe linear regression).

I don’t collect things. I collect ideas, such as a new algorithm to compute a pseudo-inverse.

Here is card #1.

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

using System;
using System.IO;

namespace MatrixPseudoInverseQRHouseholder
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nRelaxed  MP pseudo-inverse " +
        "using QR decomposition Householder algorithm ");
            
      double[][] A = new double[6][];
      A[0] = new double[] { 1, 2, 3, 4, 5 };
      A[1] = new double[] { 0, -3, 5, -7, 9 };
      A[2] = new double[] { 2, 0, -2, 0, -2 };
      A[3] = new double[] { 4, -1, 5, 6, 1 };
      A[4] = new double[] { 3, 6, 8, 2, 2 };
      A[5] = new double[] { 5, -2, 4, -4, 3 };

      Console.WriteLine("\nSource (tall) matrix A: ");
      MatShow(A, 1, 6);

      Console.WriteLine("\nComputing pinv ");
      double[][] Pinv = QRHouseholder.MatPinv(A);
      Console.WriteLine("Done ");

      Console.WriteLine("\npseudo-inverse = ");
      MatShow(Pinv, 4, 9);
   
      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main

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

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

  } // class Program

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

  public class QRHouseholder
  {
    // container class for MP pseudo-inverse via Householder
    // A = Q * R
    // pinv(A) = inv(R) * inv(Q)  note order matters
    //         = inv upper tri (easy) * transpose (easy)

    public static double[][] MatPinv(double[][] M)
    {
      double[][] Q; double[][] R;
      MatDecompQR(M, out Q, out R);  // Householder
      double[][] Ri = MatInvUpperTri(R);
      double[][] Qi = MatTranspose(Q);
      double[][] result = MatProduct(Ri, Qi);
      return result;
    }

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

    public static double[][] MatInvUpperTri(double[][] U)
    {
      int n = U.Length;  // must be square matrix

      double[][] result = MatMake(n, n);
      for (int i = 0; i "lt" n; ++i)
        result[i][i] = 1.0;
      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;
    }

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

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

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

    public static double[][] MatTranspose(double[][] M)
    {
      int nRows = M.Length;
      int nCols = M[0].Length;
      double[][] result = MatMake(nCols, nRows);
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[j][i] = M[i][j];
      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 MatDecompQR(double[][] A, 
      out double[][] Q,  out double[][] R)
    {
      int m = A.Length; int n = A[0].Length;
      if (m "lt" n)
        Console.WriteLine("FATAL: nRows must be gte nCols");

      double[][] QQ = MatMake(m, m); // working full Q
      for (int i = 0; i "lt" m; ++i)
        QQ[i][i] = 1.0;  // identity matrix

      double[][] RR = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          RR[i][j] = A[i][j]; // copy of A is working R

      int k = Math.Min(m, n);  // or just use n
      for (int j = 0; j "lt" k; ++j) // main processing loop
      {
        int xn = m - j;
        double[] x = new double[xn];
        for (int i = 0; i "lt" xn; ++i)
          x[i] = RR[j + i][j];

        double ss = 0.0;
        for (int i = 0; i "lt" xn; ++i)
          ss += x[i] * x[i];
        double normX = Math.Sqrt(ss);

        // if (normX == 0.0) continue;
        if (Math.Abs(normX) "lt" 1.0e-12) continue;

        double sign;
        if (x[0] "gte" 0.0) sign = -1.0;
        else sign = 1.0; // counter-intuitive
      
        double[] u = new double[xn];
        for (int i = 0; i "lt" xn; ++i)
          u[i] = x[i] / (x[0] - sign * normX); // check div 0
        u[0] = 1.0;

        // compute scaling factor tau = 2 / (u^T * u)
        double tau = -sign * (x[0] - sign * normX) / normX;

        // dimensions for sub-matrices
        int nRowsSubR = m - j;   int nColsSubR = n - j;
        int nRowsSubQ = m;       int nColsSubQ = m - j;

        double[] vr = new double[nColsSubR];
        for (int c = 0; c "lt" nColsSubR; ++c)
        {
          double acc = 0.0;
          for (int r = 0; r "lt" nRowsSubR; ++r)
            acc += u[r] * RR[j + r][j + c];
          vr[c] = acc;
        }

        double[] vq = new double[nRowsSubQ];
        for (int r = 0; r "lt" nRowsSubQ; ++r)
        {
          double acc = 0.0;
          for (int c = 0; c "lt" nColsSubQ; ++c)
            acc += u[c] * QQ[r][j + c];
          vq[r] = acc;
        }

        // update sub-R
        for (int r = 0; r "lt" nRowsSubR; ++r)
          for (int c = 0; c "lt" nColsSubR; ++c)
            RR[j + r][j + c] -= tau * u[r] * vr[c];

        // update sub-Q
        for (int r = 0; r "lt" nRowsSubQ; ++r)
          for (int c = 0; c "lt" nColsSubQ; ++c)
            QQ[r][j + c] -= tau * vq[r] * u[c];
       
      } // j

      // extract QQ RR into out params
      Q = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          Q[i][j] = QQ[i][j];

      R = MatMake(n, n);
      for (int i = 0; i "lt" n; ++i)
        for (int j = 0; j "lt" n; ++j)
          R[i][j] = RR[i][j];

      return;
    } // MatDecompQR

  } // class QRHouseholder

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

} // ns

Python validation code:

  import numpy as np

  A = np.array([
    [ 1, 2, 3, 4, 5 ],
    [ 0, -3, 5, -7, 9 ],
    [ 2, 0, -2, 0, -2 ],
    [ 4, -1, 5, 6, 1 ],
    [ 3, 6, 8, 2, 2 ],
    [ 5, -2, 4, -4, 3 ]], dtype=np.float64)

  print("\nA = "); print(A)

  Pinv = np.linalg.pinv(A)
  print("\nNumPy pseudo-inverse = "); print(Pinv)