If you have a square matrix A, you can compute the inverse of A so that A * Ainv = I, where * is matrix multiplication and I is the identity matrix. In machine learning scenarios, the A matrix is usually a matrix of training data or a design matrix derived from the training data, and so you need to compute the pseudo-inverse of A.
There are many ways to compute a pseudo-inverse for machine learning scenarios where the number of rows of A are greater than or equal to the number of columns. Three common techniques include singular value decomposition (SVD), QR decomposition, and left pseudo-inverse via normal equations. Each technique has several variations, for example SVD Jacobi, SVD Lanczos, SVD Householder+QR, SVD Golub-Kahan, QR Householder, QR modified Gram-Schmidt, QR Givens, normal equations Cholesky, normal equations LUP, and many others.
Why are there so many techniques to compute a pseudo-inverse? First, the problem is extremely difficult and different techniques have different pros and cons related to stability, performance, and implementation complexity. Second, the problem is so rich in ideas, many mathematicians worked for many years and came up with dozens of techniques.
One afternoon during a work break, I figured I’d refactor my C# implementations of matrix pseudo-inverse using the SVD Jacobi algorithm. It is the most complex of my pseudo-inverse implementations but it is the most stable and the fastest.
To test my refactored version, I generated 5,000 random matrices, computed the pseudo-inverse and then verified that (A * Apinv) * A = A to within 1.0e-8. Each random matrix has between 100 and 1,000 rows, and between 2 and 20 columns, and each value is between -10.0 and +10.0. All 5,000 tests passed. This is by no means thorough testing, but it indicates that basic functionality is correct for most simple machine learning scenarios.
The key calling code is:
for (int i = 0; i "lt" nTrials; ++i) {
Console.WriteLine("trial " + i);
double[][] A = MakeRandomMatrix(minRows, maxRows,
minCols, maxCols, rnd);
double[][] Apinv = SVDJacobi.MatPseudoInv(A);
// check (A * Apinv) * A = A
double[][] AApinv = MatProduct(A, Apinv);
double[][] C = MatProduct(AApinv, A);
if (MatAreEqual(C, A, 1.0e-8) == true)
++nPass;
else
++nFail;
}
The static MatPseudoInv() function is in a container class named SVDJacobi.
Anyway, good fun during a work break. But now it’s Time to get back to work.
The “Lie Detector” game was first produced in 1960. There are 24 character cards, and each is both a witness and a suspect. The game starts by blindly selecting one of the 24 “guily” cards (not shown) and inserting it into the back of the lie detector machine. Then, each player takes a turn taking one of their witness cards, placing on the machine, and inserting an electric probe. The machine needle would swing to “True” or “False”.
Eventually you build up enough information — “overweight”, “wears glasses”, “blonde hair” — to name the guilty person.
I loved this game and was fascinated by the matrix of holes and the combinatorics of the primitive machine. But alas, none of the other kids in the neighborhood (Roger, Beth, Kenny, Tommy, Nancy, Rick, Bobby, Kirk, Lisa, Bill, Sally) liked the game.
Demo program. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols.
using System.IO;
namespace MatrixPseudoInverseSVDJacobi
{
internal class MatrixPseudoInverseSVDJacobi
{
static void Main(string[] args)
{
Console.WriteLine("\nBegin relaxed Moore-Penrose" +
" pseudo-inverse using SVD decomp" +
" (Jacobi) ");
int minRows = 100; int maxRows = 1000;
int minCols = 2; int maxCols = 20;
Random rnd = new Random(0);
int nTrials = 5000;
int nPass = 0; int nFail = 0;
for (int i = 0; i "lt" nTrials; ++i)
{
Console.WriteLine("\n==========");
Console.WriteLine("trial " + i);
double[][] A =
MakeRandomMatrix(minRows, maxRows,
minCols, maxCols, rnd);
double[][] Apinv = SVDJacobi.MatPseudoInv(A);
// check (A * Apinv) * A = A
double[][] AApinv = MatProduct(A, Apinv);
double[][] C = MatProduct(AApinv, A);
//Console.WriteLine("\nA = ");
//MatShow(A, 4, 9);
//Console.WriteLine("\nApinv = ");
//MatShow(Apinv, 4, 9);
//Console.WriteLine("\n(A * Apinv) * A = ");
//MatShow(C, 4, 9);
if (MatAreEqual(C, A, 1.0e-8) == true)
{
Console.WriteLine("pass");
++nPass;
}
else
{
Console.WriteLine("FAIL");
Console.WriteLine("nRows = " + A.Length +
" nCols = " + A[0].Length);
Console.ReadLine();
++nFail;
}
Console.WriteLine("==========");
// Console.ReadLine();
} // nTrials
Console.WriteLine("\nNumber pass = " + nPass);
Console.WriteLine("Number fail = " + nFail);
Console.WriteLine("\nEnd testing ");
Console.ReadLine();
} // Main
// ------------------------------------------------------
// helpers: MakeRandomMatrix(), MatProduct(),
// MatShow(), MatAreEqual()
// ------------------------------------------------------
public static double[][] MakeRandomMatrix(int minRows,
int maxRows, int minCols, int maxCols, Random rnd)
{
int nRows = rnd.Next(minRows, maxRows); // [2,maxN)
int nCols = rnd.Next(minCols, maxCols);
double[][] result = new double[nRows][];
for (int i = 0; i "lt" nRows; ++i)
result[i] = new double[nCols];
double lo = -10.0; double hi = 10.0;
for (int r = 0; r "lt" nRows; ++r)
for (int c = 0; c "lt" nCols; ++c)
result[r][c] = (hi - lo) * rnd.NextDouble() + lo;
return result;
}
// ------------------------------------------------------
public 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 = 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] += A[i][k] * B[k][j];
return result;
}
// ------------------------------------------------------
public static void MatShow(double[][] m, int dec, int wid)
{
int nRows = m.Length; int nCols = m[0].Length;
double small = 1.0 / Math.Pow(10, dec);
for (int i = 0; i "lt" nRows; ++i)
{
for (int j = 0; j "lt" nCols; ++j)
{
double v = m[i][j];
if (Math.Abs(v) "lt" small) v = 0.0;
Console.Write(v.ToString("F" + dec).
PadLeft(wid));
}
Console.WriteLine("");
}
}
// ------------------------------------------------------
public static bool MatAreEqual(double[][] A, double[][] B,
double epsilon)
{
int nRows = A.Length; int nCols = A[0].Length;
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" nCols; ++j)
if (Math.Abs( A[i][j] - B[i][j] ) "gt" epsilon)
return false;
return true;
}
// ------------------------------------------------------
} // Program
// ========================================================
public class SVDJacobi
{
// ------------------------------------------------------
// "Compact Numerical Methods for Computers" (1979).
// Nash, J.C.
// "Jacobi's Method is More Accurate than QR" (1989),
// Demmel, J. and Veselic, K.
// C language code (unknown location) by Kosowsky, A.
// (Rutgers University)
// C language code at
// github.com/ampl/gsl/blob/master/linalg by
// Jungman, G. (Los Alamos National Laboratory)
// "Algorithm for Computing the Singular Value
// Decomposition on a Vector Computer" (1989),
// De Rijk, P.P.M.
// NumPy documentation at
// numpy.org/doc/2.1/reference/generated/
// numpy.linalg.svd.html
// ------------------------------------------------------
public static double[][] MatPseudoInv(double[][] M)
{
double[][] U; double[][] Vh; double[] s;
SVD_Jacobi(M, out U, out Vh, out s);
// pinv = V * Sinv * Uh
double[][] Sinv = MatMake(s.Length, s.Length);
for (int i = 0; i "lt" Sinv.Length; ++i)
Sinv[i][i] = 1.0 / (s[i] + 1.0e-8); // no div 0
double[][] V = MatTranspose(Vh);
double[][] Uh = MatTranspose(U);
double[][] tmp = MatProduct(Sinv, Uh);
double[][] result = MatProduct(V, tmp);
return result;
}
// ------------------------------------------------------
public static void SVD_Jacobi(double[][] M,
out double[][] U, out double[][] Vh, out double[] s)
{
// see references above
double epsilon = 1.0e-15;
double[][] A = MatCopyOf(M); // working U
int m = A.Length; int n = A[0].Length;
double[][] Q = MatIdentity(n); // working V
double[] sv = new double[n]; // working s
// initialize counters
int count = 1;
int pass = 0;
double tolerance = 10 * m * epsilon; // heuristic
// always do at least 12 sweeps
int passMax = Math.Max(5 * n, 12); // heuristic
// store the column error estimates for use
// during orthogonalization
for (int j = 0; j "lt" n; ++j)
{
double[] cj = MatGetColumn(A, j);
double sj = VecNorm(cj);
sv[j] = epsilon * sj;
}
// orthogonalize A using plane rotations
while (count "gt" 0 && pass "lte" passMax)
{
// initialize rotation counter
count = n * (n - 1) / 2;
for (int j = 0; j "lt" n - 1; ++j)
{
for (int k = j + 1; k "lt" n; ++k)
{
double cosine, sine;
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);
// test for columns j,k orthogonal,
// or dominant errors
double abserr_a = sv[j];
double abserr_b = sv[k];
bool sorted = (a "gte" b);
bool orthog = (Math.Abs(p) "lte"
tolerance * (a * b));
bool bada = (a "lt" abserr_a);
bool badb = (b "lt" abserr_b);
if (sorted == true && (orthog == true ||
bada == true || badb == true))
{
--count;
continue;
}
// calculate rotation angles
if (v == 0 || sorted == false)
{
cosine = 0.0; sine = 1.0;
}
else
{
cosine = Math.Sqrt((v + q) / (2.0 * v));
sine = p / (2.0 * v * cosine);
}
// apply rotation to A (working U)
for (int i = 0; i "lt" m; ++i)
{
double Aik = A[i][k];
double Aij = A[i][j];
A[i][j] = Aij * cosine + Aik * sine;
A[i][k] = -Aij * sine + Aik * cosine;
}
// update singular values
sv[j] = Math.Abs(cosine) * abserr_a +
Math.Abs(sine) * abserr_b;
sv[k] = Math.Abs(sine) * abserr_a +
Math.Abs(cosine) * abserr_b;
// apply rotation to Q (working V)
for (int i = 0; i "lt" n; ++i)
{
double Qij = Q[i][j];
double Qik = Q[i][k];
Q[i][j] = Qij * cosine + Qik * sine;
Q[i][k] = -Qij * sine + Qik * cosine;
} // 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);
// determine if singular value is zero
if (norm == 0.0 || prevNorm == 0.0
|| (j "gt" 0 &&
norm "lte" tolerance * prevNorm))
{
sv[j] = 0.0;
for (int i = 0; i "lt" m; ++i)
A[i][j] = 0.0;
prevNorm = 0.0;
}
else
{
sv[j] = norm;
for (int i = 0; i "lt" m; ++i)
A[i][j] = A[i][j] * 1.0 / norm;
prevNorm = norm;
}
}
if (count "gt" 0)
{
Console.WriteLine("Jacobi iterations did not" +
" converge");
}
U = A;
Vh = MatTranspose(Q);
s = sv;
// to sync with np.linalg.svd() full_matrices=False
// shapes and allow U*S*Vh:
// if m "lt" n, extract 1st m columns of U
// extract 1st m values of s
// extract 1st m rows of Vh
if (m "lt" n)
{
U = MatExtractFirstColumns(U, m);
s = VecExtractFirst(s, m);
Vh = MatExtractFirstRows(Vh, m);
}
} // SVD_Jacobi()
// === helper functions =================================
//
// MatMake, MatCopy, MatIdentity, MatGetColumn,
// MatExtractFirstColumns, MatExtractFirstRows,
// MatTranspose, MatProduct, VecNorm, VecDot,
// Hypot, VecExtractFirst
//
// ======================================================
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[][] MatCopyOf(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[][] 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[] MatGetColumn(double[][] M, int j)
{
int nRows = M.Length;
double[] result = new double[nRows];
for (int i = 0; i "lt" nRows; ++i)
result[i] = M[i][j];
return result;
}
// ------------------------------------------------------
private static double[][]
MatExtractFirstColumns(double[][] M, int n)
{
int nRows = M.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] = M[i][j];
return result;
}
// ------------------------------------------------------
private static double[][]
MatExtractFirstRows(double[][] M, int n)
{
int nCols = M[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] = M[i][j];
return result;
}
// ------------------------------------------------------
private static double[][] MatTranspose(double[][] M)
{
int nRows = M.Length;
int nCols = M[0].Length;
double[][] result = MatMake(nCols, nRows);
for (int i = 0; i "lt" nRows; ++i)
for (int j = 0; j "lt" nCols; ++j)
result[j][i] = 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)
for (int j = 0; j "lt" bCols; ++j)
for (int k = 0; k "lt" aCols; ++k)
result[i][j] += A[i][k] * B[k][j];
return result;
}
// ------------------------------------------------------
private static double VecNorm(double[] v)
{
double sum = 0.0;
int n = v.Length;
for (int i = 0; i "lt" n; ++i)
sum += v[i] * v[i];
return Math.Sqrt(sum);
}
// ------------------------------------------------------
private 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;
}
// ------------------------------------------------------
private static double Hypot(double x, double y)
{
// fancy sqrt(x^2 + y^2) std technique
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);
}
// ------------------------------------------------------
private static double[] VecExtractFirst(double[] v,
int n)
{
double[] result = new double[n];
for (int i = 0; i "lt" n; ++i)
result[i] = v[i];
return result;
}
// ------------------------------------------------------
} // class SVDJacobi
// ========================================================
} // 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.