If you have a matrix, there are several ways to decompose it: LUP decomposition, SVD decomposition, QR decomposition, Cholesky decomposition (only for positive-definite). The QR decomposition can be used for many algorithms, such as computing a matrix inverse, or pseudo-inverse.

If you have a matrix A, the QR decomposition gives you a matrix Q and a matrix R such that A = Q * R (where the * is matrix multiplication).

There are three main algorithms to compute a QR decomposition: Householder, modified Gram-Schmidt, Givens. Each of these three algorithms has several variations. Householder is by far the most complicated of the three, but gives the best answer precision. Modified Gram-Schmidt is dramatically simpler than Householder, but not as precise for pathological matrices. Givens is roughly comparable to modified Gram-Schmidt in terms of complexity and precision.

One morning, I was waiting for the rain to stop so I could walk my dogs Kevin and Riley. I figured I’d make use of the time by implementing QR decomposition using the modified Gram-Schmidt algorithm, using the C# language.

The output of my demo:

Begin QR decomposition using modified Gram-Schmidt with C#
Less precise than QR-Householder but dramatically simpler

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 QR
Done

Q =
   0.134840   0.258894   0.147094   0.310903   0.869379
   0.000000  -0.410745   0.759588  -0.274531   0.180705
   0.269680  -0.029872  -0.526675  -0.219131   0.300339
   0.539360  -0.196660   0.116670   0.738280  -0.260150
   0.404520   0.776682   0.316260  -0.255197  -0.226191
   0.674200  -0.348511  -0.101841  -0.412033   0.049823

R =
   7.416198   0.809040   8.494918   1.887760   3.505839
   0.000000   7.303797   2.618812   5.678242  -2.031322
   0.000000   0.000000   7.998637  -2.988836   9.068775
   0.000000   0.000000   0.000000   8.732743  -1.486219
   0.000000   0.000000   0.000000   0.000000   4.809501

End demo

I verified the results by using the NumPy np.linalg.qr() library function, and I got the same results. The demo implementation returns “reduced” Q and R, which is standard. The demo implementation works only for “tall” matrices where the number of rows is greater than or equal to the number of columns, which is standard for machine learning scenarios.

Good fun on a rainy morning in the Pacific Northwest.

One of the advantages of implementing a function from scratch instead of using a library function is that from-scratch allows you complete visibility into the function.

Right: “The Outer Limits” TV series ran two seasons from 1963 to 1965. In Season 2, Episode 17 (the last episode of the series) “The Probe”, four plane crash survivors find themselves in an automated alien space probe. They are placed in disinfecting tubes. Decent episode (A- grade). One of the best sci-fi TV series of all time (solid A grade).

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

using System;
using System.IO;

namespace MatrixDecompQRGramSchmidt
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin QR decomposition using" +
        " modified Gram-Schmidt with C# ");
      Console.WriteLine("Less precise than QR-Householder" +
        " but dramatically simpler ");
      
      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 QR ");
      double[][] Q;
      double[][] R;
      MatDecompQR(A, out Q, out R);
      Console.WriteLine("Done ");

      Console.WriteLine("\nQ = ");
      MatShow(Q, 6, 11);
      Console.WriteLine("\nR = ");
      MatShow(R, 6, 11);

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

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

    static void MatDecompQR(double[][] A, out double[][] Q,
      out double[][] R)
    {
      // QR decomposition using 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;
    }

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

  } // Program

} // ns