In the back of my mind, I wasn’t entirely happy with my current version of matrix QR decomposition using the Householder algorithm. The current version worked fine, but the underlying code just didn’t feel quite right in some way that I couldn’t articulate.

So I decided to take another stab at matrix decomposition using Householder. For me, the main purpose of QR decomposition is to compute the relaxed Moore-Penrose pseudo-inverse, which in turn is used to train a linear regression model or a quadratic regression model.

After a day or so of work, I got a new version of QR-Householder decomposition up and running to my satisfaction. The output of a demo run:

Begin QR decomposition using Householder algorithm
More precise than QR w/ modified Gram-Schmidt and QR w/
Givens but more complicated.
Latest version (July 2026)

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 validated my QR decomp implementation by sending the same input matrix to the Python language np.linalg.qr() function, to make sure the results were the same.

There are a ton of details. My implementation workks only for matrices that have more rows than columns — such as training data. My implementation returns a ‘reduced’ Q and a ‘reduced’ R (needed for pseudo-inverse) rather than ‘full’ Q and R matrices.

OK. Good fun. Next, I’ll need to implement a relaxed Moore-Penrose pseudo-inverse using my new QR decomp code. And then after that, I’ll need to use the pseudo-inverse to train a linear regression model.

I’m not very good at articulating subjective things. Here are two screen captures from a short AI-generated video titled “Oracle”, from a guy called Anglomangler. I can’t explain why or how, but the nightmarish quality of the video really makes an impact on me.

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 MatrixDecompQRHouseholder
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin QR decomposition using" +
        " Householder algorithm ");
      Console.WriteLine("More precise than QR w/ modified" +
        " Gram-Schmidt and QR w/ Givens but more " +
        "complicated. ");
      Console.WriteLine("Latest version (June 2026) ");
      
      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;
      QRHouseholder.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("");
      }
    }

  } // class Program

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

  public class QRHouseholder
  {
    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

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

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

  } // class QR_Householder

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

} // ns

Python language validation program:

import numpy as np

def main():
  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)

  # call NumPy function
  print("\nNumPy QR: ")
  Q, R = np.linalg.qr(A, mode='reduced')
  print("\nNumPy Q = "); print(Q)
  print("\nNumPy R = "); print(R)

if name == "__main__":
  main()