A lot of things require practice, practice, and more practice. Coding is one of these things.

To keep my numerical programming skills sharp, and to prepare for an upcoming training class I will present at my workplace, I decided to implement matrix QR decomposition from scratch. This is a very difficult problem. And, before I forget to mention it, my final implementation is mildly numerically unstable meaning that numerical round off errors can easily accumulate and end up giving a bad answer. So, my code is useful for lightweight tasks but not for most production use.

Unfortunately for me, once I get a technical challenge idea in my head, I can’t let it go. The result was that I spent an entire weekend working on the problem — I just couldn’t rest until I finished. But about 200 lines of code and quite a few hours later, I got a working example that I was quite pleased with.

Suppose you have a regular number, like x = 24. A decomposition of x is two numbers that multiply to give the original number. For example, a = 3 and b = 8 is a decomposition because a * b = x.

Matrix decomposition extends this idea to breaking down a matrix A into two matrices that multiply to give A. There are many kinds of matrix decomposition. In QR decomposition, A = Q * R where matrix Q is an orthogonal matrix (orthogonality is a big sub-topic in itself) and R is an upper triangular matrix (0s below the main diagonal).

Matrix decomposition doesn’t have any direct use, but decomposition is used by dozens of important data science and machine learning algorithms.

There are several algorithms for QR decomposition. The simplest (where “simplest” is relative — all the QR implementation algorithms are very complicated) is called the Gram–Schmidt algorithm. (The main alternative is called the Householder algorithm).

Left: The Wikipedia article on QR decomposition that I used as my reference. Right: I used Excel to convert the Wikipedia example results from fractions to decimals so I could compare them with results of my demo program.

I used the Wikipedia entry on QR decomposition as my source. The example given there is:

A =
     12    -51      4
      6    167    -68
     -4     24    -41

Q =
     6/7      21     -14
     3/7   158/175   -70
    -2/7     6/35   -33/35

  = (in decimal form)
     0.8571  -0.3943  -0.3314
     0.4286   0.9029   0.0343
    -0.2857   0.1714  -0.9429

R =
     14     21    -14
      0    175    -70
      0      0     35

Coding this demo was quite tricky. If it wasn’t for my obsessive personality when it comes to solving technical problems, I might have been tempted to give up.

My demo code is quite long but I’ve pasted all of it below.

Three images from an Internet search for “obsessive”. Left: By artist Georges Clairin (1843 – 1919). Center: By contemporary illustrator ‘Antartis’. Right: By artist Zhao Chun (b. 1970).

using System;
namespace MatLib
{
  class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin QR decomp demo");

      double[][] A = new double[3][];
      A[0] = new double[] { 12, -51, 4 };
      A[1] = new double[] { 6, 167, -68 };
      A[2] = new double[] { -4, 24, -41 };

      Console.WriteLine("\nOriginal matrix A =");
      MatShow(A, 4, 12);

      double[][] Q = null;
      double[][] R = null;

      int dummy = MatDecompQR(A, out Q, out R);
      Console.WriteLine("\nresult Q = ");
      MatShow(Q, 4, 12);
      Console.WriteLine("\nresult R = ");
      MatShow(R, 4, 12);

      double[][] qr = MatProduct(Q, R);
      Console.WriteLine("\nproduct of Q*R =");
      MatShow(qr, 4, 12);

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

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

    static int MatDecompQR(double[][] A,
      out double[][] Q, out double[][] R)
    {
      // work with rows of the transpose
      // return another transpose at end
      double[][] a = MatTranspose(A); 
      double[][] u = MatDuplicate(a);
      int rows = a.Length;  // of the transpose
      int cols = a[0].Length;

      Q = MatCreate(cols, rows);
      R = MatCreate(cols, rows);

      // first row of a (first col of M)
      for (int j = 0; j "less-than" cols; ++j)
        u[0][j] = a[0][j];

      double[] accum = new double[cols];
      // remaining rows of a
      for (int i = 1; i "less-than" rows; ++i)  
      {
        for (int j = 0; j "less-than" cols; ++j)
        {
          // accumulate projections
          accum = new double[cols];
          for (int t = 0; t "less-than" i; ++t)
          {
            double[] proj = VecProjection(u[t], a[i]);
            for (int k = 0; k "less-than" cols; ++k)
              accum[k] += proj[k];
          }
        }
        for (int k = 0; k "less-than" cols; ++k)
          u[i][k] = a[i][k] - accum[k];
      }

      for (int i = 0; i "less-than" rows; ++i)
      {
        double norm = VecNorm(u[i]);
        for (int j = 0; j "less-than" cols; ++j)
          u[i][j] = u[i][j] / norm;
      }
      // at this point u is Q(trans)

      double[][] q = MatTranspose(u);
      for (int i = 0; i "less-than" q.Length; ++i)
        for (int j = 0; j "less-than" q[0].Length; ++j)
          Q[i][j] = q[i][j];

      double[][] r = MatProduct(u, A);
      for (int i = 0; i "less-than" r.Length; ++i)
        for (int j = 0; j "less-than" r[0].Length; ++j)
          R[i][j] = r[i][j];

      return 0;
    } // 

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

    static double[][] MatDuplicate(double[][] m)
    {
      int rows = m.Length;
      int cols = m[0].Length;
      double[][] result = MatCreate(rows, cols);
      for (int i = 0; i "less-than" rows; ++i)
        for (int j = 0; j "less-than" cols; ++j)
          result[i][j] = m[i][j];
      return result;
    }

    static double VecDotProd(double[] u, double[] v)
    {
      double result = 0.0;
      for (int i = 0; i "less-than" u.Length; ++i)
        result += u[i] * v[i];
      return result;
    }

    static double[] VecProjection(double[] u, double[] a)
    {
      // proj(u, a) = (inner(u,a) / inner(u, u)) * u
      // u cannot be all 0s
      int n = u.Length;
      double dotUA = VecDotProd(u, a);
      double dotUU = VecDotProd(u, u);
      double[] result = new double[n];
      for (int i = 0; i "less-than" n; ++i)
        result[i] = (dotUA / dotUU) * u[i];
      return result;
    }

    static double[][] MatCreate(int rows, int cols)
    {
      double[][] result = new double[rows][];
      for (int i = 0; i "less-than" rows; ++i)
        result[i] = new double[cols];
      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 = MatCreate(aRows, bCols);

      for (int i = 0; i "less-than" aRows; ++i) 
        for (int j = 0; j "less-than" bCols; ++j)
          for (int k = 0; k "less-than" aCols; ++k) 
            result[i][j] += matA[i][k] * matB[k][j];

      return result;
    }

    static void MatShow(double[][] m, int dec, int wid)
    {
      for (int i = 0; i "less-than" m.Length; ++i)
      {
        for (int j = 0; j "less-than" m[0].Length; ++j)
        {
          double v = m[i][j];
          if (Math.Abs(v) "less-than" 1.0e-5)
            v = 0.0;  // avoid "-0.00"
          Console.Write(v.ToString("F" + dec).PadLeft(wid));
        }
        Console.WriteLine("");
      }
    }

    static double[][] MatTranspose(double[][] m)
    {
      int nr = m.Length;
      int nc = m[0].Length;
      double[][] result = MatCreate(nc, nr);  // note
      for (int i = 0; i "less-than" nr; ++i)
        for (int j = 0; j "less-than" nc; ++j)
          result[j][i] = m[i][j];
      return result;
    }
  } // Program
} // ns