If you have a matrix A and apply QR decomposition, you get matrices Q and R, where Q * R = A. One common use of QR decomposition is to find the inverse of A. This assumes A is a square matrix, in which case it’s possible to write a simplified version of QR decomposition.

I decided to implement a version of QR decomposition using C#, where the source matrix A is not necessarily square. This is useful to implement a pseudo-inverse function. Before I go any further, let me point out that there are many very subtle details that this blog post doesn’t address.

There are several algorithms to implement QR decomposition. The two that I come across most often are the Householder algorithm and the Gram-Schmidt algorithm. I prefer the Householder algorithm — much simpler to implement in my opinion.

One of the many tricky details with QR decomposition is that it depends entirely on whether the number of rows in the source matrix is less than the number of columns, or the number rows is greater than or equal number columns. In machine learning, it’s far more common for number of rows to be greater than or equal to number columns (think training data). This geometry is sometimes called a “tall and skinny” matrix.

And yet another tricky detail is that QR decomposition can return “complete” versions of Q and R, or return “reduced” versions. The idea is best explained by the output of my demo:

Begin QR decomposition 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

Applying QR (complete)

Q:
  -0.3636   0.5998   0.6426  -0.1606  -0.2634
  -0.5455  -0.2854   0.2952   0.5362   0.4964
  -0.7273  -0.2677  -0.3760  -0.4394  -0.2549
  -0.1818   0.4692  -0.5091   0.6279  -0.3055
  -0.0909   0.5167  -0.3153  -0.3152   0.7252

R:
 -11.0000  -4.6364 -10.0000
  -0.0000   8.8603   2.2162
  -0.0000  -0.0000  -6.1716
  -0.0000  -0.0000   0.0000
  -0.0000  -0.0000   0.0000

Applying QR (reduced)

Q:
  -0.3636   0.5998   0.6426
  -0.5455  -0.2854   0.2952
  -0.7273  -0.2677  -0.3760
  -0.1818   0.4692  -0.5091
  -0.0909   0.5167  -0.3153

R:
 -11.0000  -4.6364 -10.0000
  -0.0000   8.8603   2.2162
  -0.0000  -0.0000  -6.1716

End demo

The two result types, complete or reduced, are needed by different algorithms. For example, when computing a pseudo-inverse, the reduced results are needed.

Anyway, a fun exploration for me.

There are many different algorithms for matrix decomposition. I’m a big fan of 1950s science fiction movies. Several of these movies featured ray guns that could decompose the molecular structure of their target.

Right: In “Earth vs. the Flying Saucers” (1956), aliens with formidable ray weapons arrive. They demand surrender. Earth fights back and devises a sonic weapon to defeat the aliens.

Demo code. 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 QR decomposition C# ");

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

      double[][] Q; double[][] R;

      Console.WriteLine("\nApplying QR (complete) ");
      MatDecomposeQR(A, out Q, out R, reduced: false);
      Console.WriteLine("\nQ: ");
      MatShow(Q, 4, 9);
      Console.WriteLine("\nR: ");
      MatShow(R, 4, 9);

      Console.WriteLine("\nApplying QR (reduced) ");
      MatDecomposeQR(A, out Q, out R, reduced: true);
      Console.WriteLine("\nQ: ");
      MatShow(Q, 4, 9);
      Console.WriteLine("\nR: ");
      MatShow(R, 4, 9);

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

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

    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[][] 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);
      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;
    }

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

    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 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 = MatProduct(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 = MatProduct(QQ, H);
        RR = MatProduct(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: MatMake, MatIdentity, MatCopy,
      // 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 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()

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

  } // Program

} // ns