One evening, after taking my two dogs out for a walk, I decided to tidy up my C# implementation of matrix inverse using QR decomposition. My main change was to refactor the MatInverseQR() function so that all the helper functions (10 of them) were nested inside, rather than existing in a utility library.
Computing the inverse of a matrix is extremely difficult. There are many different algorithms (QR decomposition, SVD decomposition, LUP decomposition, Cholesky decomposition, Gauss-Jordan elimination, etc.) and each algorithm has several variations. Three common variations for QR decomposition are Gram-Schmidt, Givens, and Householder. My matrix inverse uses QR-Householder.
My demo sets up a small 4-by-4 matrix A:
4.00 7.00 1.00 2.00 6.00 0.00 3.00 5.00 8.00 1.00 9.00 2.00 2.00 5.00 6.00 -3.00
The demo computes the inverse Ai:
0.5735 -1.2426 1.0221 -1.0074 0.0000 0.2500 -0.2500 0.2500 -0.4118 0.8088 -0.5735 0.6912 -0.4412 1.2059 -0.8824 0.7941
The demo concludes by computing Ai * A to verify the result is the Identity matrix:
1.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 -0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 1.0000
The -0.0000 at cell [1][3] indicates the value there is a very small negative number, not a special negative-zero value of some sort.
Good fun!
The films of movie director Alfred Hitchcock (1899-1980) often featured plot twists/inversions. Here are two of my favorites.
Left: In “To Catch a Thief” (1955), retired jewel thief John Robie (actor Cary Grant) is framed by an unknown criminal for thefts at a hotel on the French Riviera. Robie runs into several people, including a rich American widow Jessie and her daughter Frances (actress Grace Kelly), and the hotel wine steward Foussard and his pretty daughter Danielle. Robie discovers Danielle is the thief.
Right: In “North by Northwest” (1959), advertising executive Roger Thornhill (actor Cary Grant again) is mistaken by enemy agents as a U.S. spy. The enemy agents botch an attempt to kill Thornhill but frame him for a murder and so everyone is now after him. Thornhill meets pretty Eve Kendall (actress Eva Marie Saint) on a train. It turns out that Eve is actually the U.S. spy and the enemy agents are captured or killed in a chase scene at Mount Rushmore.
Demo code. Replace “lt” (less than), “gt”, “lte”, “gte”, “and” with Boolean operator symbols. My lame blog editor sometimes chokes on symbols.
using System;
using System.IO;
namespace MatrixInverseQR
{
internal class MatrixInverseQRProgram
{
static void Main(string[] args)
{
Console.WriteLine("\nBegin matrix inverse using QR ");
double[][] A = new double[4][];
A[0] = new double[] { 4.0, 7.0, 1.0, 2.0 };
A[1] = new double[] { 6.0, 0.0, 3.0, 5.0 };
A[2] = new double[] { 8.0, 1.0, 9.0, 2.0 };
A[3] = new double[] { 2.0, 5.0, 6.0, -3.0 };
Console.WriteLine("\nSource matrix A: ");
MatShow(A, 2, 6);
double[][] Ai = MatInverseQR(A);
Console.WriteLine("\nInverse matrix Ai: ");
MatShow(Ai, 4, 9);
double[][] AiA = MatProd(Ai, A);
Console.WriteLine("\nAi * A: ");
MatShow(AiA, 4, 9);
Console.WriteLine("\nEnd demo ");
Console.ReadLine();
} // Main
static double[][] MatProd(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 = MatMake(aRows, bCols);
double[][] result = new double[aRows][];
for (int i = 0; i "lt" aRows; ++i)
result[i] = new double[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] += matA[i][k] * matB[k][j];
return result;
}
// ======================================================
static double[][] MatInverseQR(double[][] M)
{
// inverse using QR decomp (Householder algorithm)
double[][] Q;
double[][] R;
MatDecomposeQR(M, out Q, out R);
// TODO: check if determinant is zero (no inverse)
// double absDet = 1.0;
// for (int i = 0; i "lt" M.Length; ++i)
// absDet *= R[i][i];
double[][] Rinv = MatInverseUpperTri(R); // std algo
double[][] Qtrans = MatTranspose(Q); // is inv(Q)
return MatProduct(Rinv, Qtrans);
// ----------------------------------------------------
// helpers: MatDecomposeQR, MatMake, MatTranspose,
// MatInverseUpperTri, MatIdentity, MatCopy, VecNorm,
// VecDot, VecToMat, MatProduct
// ----------------------------------------------------
static void MatDecomposeQR(double[][] mat,
out double[][] q, out double[][] r)
{
// QR decomposition, Householder algorithm.
// https://rosettacode.org/wiki/QR_decomposition
int m = mat.Length; // assumes mat is nxn
int n = mat[0].Length; // check m == n
double[][] Q = MatIdentity(m);
double[][] R = MatCopy(mat);
int end = n-1;
// 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[n - i];
int k = 0;
for (int ii = i; ii "lt" n; ++ii)
a[k++] = R[ii][i];
double normA = VecNorm(a);
if (a[0] "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;
double[][] h = MatIdentity(a.Length);
double vvDot = VecDot(v, v);
double[][] alpha = VecToMat(v, v.Length, 1);
double[][] beta = VecToMat(v, 1, v.Length);
double[][] aMultB = MatProduct(alpha, beta);
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) * aMultB[ii][jj];
// copy h into lower right of H
int d = n - h.Length;
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];
Q = MatProduct(Q, H);
R = MatProduct(H, R);
} // i
q = Q;
r = R;
} // QR decomposition
static double[][] MatInverseUpperTri(double[][] U)
{
int n = U.Length; // must be square matrix
double[][] result = MatIdentity(n);
for (int k = 0; k "lt" n; ++k)
{
for (int j = 0; j "lt" n; ++j)
{
for (int i = 0; i "lt" k; ++i)
{
result[j][k] -= result[j][i] * U[i][k];
}
result[j][k] /= U[k][k];
}
}
return result;
}
static double[][] MatTranspose(double[][] m)
{
int nr = m.Length;
int nc = m[0].Length;
double[][] result = MatMake(nc, nr); // note
for (int i = 0; i "lt" nr; ++i)
for (int j = 0; j "lt" nc; ++j)
result[j][i] = m[i][j];
return result;
}
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;
}
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;
}
static double[][] MatCopy(double[][] m)
{
int nRows = m.Length; int nCols = m[0].Length;
double[][] result = MatMake(nRows, nCols);
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" nCols; ++j)
result[i][j] = m[i][j];
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 = 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] += matA[i][k] * matB[k][j];
return result;
}
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;
}
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);
}
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;
}
} // MatInverseQR
// ======================================================
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 double[][] MatLoad(string fn, int[] usecols,
char sep, string comment)
{
// self-contained version
// first, count number of non-comment lines
int nRows = 0;
string line = "";
FileStream ifs = new FileStream(fn, FileMode.Open);
StreamReader sr = new StreamReader(ifs);
while ((line = sr.ReadLine()) != null)
if (line.StartsWith(comment) == false)
++nRows;
sr.Close(); ifs.Close();
int nCols = usecols.Length;
double[][] result = new double[nRows][];
for (int r = 0; r "lt" nRows; ++r)
result[r] = new double[nCols];
line = "";
string[] tokens = null;
ifs = new FileStream(fn, FileMode.Open);
sr = new StreamReader(ifs);
int i = 0;
while ((line = sr.ReadLine()) != null)
{
if (line.StartsWith(comment) == true)
continue;
tokens = line.Split(sep);
for (int j = 0; j "lt" nCols; ++j)
{
int k = usecols[j]; // into tokens
result[i][j] = double.Parse(tokens[k]);
}
++i;
}
sr.Close(); ifs.Close();
return result;
}
} // Program
} // ns
.NET Test Automation Recipes
Software Testing
SciPy Programming Succinctly
Keras Succinctly
R Programming
2024 Visual Studio Live Conference
2024 Predictive Analytics World
2024 DevIntersection Conference
2023 Fall MLADS Conference
2022 Money 20/20 Conference
2022 DEFCON Conference
2022 G2E Conference
2023 ICGRT Conference
2024 CEC eSports Conference
2024 ISC West Conference
James D. McCaffrey 
You must be logged in to post a comment.