I cleaned up an old implementation of a C# function to compute the Moore-Penrose pseudo-inverse of a matrix, via singular value decomposition (SVD). A regular inverse works only for square matrices with the same number of rows and columns. But in machine learning there are many problems where you want to compute an inverse of a non-square matrix. One example is computing a closed-form solution for the weights of a linear regression model.
Computing the Moore-Penrose pseudo-inverse is simultaneously easy and extremely difficult — you apply singular value decomposition (SVD) to the source matrix, take the three resulting matrices, compute their inverses (which is very easy) and then multiply the three inverses together. Easy. Unfortunately, computing the singular value decomposition of a matrix is one of the most difficult problems in numerical programming.
Before I go any further, let me point out that after a lot of experimentation, I have concluded that computing the pseudo-inverse of a matrix using from-scratch SVD-Jacobi is not as numerically stable as computing the pseudo-inverse using from scratch QR-Householder. Therefore, the code in this blog post is mostly for “historical” (to me) interest. See jamesmccaffreyblog.com/2025/11/21/matrix-pseudo-inverse-using-qr-householder-with-csharp/ for an example of computing matrix pseudo-inverse from scratch using QR-Householder.
To compute SVD, I used the Jacobi algorithm. I based my C# version on a GNU scientific library C language implementation. It took many hours.
One practical problem I faced was the need to implement a whole slew of helper functions such as MatMake(), MatTranpose(), MatProduct(), and so on. To manage the complexity of all the functions, I initially put the pseudo-inverse function, and the main worker SVD function, and all the helper functions, in a PseudoInverse class. But that design just didn’t feel right. So I refactored again to a standalone MatDecompSVD() function without any of the many helper functions, and a standalone MatPinv() function.
The Pinv() function is deceptively simple-looking:
public static double[][] MatPinv(double[][] M)
{
// Moore-Penrose pseudo-inverse from SVD (Jacobi)
double[][] U; double[] s; double[][] Vh;
MatDecompSVD(M, out U, out s, out Vh);
double[][] Sinv = MatMake(s.Length, s.Length);
for (int i = 0; i < Sinv.Length; ++i)
Sinv[i][i] = 1.0 / (s[i] + 1.0e-12); // avoid 0
// pinv = V * Sinv * Uh
double[][] V = MatTranspose(Vh);
double[][] Uh = MatTranspose(U);
double[][] tmp = MatProduct(Sinv, Uh);
double[][] result = MatProduct(V, tmp);
return result;
// --------------------------------------------------------
// helpers MatMake(), MatTranspose(), MatProduct() here
// --------------------------------------------------------
}
The call to MatDecompSVD() returns a matrix U, a vector s (the singular values), and a vector Vh. The inverse of U, Uh, is just the transpose of U. The matrix inverse of vector s is a square matrix with the reciprocals of the values of s on the diagonal. The inverse of Vh, V, is just the transpose of Vh. The pseudo-inverse is V * Sinv * Uh.
The output of my demo is:
Begin matrix pseudo-inverse using C# A: 1.0000 2.0000 3.0000 5.0000 0.0000 2.0000 8.0000 5.0000 4.0000 1.0000 0.0000 9.0000 Apinv: -0.0784 0.1706 0.0314 -0.0257 0.1568 -0.2390 0.1373 -0.0602 0.0287 -0.0091 -0.0115 0.1087 B: -1.0000 2.0000 3.0000 9.0000 5.0000 0.0000 -2.0000 4.0000 8.0000 -5.0000 4.0000 7.0000 Bpinv: -0.0572 0.0945 0.0363 0.0489 0.0838 -0.0832 0.0269 -0.1680 0.0815 0.0849 0.0479 -0.0046 End demo
I validated my implementation using the Python NumPy library linalg.pinv() function:
C:\VSM\PseudoInverse: python pseudoinv_numpy.py A: [[1 2 3] [5 0 2] [8 5 4] [1 0 9]] Apinv: [[-0.0784 0.1706 0.0314 -0.0257] [ 0.1568 -0.239 0.1373 -0.0602] [ 0.0287 -0.0091 -0.0115 0.1087]] B: [[-1 2 3 9] [ 5 0 -2 4] [ 8 -5 4 7]] Bpinv: [[-0.0572 0.0945 0.0363] [ 0.0489 0.0838 -0.0832] [ 0.0269 -0.168 0.0815] [ 0.0849 0.0479 -0.0046]] C:\VSM\PseudoInverse:
Very interesting, lots of fun, very satisfying.
Algorithms like the Jacobi Moore-Penrose pseudo-inverse are old but still effective and useful. I have always been absolutely fascinated by old electro-mechanical baseball games. Hundreds of variations from many (I'm guessing about two dozen) manufacturers were made from the late 1920s through the late 1970s. The machines were marvels of design and are still fun to play even today.
Demo program. Replace "lt" (less than), "gt", "lte", "gte", "and" with Boolean operator symbols. My blog editor often chokes on symbols.
using System;
using System.IO;
using System.Collections.Generic;
namespace PseudoInverse
{
internal class PseudoInverseProgram
{
static void Main(string[] args)
{
Console.WriteLine("\nBegin matrix pseudo-inverse" +
" using C# ");
// 1. m (rows) greater than n (cols)
double[][] A = new double[4][];
A[0] = new double[] { 1, 2, 3 };
A[1] = new double[] { 5, 0, 2 };
A[2] = new double[] { 8, 5, 4 };
A[3] = new double[] { 1, 0, 9 };
Console.WriteLine("\nA: ");
for (int i = 0; i "lt" A.Length; ++i)
VecShow(A[i], 4, 8);
double[][] Apinv = MatPinv(A);
Console.WriteLine("\nApinv: ");
for (int i = 0; i "lt" Apinv.Length; ++i)
VecShow(Apinv[i], 4, 8);
// 2. m less than n
double[][] B = new double[3][];
B[0] = new double[] { -1, 2, 3, 9 };
B[1] = new double[] { 5, 0, -2, 4 };
B[2] = new double[] { 8, -5, 4, 7 };
Console.WriteLine("\nB: ");
for (int i = 0; i "lt" B.Length; ++i)
VecShow(B[i], 4, 8);
double[][] Bpinv = MatPinv(B);
Console.WriteLine("\nBpinv: ");
for (int i = 0; i "lt" Bpinv.Length; ++i)
VecShow(Bpinv[i], 4, 8);
Console.WriteLine("\nEnd demo ");
Console.ReadLine();
} // Main
// ------------------------------------------------------
// helpers for Main(): MatLoad(), VecShow()
static double[][] MatLoad(string fn, int[] usecols,
char sep, string comment) // not used this demo
{
List"lt"double[]"gt" result = new List"lt"double[]"gt"();
string line = "";
FileStream ifs = new FileStream(fn, FileMode.Open);
StreamReader sr = new StreamReader(ifs);
while ((line = sr.ReadLine()) != null)
{
if (line.StartsWith(comment) == true)
continue;
string[] tokens = line.Split(sep);
List"lt"double"gt" lst = new List"lt"double"gt"();
for (int j = 0; j "lt" usecols.Length; ++j)
lst.Add(double.Parse(tokens[usecols[j]]));
double[] row = lst.ToArray();
result.Add(row);
}
sr.Close(); ifs.Close();
return result.ToArray();
}
// ------------------------------------------------------
static void VecShow(double[] vec, int dec, int wid)
{
double small = 1.0 / Math.Pow(10, dec);
for (int i = 0; i "lt" vec.Length; ++i)
{
double v = vec[i];
if (Math.Abs(v) "lt" small) v = 0.0;
Console.Write(v.ToString("F" + dec).PadLeft(wid));
}
Console.WriteLine("");
}
// ------------------------------------------------------
public static double[][] MatPinv(double[][] M)
{
// Moore-Penrose pseudo-inverse from SVD (Jacobi)
double[][] U; double[] s; double[][] Vh;
MatDecompSVD(M, out U, out s, out Vh);
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-12); // avoid 0
// pinv = V * Sinv * Uh
double[][] V = MatTranspose(Vh);
double[][] Uh = MatTranspose(U);
double[][] tmp = MatProduct(Sinv, Uh);
double[][] result = MatProduct(V, tmp);
return result;
// ----------------------------------------------------
// helpers MatMake(), MatTranspose(), MatProduct()
// ----------------------------------------------------
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[][] 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;
}
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;
}
} // MatPinv()
// ------------------------------------------------------
public static void MatDecompSVD(double[][] M,
out double[][] U, out double[] s, out double[][] Vh)
{
// Jacobi algorithm
// see github.com/ampl/gsl/blob/master/linalg/svd.c
// "full_matrices=False" version
double DBL_EPSILON = 1.0e-15;
int m = M.Length; // nRows
int n = M[0].Length; // nCols
// A = MatCopy(M);
double[][] A = new double[m][]; // working U
for (int i = 0; i "lt" m; ++i)
A[i] = new double[n];
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" n; ++j)
A[i][j] = M[i][j];
// Q = MatIdentity(n);
double[][] Q = new double[n][]; // working V
for (int i = 0; i "lt" n; ++i)
Q[i] = new double[n];
for (int i = 0; i "lt" n; ++i)
Q[i][i] = 1.0;
double[] t = new double[n]; // working s
// initialize counters
int count = 1;
int sweep = 0;
double tol = 10 * m * DBL_EPSILON; // heuristic
// Always do at least 12 sweeps
int sweepmax = Math.Max(5 * n, 12); // heuristic
// store the column error estimates in St for use
// during orthogonalization
for (int j = 0; j "lt" n; ++j)
{
// double[] cj = MatGetColumn(A, j);
double[] cj = new double[m];
for (int i = 0; i "lt" m; ++i)
cj[i] = A[i][j];
// double sj = VecNorm(cj);
double sj = 0.0;
for (int i = 0; i "lt" cj.Length; ++i)
sj += cj[i] * cj[i];
sj = Math.Sqrt(sj);
t[j] = DBL_EPSILON * sj;
} // j
// orthogonalize A by plane rotations
while (count "gt" 0 "and" sweep "lte" sweepmax)
{
// 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[] cj = new double[m];
for (int i = 0; i "lt" m; ++i)
cj[i] = A[i][j];
// double[] ck = MatGetColumn(A, k);
double[] ck = new double[m];
for (int i = 0; i "lt" m; ++i)
ck[i] = A[i][k];
// double p = 2.0 * VecDot(cj, ck);
double p = 0.0;
for (int i = 0; i "lt" cj.Length; ++i)
p += cj[i] * ck[i];
p *= 2.0;
// double a = VecNorm(cj);
double a = 0.0;
for (int i = 0; i "lt" cj.Length; ++i)
a += cj[i] * cj[i];
a = Math.Sqrt(a);
//double b = VecNorm(ck);
double b = 0.0;
for (int i = 0; i "lt" ck.Length; ++i)
b += ck[i] * ck[i];
b = Math.Sqrt(b);
double q = a * a - b * b;
// double v = Hypot(p, q); // sqrt(x^2 + y^2)
double v = Math.Sqrt((p * p) + (q * q));
// test for columns j,k orthogonal,
// or dominant errors
double abserr_a = t[j];
double abserr_b = t[k];
bool sorted = (a "gte" b);
bool orthog = (Math.Abs(p) "lte"
tol * (a * b));
bool noisya = (a "lt" abserr_a);
bool noisyb = (b "lt" abserr_b);
if (sorted "and" (orthog ||
noisya || noisyb))
{
--count; continue;
}
// calculate rotation angles
if (v == 0 || !sorted)
{
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 (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
t[j] = Math.Abs(cosine) * abserr_a +
Math.Abs(sine) * abserr_b;
t[k] = Math.Abs(sine) * abserr_a +
Math.Abs(cosine) * abserr_b;
// 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 * cosine + Qik * sine;
Q[i][k] = -Qij * sine + Qik * cosine;
} // i
} // k
} // j
++sweep;
} // while
// compute singular values
double prevNorm = -1.0;
for (int j = 0; j "lt" n; ++j)
{
// double[] column = MatGetColumn(A, j);
double[] column = new double[m];
for (int i = 0; i "lt" m; ++i)
column[i] = A[i][j];
// double norm = VecNorm(column);
double norm = 0.0;
for (int i = 0; i "lt" column.Length; ++i)
norm += column[i] * column[i];
norm = Math.Sqrt(norm);
// determine 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 (count "gt" 0)
{
Console.WriteLine("Jacobi iterations did not" +
" converge");
}
// U = A;
U = new double[m][];
for (int i = 0; i "lt" m; ++i)
U[i] = new double[n];
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" n; ++j)
U[i][j] = A[i][j];
// s = t;
s = new double[n];
for (int i = 0; i "lt" n; ++i)
s[i] = t[i];
// Vh = MatTranspose(Q); // Q is nxn
Vh = new double[n][];
for (int i = 0; i "lt" n; ++i)
Vh[i] = new double[n];
for (int i = 0; i "lt" n; ++i)
for (int j = 0; j "lt" n; ++j)
Vh[i][j] = Q[j][i];
// to sync with np.linalg.svd() shapes
// for full_matrices=False (not the default)
// 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);
double[][] newU = new double[U.Length][]; // all rows
for (int i = 0; i "lt" U.Length; ++i)
newU[i] = new double[m]; // just m cols
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" m; ++j) // m cols
newU[i][j] = U[i][j];
U = newU;
// s = VecExtractFirst(s, m);
double[] newS = new double[m]; // just m vals
for (int i = 0; i "lt" m; ++i)
newS[i] = s[i];
s = newS;
// Vh = MatExtractFirstRows(Vh, m);
double[][] newVh = new double[m][]; // just m rows
for (int i = 0; i "lt" m; ++i)
newVh[i] = new double[Vh[0].Length]; // all cols
for (int i = 0; i "lt" m; ++i)
for (int j = 0; j "lt" Vh[0].Length; ++j)
newVh[i][j] = Vh[i][j];
Vh = newVh;
} // m "lt" n
} // MatDecompSVD()
// ------------------------------------------------------
} // class Program
// ========================================================
} // ns
The Python validation program:
# pseudoinv_numpy.py
import numpy as np
np.set_printoptions(precision=4, suppress=True)
A = np.array([
[1, 2, 3],
[5, 0, 2],
[8, 5, 4],
[1, 0, 9]])
print("\nA: ")
print(A)
Apinv = np.linalg.pinv(A)
print("\nApinv: ")
print(Apinv)
B = np.array([
[-1, 2, 3, 9],
[5, 0, -2, 4],
[8, -5, 4, 7]])
print("\nB: ")
print(B)
Bpinv = np.linalg.pinv(B)
print("\nBpinv: ")
print(Bpinv)
.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.