Note: I wasn’t satisfied with the SVD code implementation in this blog post. I did a much better implementation in “Singular Value Decomposition (SVD) from Scratch with C# Using the Jacobi Algorithm”. See https://jamesmccaffreyblog.com/2023/11/24/singular-value-decomposition-svd-from-scratch-with-csharp-using-the-jacobi-algorithm/.
The idea of singular value decomposition (SVD) is to take a source matrix A and decompose it into a U matrix, an s vector, and a Vh matrix (the h in Vh stands for Hermitian transpose). For example, if a source matrix A is:
[1 2 3] [5 0 2] [8 5 4] [1 0 9]
Then the SVD produces:
U = [-0.2581 -0.12 -0.4902] [-0.3515 0.2318 0.811 ] [-0.7265 0.5233 -0.2994] [-0.531 -0.8112 0.1111] s = [13.0781 7.1542 2.7892] Vh = [-0.6392 -0.3172 -0.7006] [ 0.6171 0.3322 -0.7134] [ 0.459 -0.8883 -0.0166]
Now if you convert the s (“singular values”) vector to a diagonal matrix S:
diag(S) = [13.0781 0 0 ] [ 0 7.1542 0 ] [ 0 0 2.7892]
and then perform matrix multiplication:
U * diag(S) * Vh
the result is the original source matrix A.
The left shell shows my from-scratch C# version. The shell on the right is the NumPy library version. Results are identical.
There are several algorithms to compute SVD. One weekend, I decided to implement SVD from scratch using the C# language. What an adventure it turned out to be. I eventually determined that I needed to leverage existing resources. The algorithm for SVD that I kludged together from several sources is extremely complex — and extremely ugly.
I used the numpy.linalg.svd() library function to guide my interface design. After a lot of work — and I mean many, many hours over the course of five days — I had a preliminary demo up and running. Calling the function looks like:
double[][] A = new double[4][]; // 4x3
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("Source A = ");
MatShow(A, 1, 5); // program-defined helper
double[][] U;
double[] S;
double[][] Vh;
MatDecompSVD(A, out U, out S, out Vh);
Console.WriteLine("U = ");
MatShow(U, 4, 9);
Console.WriteLine("S = ");
VecShow(S, 4, 9);
Console.WriteLine("Vh = ");
MatShow(Vh, 4, 9);
double[][] USVh =
MatProduct(U, MatProduct(MatDiag(S), Vh));
Console.WriteLine("U * diag(S) * Vh = ");
MatShow(USVh, 1, 5);
My from-scratch SVD function returns Vh, the transpose of V, because that’s what the numpy.linalg.svd() function does. Some library implementations, notably MatLab, return V rather than the transpose of V. Similarly, some library SVD functions return S as a diagonal matrix, rather than a vector of the diagonal values.
Using SVD on a matrix has a lot of important uses in machine learning. You can use the U, S, Vh results to compute the inverse of a source matrix A if A is square. And you can use the SVD results to compute a pseudo-inverse of A if A is not square. Anyway, creating an SVD function from scratch using C# was a lot of fun and very interesting.
I am not at all satisfied with this implementation, but at least it works. When I get a chance, I intend to try an implementation that uses what is called the Jacobi algorithm. Note: I did so. See note below.
Matrix SVD decomposition has many practical uses in machine learning algorithms. A matrix decomposition is either correct or incorrect (subject to numerical rounding and precision issues). Photo decomposition in the form of a dispersion effect is subjective in the sense that an attractive image is difficult to define. Here are three examples.
Note: This version of SVD has been replaced with a better one. See https://jamesmccaffreyblog.com/2023/11/24/singular-value-decomposition-svd-from-scratch-with-csharp-using-the-jacobi-algorithm/.
Demo code. Very long. Very complex. Very ugly. Not thoroughly tested. Replace “lt”, “gt”, “lte”, “gte”, “and” with Boolean operator symbols. The code for the NumPy library version is at the bottom of this post.
using System;
using System.IO;
// kludged together from many sources, but especially:
// GNU scientific library
// Accord.Net library
// Lutz Roeder Mapack library
// numpy.linalg library
// Numerical Recipes in C library
// Java Jama library
namespace SingularValueDecomposition
{
internal class SVDProgram
{
static void Main(string[] args)
{
Console.WriteLine("\nSingular value decomp with C ");
double[][] A = new double[4][]; // 4x3
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("\nSource A = ");
MatShow(A, 1, 5);
double[][] U;
double[] S;
double[][] Vh;
MatDecompSVD(A, out U, out S, out Vh);
Console.WriteLine("\nU = ");
MatShow(U, 4, 9);
Console.WriteLine("\nS = ");
VecShow(S, 4, 9);
Console.WriteLine("\nVh = ");
MatShow(Vh, 4, 9);
double[][] USVh =
MatProduct(U, MatProduct(MatDiag(S), Vh));
Console.WriteLine("\nU * diag(S) * Vh = ");
MatShow(USVh, 1, 5);
Console.WriteLine("\nEnd C 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];
if (Math.Abs(v) "lt" 1.0e-8) v = 0.0; // hack
Console.Write(v.ToString("F" + dec).
PadLeft(wid));
}
Console.WriteLine("");
}
}
// ------------------------------------------------------
static void VecShow(double[] vec,
int dec, int wid)
{
for (int i = 0; i "lt" vec.Length; ++i)
{
double x = vec[i];
if (Math.Abs(x) "lt" 1.0e-8) x = 0.0;
Console.Write(x.ToString("F" + dec).
PadLeft(wid));
}
Console.WriteLine("");
}
// ------------------------------------------------------
static double[][] MatDiag(double[] vec)
{
int n = vec.Length;
double[][] result = MatCreate(n, n);
for (int i = 0; i "lt" n; ++i)
result[i][i] = vec[i];
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 = MatCreate(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;
}
// ------------------------------------------------------
static double[][] MatCreate(int rows, int cols)
{
double[][] result = new double[rows][];
for (int i = 0; i "lt" rows; ++i)
result[i] = new double[cols];
return result;
}
// ------------------------------------------------------
static void MatDecompSVD(double[][] A,
out double[][] U, out double[] S, out double[][] Vh)
{
double[][] a = LocalMatCopy(A); // local function below
int m = a.Length;
int n = a[0].Length;
int mn = Math.Min(m, n);
double[] s = new double[Math.Min(m + 1, n)];
double[][] u = LocalMatCreate(m, mn);
double[][] v = LocalMatCreate(n, n);
double[] z = new double[n];
double[] scratch = new double[m];
// transform A to bidiagonal form
// store the diagonal elements in s
// store the super-diagonal elements in z
int nct = Math.Min(m - 1, n);
int nrt = Math.Max(0, Math.Min(n - 2, m));
for (int k = 0; k "lt" Math.Max(nct, nrt); ++k)
{
if (k "lt" nct)
{
// compute transformation for the
// kth column and place it in s[k].
// compute two-norm of kth column without
// under/overflow.
s[k] = 0;
for (int i = k; i "lt" m; ++i)
s[k] = LocalHypot(s[k], a[i][k]);
if (s[k] != 0.0)
{
if (a[k][k] "lt" 0.0)
s[k] = -s[k];
for (int i = k; i "lt" m; ++i)
a[i][k] /= s[k];
a[k][k] += 1.0;
}
s[k] = -s[k];
}
for (int j = k + 1; j "lt" n; ++j)
{
if ((k "lt" nct) "and" (s[k] != 0.0)) // fixed
{
double t = 0;
for (int i = k; i "lt" m; ++i)
t += a[i][k] * a[i][j];
t = -1 * t / a[k][k];
for (int i = k; i "lt" m; ++i)
a[i][j] += t * a[i][k];
}
// Place the kth row of A into z for
z[j] = a[k][j];
}
if (k "lt" nct)
{
// Place the transformation in U
for (int i = k; i "lt" m; ++i)
u[i][k] = a[i][k];
}
if (k "lt" nrt)
{
// compute the kth row transformation
// and place the kth super-diagonal in z[k].
// compute two-norm without under/overflow.
z[k] = 0;
for (int i = k + 1; i "lt" n; ++i)
{
z[k] = LocalHypot(z[k], z[i]);
}
if (z[k] != 0.0)
{
if (z[k + 1] "lt" 0.0)
z[k] = -1 * z[k];
for (int i = k + 1; i "lt" n; ++i)
z[i] /= z[k];
z[k + 1] += 1.0;
}
z[k] = -1 * z[k];
if ((k + 1 "lt" m) "and" (z[k] != 0.0)) // bug fixed
{
// Apply the transformation.
for (int i = k + 1; i "lt" m; ++i)
scratch[i] = 0.0;
for (int j = k + 1; j "lt" n; ++j)
for (int i = k + 1; i "lt" m; ++i)
scratch[i] += z[j] * a[i][j];
for (int j = k + 1; j "lt" n; ++j)
{
double t = -1 * z[j] / z[k + 1];
for (int i = k + 1; i "lt" m; ++i)
a[i][j] += t * scratch[i];
}
}
// place the transformation into V
for (int i = k + 1; i "lt" n; ++i)
v[i][k] = z[i];
} // if k "lt" nrt
} // for k
// bidiagonal matrix of p
int p = Math.Min(n, m + 1);
if (nct "lt" n) s[nct] = a[nct][nct];
if (m "lt" p) s[p - 1] = 0.0;
if (nrt + 1 "lt" p) z[nrt] = a[nrt][p - 1];
z[p - 1] = 0.0;
for (int j = nct; j "lt" mn; ++j)
{
for (int i = 0; i "lt" m; ++i)
u[i][j] = 0.0;
u[j][j] = 1.0;
}
for (int k = nct - 1; k "gte" 0; --k)
{
if (s[k] != 0.0)
{
for (int j = k + 1; j "lt" mn; ++j)
{
double t = 0;
for (int i = k; i "lt" m; ++i)
t += u[i][k] * u[i][j];
t = -1 * t / u[k][k];
for (int i = k; i "lt" m; ++i)
u[i][j] += t * u[i][k];
}
for (int i = k; i "lt" m; ++i)
u[i][k] = -1 * u[i][k];
u[k][k] = 1 + u[k][k];
for (int i = 0; i "lt" k - 1; ++i)
u[i][k] = 0.0;
}
else
{
for (int i = 0; i "lt" m; ++i)
u[i][k] = 0.0;
u[k][k] = 1.0;
}
}
for (int k = n - 1; k "gte" 0; --k)
{
if ((k "lt" nrt) "and" (z[k] != 0.0)) // bug fixed
{
for (int j = k + 1; j "lt" n; ++j)
{
double t = 0;
for (int i = k + 1; i "lt" n; ++i)
t += v[i][k] * v[i][j];
t = -1 * t / v[k + 1][k];
for (int i = k + 1; i "lt" n; ++i)
v[i][j] += t * v[i][k];
}
}
for (int i = 0; i "lt" n; ++i)
v[i][k] = 0.0;
v[k][k] = 1.0;
}
// compute singular values
int pp = p - 1;
int iter = 0;
double eps = 1.0e-32;
while (p "gt" 0)
{
if (iter "gte" 100)
Console.WriteLine("Warn: no convergence ");
int k; int branch;
// branch = 1 if s(p) and z[k-1] negligible, k "lt" p
// branch = 2 if s(k) is negligible and k "lt" p
// branch = 3 if z[k-1] is negligible, k "lt" p,
// and s(k), . . . s(p) not negligible (qr step).
// branch = 4 if z(p-1) is negligible (convergence)
for (k = p - 2; k "gte" -1; --k)
{
if (k == -1)
break;
if (Math.Abs(z[k]) "lte" eps * (Math.Abs(s[k]) +
Math.Abs(s[k + 1])))
{
z[k] = 0.0;
break;
}
}
if (k == p - 2)
branch = 4;
else
{
int ks;
for (ks = p - 1; ks "gte" k; --ks)
{
if (ks == k)
break; // out of k loop
double t; double tmp;
if (ks != p)
tmp = Math.Abs(z[ks]);
else
tmp = 0.0;
if (ks != k + 1)
t = tmp + Math.Abs(z[ks - 1]);
else
t = tmp + 0.0;
if (Math.Abs(s[ks]) "lte" eps * t)
{
s[ks] = 0.0;
break; // out of k loop
}
}
if (ks == k)
branch = 3;
else if (ks == p - 1)
branch = 1;
else
{
branch = 2; k = ks;
}
}
++k;
// take appropriate action based on branch
if (branch == 1)
{
double f = z[p - 2];
z[p - 2] = 0.0;
for (int j = p - 2; j "gte" k; --j)
{
double t = LocalHypot(s[j], f);
double cs = s[j] / t;
double sn = f / t;
s[j] = t;
if (j != k)
{
f = -1 * sn * z[j - 1];
z[j - 1] = cs * z[j - 1];
}
for (int i = 0; i "lt" n; ++i)
{
t = cs * v[i][j] + sn * v[i][p - 1];
v[i][p - 1] =
-1 * sn * v[i][j] + cs * v[i][p - 1];
v[i][j] = t;
}
}
} // branch 1
else if (branch == 2)
{
double f = z[k - 1];
z[k - 1] = 0.0;
for (int j = k; j "lt" p; ++j)
{
double t = LocalHypot(s[j], f);
double cs = s[j] / t;
double sn = f / t;
s[j] = t;
f = -sn * z[j];
z[j] = cs * z[j];
for (int i = 0; i "lt" m; ++i)
{
t = cs * u[i][j] + sn * u[i][k - 1];
u[i][k - 1] = -1 * sn * u[i][j] +
cs * u[i][k - 1];
u[i][j] = t;
}
}
}
else if (branch == 3)
{
// calculate scale and shift
double max1 = Math.Max(Math.Abs(s[p - 1]),
Math.Abs(s[p - 2]));
double max2 = Math.Max(max1, Math.Abs(z[p - 2]));
double max3 = Math.Max(max2, Math.Abs(s[k]));
double max4 = Math.Max(max3, Math.Abs(z[k]));
double scale = max4;
double sp = s[p - 1] / scale;
double spm1 = s[p - 2] / scale;
double epm1 = z[p - 2] / scale;
double sk = s[k] / scale;
double zk = z[k] / scale;
double b = ((spm1 + sp) * (spm1 - sp) +
epm1 * epm1) / 2.0;
double c = (sp * epm1) * (sp * epm1);
double shift = 0.0;
if ((b != 0.0) || (c != 0.0)) // bug fixed
{
shift = Math.Sqrt(b * b + c);
if (b "lt" 0.0)
shift = -1 * shift;
shift = c / (b + shift);
}
double f = (sk + sp) * (sk - sp) + shift;
double g = sk * zk;
// chase out intermediate and trailing zeros
for (int j = k; j "lt" p - 1; ++j)
{
double t = LocalHypot(f, g);
double cs = f / t;
double sn = g / t;
if (j != k)
z[j - 1] = t;
f = cs * s[j] + sn * z[j];
z[j] = cs * z[j] - sn * s[j];
g = sn * s[j + 1];
s[j + 1] = cs * s[j + 1];
for (int i = 0; i "lt" n; ++i)
{
t = cs * v[i][j] + sn * v[i][j + 1];
v[i][j + 1] = -1 * sn * v[i][j] +
cs * v[i][j + 1];
v[i][j] = t;
}
t = LocalHypot(f, g);
cs = f / t;
sn = g / t;
s[j] = t;
f = cs * z[j] + sn * s[j + 1];
s[j + 1] = -1 * sn * z[j] + cs * s[j + 1];
g = sn * z[j + 1];
z[j + 1] = cs * z[j + 1];
if (j "lt" m - 1)
{
for (int i = 0; i "lt" m; ++i)
{
t = cs * u[i][j] + sn * u[i][j + 1];
u[i][j + 1] = -1 * sn * u[i][j] +
cs * u[i][j + 1];
u[i][j] = t;
}
}
} // for j
z[p - 2] = f;
++iter;
} // branch == 3
else if (branch == 4)
{
// make singular values positive
if (s[k] "lte" 0.0)
{
if (s[k] "lt" 0.0)
s[k] = -1 * s[k];
else
s[k] = 0.0;
for (int i = 0; i "lte" pp; ++i)
v[i][k] = -1 * v[i][k];
}
// sort singular values
while (k "lt" pp)
{
if (s[k] "gte" s[k + 1])
break;
double t = s[k];
s[k] = s[k + 1];
s[k + 1] = t;
if (k "lt" n - 1)
for (int i = 0; i "lt" n; ++i)
{
t = v[i][k + 1];
v[i][k + 1] = v[i][k];
v[i][k] = t;
}
if (k "lt" m - 1)
for (int i = 0; i "lt" m; ++i)
{
t = u[i][k + 1];
u[i][k + 1] = u[i][k];
u[i][k] = t;
}
++k;
} // while
iter = 0;
--p;
} // branch == 4
else
{
throw new Exception("Unhandled logic branch");
}
} // master while-loop
// return the computed items as out parameters
U = LocalMatNegate(u); // sync with NumPy output
S = s;
Vh = LocalMatNegate(LocalMatTrans(v)); // sync
// ----------------------------------------------------
// local nested functions
// ----------------------------------------------------
static double[][] LocalMatCreate(int rows, int cols)
{
double[][] result = new double[rows][];
for (int i = 0; i "lt" rows; ++i)
result[i] = new double[cols];
return result;
}
// ----------------------------------------------------
static double[][] LocalMatCopy(double[][] mat)
{
int rows = mat.Length; int cols = mat[0].Length;
double[][] result = LocalMatCreate(rows, cols);
for (int i = 0; i "lt" rows; ++i)
for (int j = 0; j "lt" cols; ++j)
result[i][j] = mat[i][j];
return result;
}
// ----------------------------------------------------
static double LocalHypot(double a, double b) // local
{
// wacky sqrt(a^2 + b^2) used by most svd() funcs
double xabs = Math.Abs(a); double yabs = Math.Abs(b);
double min, max;
if (xabs "lt" yabs) {
min = xabs; max = yabs;
}
else {
min = yabs; max = xabs;
}
if (min == 0)
return max;
else {
double u = min / max;
return max * Math.Sqrt(1 + u * u);
}
}
// ----------------------------------------------------
static double[][] LocalMatTrans(double[][] m)
{
int r = m.Length; int c = m[0].Length;
double[][] result = LocalMatCreate(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[][] LocalMatNegate(double[][] m)
{
int r = m.Length; int c = m[0].Length;
double[][] result = LocalMatCreate(r, c);
for (int i = 0; i "lt" r; ++i)
for (int j = 0; j "lt" c; ++j)
result[i][j] = -1 * m[i][j];
return result;
}
// ----------------------------------------------------
} // MatDecompSVD()
} // Program
} // ns
The NumPy version:
# svd_demo.py
import numpy as np
print("Singular value decomp with np.linalg.svd")
np.set_printoptions(precision=4, suppress=True)
A = np.array([[1,2,3],
[5,0,2],
[8,5,4],
[1,0,9]])
print("Source A: ")
print(A)
U, S, Vh = np.linalg.svd(A, full_matrices=False)
print("U: ")
print(U)
print("S: ")
print(S)
print("Vh: ")
print(Vh)
USVh = np.matmul(U, np.matmul(np.diag(S), Vh))
print("U * diag(S) * Vh: ")
print(USVh)
print("End NumPy demo ")
.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.