One early morning, at about 3:00 AM, I figured I’d take a look at an idea that had been percolating in my mind for a while. I have two different C# implementations of a function to compute the QR decomposition of a matrix. Both use the Householder algorithm, as opposed to the Gram-Schmidt algorithm, or the Givens algorithm.
The first QR-Householder implementation is relatively simpler (but still extremely complex) but has accuracy to approximately 1.0e-8. I use it for several machine learning functions (Moore-Penrose pseudo-inverse in particular).
The second QR-Householder implementation is relatively more complicated (and extremely complex) and has accuracy to approximately 1.0e-12. I use it primarily as a helper function for the extraordinarily complicated SVD (singular value decomposition) decomposition.
Both QR-Householder versions work only for matrices that have number of rows greater than or equal to number of columns (standard and adequate for nearly all machine learning scenarios).
The first version has a ‘reduced’ parameter to allow full or reduced results. But in machine learning, the reduced result is always used. The second version computes only a reduced result because that’s required for SVD.
The first QR-Householder medium-accuracy, medium-complexity version has 91 lines of code and 7 helper functions. The second higher-accuracy, higher-complexity version has 50 lines of code but 16 helper functions.
The first version looks like:
public class MatDecomQR1 // older version based on rosetta
{
public static void MatDecomposeQR(double[][] M,
out double[][] Q, out double[][] R, bool reduced)
{
. . . // 91 LOC
}
// ------------------------------------------------------
// 7 helpers:
private static double[][] MatMake(int nRows,
int nCols) { . . }
private static double[][] MatIdentity(int n) { . . }
private static double[][] MatCopy(double[][] M) { . . }
private static double[][] MatProduct(double[][] A,
double[][] B) { . . }
private static double[][] VecToMat(double[] vec,
int nRows, int nCols) { . . }
private static double VecNorm(double[] vec) { . . }
private static double VecDot(double[] v1,
double[] v2) { . . }
}
The second version looks like:
public class MatDecomQR2 // newer version for SVD
{
public static void MatDecompQR(double[][] A,
out double[][] Q, out double[][] R)
{
. . . // 50 LOC
}
// ------------------------------------------------------
// 16 helpers:
private static double[][] MatMake(int nRows,
int nCols) { . . }
private static double[][] MatIdentity(int n) { . . }
private static double[][] MatCopy(double[][] M) { . . }
private static double[] MatColToEnd(double[][] A,
int k) { . . }
private static double[][] MatRowColToEnd(double[][] A,
int k) { . . }
private static double[] VecCopy(double[] v) { . . }
private static double[] VecMatProduct(double[] v,
double[][] A) { . . }
private static double[][] VecOuter(double[] v1,
double[] v2) { . . }
private static double[][] MatScalarMult(double x,
double[][] A) { . . }
private static void MatSubtractCorner(double[][] A,
int k, double[][] C) { . . }
private static double[][]
MatExtractAllRowsColToEnd(double[][] A, int k) { . . }
private static double[] MatVecProduct(double[][] A,
double[] v) { . . }
private static void MatSubtractRows(double[][] A,
int k, double[][] C) { . . }
private static double[][]
MatExtractFirstCols(double[][] A, int k) { . . }
private static double[][]
MatExtractFirstRows(double[][] A, int k) { . . }
private static double VecNorm(double[] vec) { . . }
} // MatDecomQR2
After exploring both versions, I didn’t come to any surprising conclusions. The first QR-Householder implementation is suitable for general purpose machine learning uses. The second QR-Householder implementation is suitable as a helper for SVD.
Comparing two different versions of QR decomposition was an interesting exploration. Researchers have compared the average IQ of people in different countries. Left: The United States has an overage average of approximately 100. Right: But there is a lot of variation in IQ in the U.S. Here is a woman who didn’t fully grasp the fact that she is subject to traffic laws like everyone else. If you compute the IQ in the United States without African Americans, the average IQ increases to about 104.
Demo program. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols (my blog editor chokes on symbols).
using System;
using System.IO;
namespace MatDecompQR
{
internal class Program
{
static void Main(string[] args)
{
Console.WriteLine("\nBegin matrix QR decomposition" +
" Householder variations ");
Console.WriteLine("\nKey takeaway: standard QR" +
" decomp algorithms work only for nRows gte nCols ");
double[][] Q;
double[][] R;
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);
double[][] B = new double[3][];
B[0] = new double[] { 3, 1, 1, 0, 5 };
B[1] = new double[] { -1, 3, 1, -2, 4 };
B[2] = new double[] { 0, 2, 2, 1, -3 };
Console.WriteLine("\nSource (wide) matrix B: ");
MatShow(B, 1, 6);
Console.WriteLine("\nComputing QR for A using" +
" older rosetta version ");
MatDecomQR1.MatDecomposeQR(A, out Q, out R,
reduced: true);
Console.WriteLine("\nQ = ");
MatShow(Q, 6, 11);
Console.WriteLine("\nR = ");
MatShow(R, 6, 11);
//// the rosetta version fails for wide matrices
//Console.WriteLine("\nThe old rosetta version" +
// " fails for wide matrices ");
Console.WriteLine("\n=========================== ");
Console.WriteLine("\nComputing QR for A using newer" +
" hybrid version for SVD ");
MatDecomQR2.MatDecompQR(A, out Q, out R);
Console.WriteLine("\nQ = ");
MatShow(Q, 6, 11);
Console.WriteLine("\nR = ");
MatShow(R, 6, 11);
//Console.WriteLine("\nComputing QR for B using newer" +
// " AI hybrid version ");
//MatDecomQR2.MatDecompQR(B, out Q, out R);
//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("");
}
}
} // Program
public class MatDecomQR1 // older version based on rosetta
{
public 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("No 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;
}
} // MatDecomposeQR()
// ------------------------------------------------------
private 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;
}
// ------------------------------------------------------
private 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;
}
// ------------------------------------------------------
private 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;
}
// ------------------------------------------------------
private 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);
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;
}
// ------------------------------------------------------
private 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;
}
// ------------------------------------------------------
private 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);
}
// ------------------------------------------------------
private 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;
}
// ------------------------------------------------------
} // MatDecomQR1
public class MatDecomQR2 // newer hybrid version from AI
{
public static void MatDecompQR(double[][] A,
out double[][] Q, out double[][] R)
{
// reduced QR decomp using Householder reflections
// this version used by SVD
int m = A.Length; int n = A[0].Length;
//if (m "lt" n)
// throw new Exception("No rows less than cols");
double[][] QQ = MatIdentity(m); // working Q
double[][] RR = MatCopy(A); // working R
for (int k = 0; k "lt" n; ++k)
{
double[] x = MatColToEnd(RR, k); // RR[k:, k]
double normx = VecNorm(x);
if (Math.Abs(normx) "lt" 1.0e-12) continue;
double[] v = VecCopy(x);
double sign;
if (x[0] "gte" 0.0) sign = 1.0; else sign = -1.0;
v[0] += sign * normx;
double normv = VecNorm(v);
for (int i = 0; i "lt" v.Length; ++i)
v[i] /= normv;
// apply reflection to R
double[][] tmp1 =
MatRowColToEnd(RR, k); // RR[k:, k:]
double[] tmp2 = VecMatProduct(v, tmp1);
double[][] tmp3 = VecOuter(v, tmp2);
double[][] B = MatScalarMult(2.0, tmp3);
MatSubtractCorner(RR, k, B); // R[k:, k:] -= B
// accumulate Q
double[][] tmp4 =
MatExtractAllRowsColToEnd(QQ, k); // QQ[:, k:]
double[] tmp5 = MatVecProduct(tmp4, v);
double[][] tmp6 = VecOuter(tmp5, v);
double[][] C = MatScalarMult(2.0, tmp6);
MatSubtractRows(QQ, k, C); // QQ[:, k:] -= C
} // for k
// return parts of QQ and RR
Q = MatExtractFirstCols(QQ, n); // Q[:, :n]
R = MatExtractFirstRows(RR, n); // R[:n, :]
} // MatDecompQR()
// ------------------------------------------------------
private 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;
}
// ------------------------------------------------------
private 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;
}
// ------------------------------------------------------
private 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;
}
// ------------------------------------------------------
private static double[] MatColToEnd(double[][] A,
int k)
{
// last part of col k, from [k,k] to end
int m = A.Length; int n = A[0].Length;
double[] result = new double[m - k];
for (int i = 0; i "lt" m - k; ++i)
result[i] = A[i + k][k];
return result;
}
// ------------------------------------------------------
private static double[][] MatRowColToEnd(double[][] A,
int k)
{
// block of A, row k to end, col k to end
// A[k:, k:]
int m = A.Length; int n = A[0].Length;
double[][] result = MatMake(m - k, n - k);
for (int i = 0; i "lt" m - k; ++i)
for (int j = 0; j "lt" n - k; ++j)
result[i][j] = A[i + k][j + k];
return result;
}
// ------------------------------------------------------
private static double[] VecCopy(double[] v)
{
double[] result = new double[v.Length];
for (int i = 0; i "lt" v.Length; ++i)
result[i] = v[i];
return result;
}
// ------------------------------------------------------
private static double[] VecMatProduct(double[] v,
double[][] A)
{
// v * A
int m = A.Length; int n = A[0].Length;
double[] result = new double[n];
for (int j = 0; j "lt" n; ++j)
for (int i = 0; i "lt" m; ++i)
result[j] += v[i] * A[i][j];
return result;
}
// ------------------------------------------------------
private static double[][] VecOuter(double[] v1,
double[] v2)
{
// vector outer product
double[][] result = MatMake(v1.Length, v2.Length);
for (int i = 0; i "lt" v1.Length; ++i)
for (int j = 0; j "lt" v2.Length; ++j)
result[i][j] = v1[i] * v2[j];
return result;
}
// ------------------------------------------------------
private static double[][] MatScalarMult(double x,
double[][] A)
{
// x * A element-wise
int m = A.Length; int n = A[0].Length;
double[][] result = MatMake(m, n);
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" n; ++j)
result[i][j] = x * A[i][j];
return result;
}
// ------------------------------------------------------
private static void MatSubtractCorner(double[][] A,
int k, double[][] C)
{
// subtract elements in C from lower right A
// A[k:, k:] -= C
int m = A.Length; int n = A[0].Length;
for (int i = 0; i "lt" m - k; ++i)
for (int j = 0; j "lt" n - k; ++j)
A[i + k][j + k] -= C[i][j];
return; // modify A in-place
}
// ------------------------------------------------------
private static double[][]
MatExtractAllRowsColToEnd(double[][] A, int k)
{
// A[:, k:]
int m = A.Length; int n = A[0].Length;
double[][] result = MatMake(m, n - k);
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" n - k; ++j)
result[i][j] = A[i][j + k];
return result;
}
// -----------------------------------------------------
private static double[] MatVecProduct(double[][] A,
double[] v)
{
// A * v
int m = A.Length; int n = A[0].Length;
double[] result = new double[m];
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" n; ++j)
result[i] += A[i][j] * v[j];
return result;
}
// ------------------------------------------------------
private static void MatSubtractRows(double[][] A,
int k, double[][] C)
{
// A[:, k:] -= C
int m = A.Length; int n = A[0].Length;
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" n - k; ++j)
A[i][j + k] -= C[i][j];
return; // modify A in-place
}
// ------------------------------------------------------
private static double[][]
MatExtractFirstCols(double[][] A, int k)
{
// A[:, :n]
int m = A.Length; int n = A[0].Length;
double[][] result = MatMake(m, k);
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" k; ++j)
result[i][j] = A[i][j];
return result;
}
// ------------------------------------------------------
private static double[][]
MatExtractFirstRows(double[][] A, int k)
{
// A[:n, :]
int m = A.Length; int n = A[0].Length;
double[][] result = MatMake(k, n);
for (int i = 0; i "lt" k; ++i)
for (int j = 0; j "lt" n; ++j)
result[i][j] = A[i][j];
return result;
}
// ------------------------------------------------------
private 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);
}
// ------------------------------------------------------
} // MatDecomQR2
} // ns
.NET Test Automation Recipes
Software Testing
SciPy Programming Succinctly
Keras Succinctly
R Programming
2026 Visual Studio Live
2025 Summer MLADS Conference
2026 DevIntersection Conference
2025 Machine Learning Week
2025 Ai4 Conference
2026 G2E Conference
2026 iSC West Conference
You must be logged in to post a comment.