One morning, before work, I decided to tidy up my C# implementation of matrix inverse using singular value decomposition (SVD). My main change was to refactor the MatInverseSVD() function so that all the helper functions (11 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. Some common variations for SVD decomposition are symmetric Jacobi, asymmetric Jacobi, Householder reduction, Golub–Reinsch, Demmel–Kahan, and bidirectional bisection.
My demo uses the symmetric Jacobi version of SVD. It assumes the source matrix is square, which is OK because matrix inverse applies only to square matrices. (There is a “pseudo-inverse” for non-square matrices, but that’s a different topic).
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 values indicate a very small negative number, truncated for display, not a special negative-zero value of some sort.
Good fun!
I’m a big fan of comic books from the Silver Age (1956-1970). A common item was a ray gun that could presumably decompose matter.
Left: “The Flash” #105 in January 1959. This was the first apearance of a new Flash character who succeeded the old 1940s Flash.
Center: “Mystery in Space” #53 in August 1959. This was the first appearance of Adam Strange, an Earthman without special powers who was a hero on the planet Rann.
Right: “Adventure Comics” #325 in October 1964. This was the first time the evil Lex Luthor tangled with the Legion of Super-Heroes (which included Superboy and Supergirl).
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 MatrixInverseSVD
{
internal class MatrixInverseSVDProgram
{
static void Main(string[] args)
{
Console.WriteLine("\nMatrix inverse using SVD ");
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 = MatInverseSVD(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[][] MatInverseSVD(double[][] M)
{
// SVD Jacobi algorithm
// A = U * S * Vh
// inv(A) = tr(Vh) * inv(S) * tr(U)
double[][] U;
double[][] Vh;
double[] s;
MatDecomposeSVD(M, out U, out Vh, out s);
// TODO: check if determinant is zero (no inverse)
// double absDet = 1.0;
// for (int i = 0; i "lt" M.Length; ++i)
// absDet *= s[i];
// Console.WriteLine(absDet); Console.ReadLine();
// convert s vector to Sinv matrix
double[][] Sinv = MatInvFromVec(s);
double[][] V = MatTranspose(Vh);
double[][] Utrans = MatTranspose(U);
double[][] result = MatProduct(V,
MatProduct(Sinv, Utrans));
return result;
// ----------------------------------------------------
// 11 helpers: MatDecomposeSVD, MatMake, MatCopy,
// MatIdentity, MatGetColumn, MatTranspose,
// MatInvFromVec, MatProduct, VecNorm, VecDot, Hypot
// ----------------------------------------------------
static void MatDecomposeSVD(double[][] mat,
out double[][] U, out double[][] Vh,
out double[] s)
{
// assumes source matrix is square
double EPSILON = 1.0e-15;
double[][] A = MatCopy(mat); // working U
int m = A.Length;
int n = A[0].Length; // check m == n
double[][] Q = MatIdentity(n); // working V
double[] t = new double[n]; // working s
int ct = 1; // rotation counter
int pass = 0;
double tol = 10 * n * EPSILON; // heuristic
int passMax = 5 * n;
if (passMax "lt" 15) passMax = 15; // heuristic
// save the column error estimates
for (int j = 0; j "lt" n; ++j)
{
double[] cj = MatGetColumn(A, j);
double sj = VecNorm(cj);
t[j] = EPSILON * sj;
}
while (ct "gt" 0 "and" pass "lte" passMax)
{
ct = n * (n - 1) / 2; // rotation counter
for (int j = 0; j "lt" n - 1; ++j)
{
for (int k = j + 1; k "lt" n; ++k)
{
double sin; double cos;
double[] cj = MatGetColumn(A, j);
double[] ck = MatGetColumn(A, k);
double p = 2.0 * VecDot(cj, ck);
double a = VecNorm(cj);
double b = VecNorm(ck);
double q = a * a - b * b;
double v = Hypot(p, q);
double errorA = t[j];
double errorB = t[k];
bool sorted = false;
if (a "gte" b) sorted = true;
bool orthog = false;
if (Math.Abs(p) "lte"
tol * (a * b)) orthog = true;
bool badA = false;
if (a "lt" errorA) badA = true;
bool badB = false;
if (b "lt" errorB) badB = true;
if (sorted == true "and" (orthog == true ||
badA == true || badB == true))
{
--ct;
continue;
}
// compute rotation angles
if (v == 0 || sorted == false)
{
cos = 0.0;
sin = 1.0;
}
else
{
cos = Math.Sqrt((v + q) / (2.0 * v));
sin = p / (2.0 * v * cos);
}
// apply rotation to A (U)
for (int i = 0; i "lt" m; ++i)
{
double Aik = A[i][k];
double Aij = A[i][j];
A[i][j] = Aij * cos + Aik * sin;
A[i][k] = -Aij * sin + Aik * cos;
}
// update singular values
t[j] = Math.Abs(cos) * errorA +
Math.Abs(sin) * errorB;
t[k] = Math.Abs(sin) * errorA +
Math.Abs(cos) * errorB;
// apply rotation to Q (V)
for (int i = 0; i "lt" n; ++i)
{
double Qij = Q[i][j];
double Qik = Q[i][k];
Q[i][j] = Qij * cos + Qik * sin;
Q[i][k] = -Qij * sin + Qik * cos;
} // i
} // k
} // j
++pass;
} // while
// compute singular values
double prevNorm = -1.0;
for (int j = 0; j "lt" n; ++j)
{
double[] column = MatGetColumn(A, j);
double norm = VecNorm(column);
// check if singular value is zero
if (norm == 0.0 || prevNorm == 0.0
|| (j "gt" 0 "and" norm "lte" tol * prevNorm))
{
t[j] = 0.0;
for (int i = 0; i "lt" m; ++i)
A[i][j] = 0.0;
prevNorm = 0.0;
}
else
{
t[j] = norm;
for (int i = 0; i "lt" m; ++i)
A[i][j] = A[i][j] * 1.0 / norm;
prevNorm = norm;
}
}
if (ct "gt" 0)
{
Console.WriteLine("Jacobi iterations did not" +
" converge");
}
U = A;
Vh = MatTranspose(Q);
s = t;
// to sync with default np.linalg.svd() shapes:
// if m "lt" n, extract 1st m columns of U
// extract 1st m values of s
// extract 1st m rows of Vh
// not applicable for matrix inverse.
// if (m "lt" n)
// {
// U = MatExtractFirstColumns(U, m);
// s = VecExtractFirst(s, m);
// Vh = MatExtractFirstRows(Vh, m);
// }
} // MatDecomposeSVD
static double[][] MatMake(int r, int c)
{
double[][] result = new double[r][];
for (int i = 0; i "lt" r; ++i)
result[i] = new double[c];
return result;
}
static double[][] MatCopy(double[][] m)
{
int r = m.Length; int c = m[0].Length;
double[][] result = MatMake(r, c);
for (int i = 0; i "lt" r; ++i)
for (int j = 0; j "lt" c; ++j)
result[i][j] = m[i][j];
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[] MatGetColumn(double[][] m, int j)
{
int rows = m.Length;
double[] result = new double[rows];
for (int i = 0; i "lt" rows; ++i)
result[i] = m[i][j];
return result;
}
static double[][] MatTranspose(double[][] m)
{
int r = m.Length;
int c = m[0].Length;
double[][] result = MatMake(c, r);
for (int i = 0; i "lt" r; ++i)
for (int j = 0; j "lt" c; ++j)
result[j][i] = m[i][j];
return result;
}
static double[][] MatInvFromVec(double[] s)
{
// Sinv from s
int n = s.Length;
double[][] result = MatMake(n, n);
for (int i = 0; i "lt" n; ++i)
result[i][i] = 1.0 / s[i];
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)
for (int j = 0; j "lt" bCols; ++j)
for (int k = 0; k "lt" aCols; ++k)
result[i][j] += matA[i][k] * matB[k][j];
return result;
}
static double VecNorm(double[] vec)
{
double sum = 0.0;
int n = vec.Length;
for (int i = 0; i "lt" n; ++i)
sum += vec[i] * vec[i];
return Math.Sqrt(sum);
}
static double VecDot(double[] v1, double[] v2)
{
int n = v1.Length;
double sum = 0.0;
for (int i = 0; i "lt" n; ++i)
sum += v1[i] * v2[i];
return sum;
}
static double Hypot(double x, double y)
{
// fancy sqrt(x^2 + y^2)
double xabs = Math.Abs(x);
double yabs = Math.Abs(y);
double min, max;
if (xabs "lt" yabs)
{
min = xabs; max = yabs;
}
else
{
min = yabs; max = xabs;
}
if (min == 0)
return max;
double u = min / max;
return max * Math.Sqrt(1 + u * u);
}
// The following three helpers are used by
// SVD decomposition when m != n and so aren't
// needed for matrix inverse
//static double[][] MatExtractFirstColumns(double[][]
// src, int n)
//{
// int nRows = src.Length;
// // int nCols = src[0].Length;
// double[][] result = MatMake(nRows, n);
// for (int i = 0; i "lt" nRows; ++i)
// for (int j = 0; j "lt" n; ++j)
// result[i][j] = src[i][j];
// return result;
//}
//static double[][] MatExtractFirstRows(double[][]
// src, int n)
//{
// // int nRows = src.Length;
// int nCols = src[0].Length;
// double[][] result = MatMake(n, nCols);
// for (int i = 0; i "lt" n; ++i)
// for (int j = 0; j "lt" nCols; ++j)
// result[i][j] = src[i][j];
// return result;
//}
//static double[] VecExtractFirst(double[] vec, int n)
//{
// double[] result = new double[n];
// for (int i = 0; i "lt" n; ++i)
// result[i] = vec[i];
// return result;
//}
} // MatInverseSVD
// ======================================================
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.