If you have a matrix, there are several ways to decompose it: LUP decomposition, SVD decomposition, QR decomposition, Cholesky decomposition (only for positive-definite matrices). 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 evening, 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 Givens rotations algorithm, using the C# language.

The output of my demo:

Begin QR decomposition using Givens rotations 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 evening in the Pacific Northwest.

I’m a big fan of science fiction movies — good movies, bad movies, old movies, new movies — any kind of science fiction movies. Here are three that were OK (I give all three a “C” grade on my purely personal “A” to “F” scale) but could have been great.

Right: “John Carter” (2012) – Like many of my friends, we all grew up reading “A Princess of Mars” (1912) by ER Burroughs, and it is hands-down our favorite book of all time. When we heard that Disney was adapting the book to a movie, we were all excited. And then . . disaster. The main problem is the actress who plays Princess Dejah Thoris. Instead of a princess, she comes off as an annoying, abasive, tatoo-covered, entitled girl-thug. She single-handedly destroyed the movie.

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 MatrixDecompQRGivens
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin QR decomposition using" +
        " Givens rotations 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 Givens rotations
      // '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 ");
      int k = Math.Min(m, n);

      double[][] QQ = new double[m][];  // working Q mxm
      for (int i = 0; i "lt" m; ++i)
        QQ[i] = new double[m];
      for (int i = 0; i "lt" m; ++i)
        QQ[i][i] = 1.0;  // Identity

      double[][] RR = new double[m][];  // working R mxn
      for (int i = 0; i "lt" m; ++i)
        RR[i] = new double[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

      for (int j = 0; j "lt" k; ++j) // outer loop
      {
        for (int i = m-1; i "gt" j; --i) // inner loop
        {
          double a = RR[j][j];
          double b = RR[i][j];
          //if (b == 0.0) continue;  // risky
          if (Math.Abs(b) "lt" 1.0e-12) continue;

          double r = Math.Sqrt((a * a) + (b * b));
          double c = a / r;
          double s = -b / r;

          for (int col = j; col "lt" n; ++col)
          {
            double tmp = 
              (c * RR[j][col]) - (s * RR[i][col]);
            RR[i][col] = 
              (s * RR[j][col]) + (c * RR[i][col]);
            RR[j][col] = tmp;
          }

          for (int row = 0; row "lt" m; ++row)
          {
            double tmp = 
              (c * QQ[row][j]) - (s * QQ[row][i]);
            QQ[row][i] = 
              (s * QQ[row][j]) + (c * QQ[row][i]);
            QQ[row][j] = tmp;
          }
        } // inner i loop

      } // j loop

      // Q_reduced = Q[:, :k] // all rows, col 0 to k
      // R_reduced = R[:k, :] // rows 0 to k, all cols

      double[][] Qred = new double[m][];
      for (int i = 0; i "lt" m; ++i)
        Qred[i] = new double[k];
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" k; ++j)
          Qred[i][j] = QQ[i][j];
      Q = Qred;

      double[][] Rred = new double[k][];
      for (int i = 0; i "lt" k; ++i)
        Rred[i] = new double[n];
      for (int i = 0; i "lt" k; ++i)
        for (int j = 0; j "lt" n; ++j)
          Rred[i][j] = RR[i][j];
      R = Rred;

    } // MatDecompQR()

  } // Program

} // ns