Bottom line: I implemented SVD using the Householder + QR algorithm, with JavaScript. The implementation works but it’s not as stable as my implementation that uses the Jacobi algorithm. In other words, the code in this post is interesting, but it’s not good for use in a production environment.
If you have a matrix A, the singular value decomposition (SVD) gives you a matrix U, a vector s, and a matrix Vh, so that A = U * S * Vh. The s is a vector like [2, 4, 6], and S is a square matrix with the s values on the diagonal, 0s off the diagonal.
SVD is used in dozens of important software algorithms.
SVD is an extremely complicated topic, and computing the SVD of a matrix is one of the most challenging problems in all of numerical programming. Many researchers worked for many years on the problem, and came up with many different versions of SVD implementations that vary in accuracy, complexity, and robustness. The standard SVD implementation in LAPACK is over 10,000 lines of Fortran code.
I put together an implementation of SVD using the JavaScript language (running in a node.js environemnt). The first part of my demo output is:
Source (tall) matrix: 1.0 2.0 3.0 4.0 5.0 0.0 -3.0 5.0 -7.0 9.0 2.0 0.0 -2.0 0.0 -2.0 4.0 -1.0 5.0 6.0 1.0 3.0 6.0 8.0 2.0 2.0 5.0 -2.0 4.0 -4.0 3.0 Computing SVD Done U = 0.319267 0.295047 0.358160 0.465731 0.652040 0.625260 -0.614481 0.189087 0.174205 -0.138090 -0.129989 0.032917 -0.342862 -0.160369 0.572690 0.284056 0.494504 -0.490845 0.494308 -0.381366 0.486033 0.495867 0.263713 -0.648653 -0.103790 0.416300 -0.209426 -0.638701 -0.248874 0.267559 s = 15.967661 12.793149 6.297366 5.706879 2.478679 Vh = 0.296544 0.035215 0.708796 -0.130799 0.625557 0.217254 0.416871 0.261754 0.803400 -0.255056 -0.745282 0.555722 -0.030757 0.039094 0.365040 -0.187161 -0.609726 -0.196995 0.579569 0.467438 0.523820 0.379983 -0.623971 0.003540 0.438033
Behind the scenes, I verified my implementation by feeding the same source matrix A to the NumPy library np.linalg.svd() function to make sure I got the same results. I also constructed the S matrix from the s vector, and then verified that A = U * S * Vh.
For the second part of my demo, I used a “wide” matrix that has more columns than rows – the SVD algorithm works differently for tall matrices than it does for wide matrices:
Source (wide) matrix: 3.0 1.0 1.0 0.0 5.0 -1.0 3.0 1.0 -2.0 4.0 0.0 2.0 2.0 1.0 -3.0 Computing SVD Done U = -0.727149 0.195146 -0.658158 -0.621912 -0.593195 0.511219 0.290654 -0.781049 -0.552705 s = 7.639218 4.023860 3.232785 Vh = -0.204149 -0.263322 -0.100502 0.200868 -0.915716 0.292911 -0.781970 -0.487131 0.100734 0.235122 -0.768902 -0.071118 -0.387390 -0.487240 0.127506
My JavaScript implementation is intended mostly for machine learning scenarios where accuracy to 1.0e-8 is good enough, as opposed to scientific programming where accuracy to 1.0e-15 is typical.
I worked on this SVD implementation off-and-on for several years. I finally got some extended uninterrupted time and coded up a demo but the implementation took several weeks to write. I used over a dozen major resources. The implementation here has all normal error-checking removed to keep the number of lines of code as small as possible. The implementation is a combination of several different SVD algorithms, mostly using QR decomposition (Householder algorithm), and eigen decomposition for a symmetric algorithms (which also uses QR decomposition).
I’ve written some complex code over the years but this SVD implementation is one of the most complex functions I’ve ever implemented. I consider the code presented in this blog post to be one of the significant accomplishments of my software engineering career. Unfortunately, this implementation is not as stable as my older version the uses the Jacobi algorithm.
The SVD decomposition algorithm has many patterns that are not immediately obvious.
Every cover of Playboy magazine has an image of a representation of the bunny logo. The logo is plain and obvious on most covers, but some covers have the logo cleverly hidden. Here are three nice examples. Left: December 1966 (look in the model’s hair). Center: November 1967 (the model’s dress represents the logo — very clever). Right: May 1968 (look at the flying ribbons).
Demo program. Very long, very complicated. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols. (My blog editor chokes on symbols).
// svd.js
// singular value decomposition via QR Householder
// node.js
let FS = require("fs") // for matLoad()
class SVD
{
// --------------------------------------------------------
static matDecomp(A, tol = 1.0e-12)
{
let m = A.length; let n = A[0].length;
if (m "gte" n) {
let QR = SVD.matDecompQR(A);
let Q = QR[0]; let R = QR[1];
let RtR = SVD.matProduct(SVD.matTranspose(R), R);
let EV = SVD.matDecompEigen(RtR);
let evals = EV[0]; let V = EV[1];
let sortIdxs = SVD.vecReversed(SVD.argsort(evals));
evals = SVD.vecSortedUsingIdxs(evals, sortIdxs);
V = SVD.matSortedUsingIdxs(V, sortIdxs);
let ss = SVD.vecSqrt(SVD.vecClip(evals, 0.0));
let idxs = SVD.vecIndicesWhereNonNegative(ss, tol);
let ssNonZero = SVD.vecRemoveZeros(ss, idxs);
let VNonZero = SVD.matRemoveCols(V, idxs);
let tmp1 = SVD.matProduct(R, VNonZero);
let tmp2 = SVD.matDivVector(tmp1, ssNonZero);
let UU = SVD.matProduct(Q, tmp2);
let result = [];
result[0] = UU;
result[1] = ssNonZero;
result[2] = SVD.matTranspose(VNonZero);
return result;
}
else if (m "lt" n) {
// use the transpose trick
let UsVh = SVD.matDecomp(SVD.matTranspose(A));
let result = [];
result[0] = SVD.matTranspose(UsVh[2]);
result[1] = UsVh[1];
result[2] = SVD.matTranspose(UsVh[0]);
return result;
}
} // matDecomp()
// --------------------------------------------------------
static matDecompQR(A)
{
// reduced QR decomp using Householder reflections
let m = A.length; let n = A[0].length;
let QQ = SVD.matIdentity(m); // working Q
let RR = SVD.matCopy(A); // working R
for (let k = 0; k "lt" n; ++k) {
let x = SVD.matColToEnd(RR, k); // RR[k:, k]
let normx = SVD.vecNorm(x);
if (Math.abs(normx) "lt" 1.0e-12) continue;
let v = SVD.vecCopy(x);
let sign = 0.0;
if (x[0] "gte" 0.0) sign = 1.0; else sign = -1.0;
v[0] += sign * normx;
let normv = SVD.vecNorm(v);
for (let i = 0; i "lt" v.length; ++i)
v[i] /= normv;
// apply reflection to R
let tmp1 = SVD.matRowColToEnd(RR, k); // RR[k:, k:]
let tmp2 = SVD.vecMatProd(v, tmp1);
let tmp3 = SVD.vecOuter(v, tmp2);
let B = SVD.matScalarMult(2.0, tmp3);
SVD.matSubtractCorner(RR, k, B); // R[k:, k:] -= B
// accumulate Q
// QQ[:, k:]
let tmp4 = SVD.matExtractAllRowsColToEnd(QQ, k);
let tmp5 = SVD.matVecProd(tmp4, v);
let tmp6 = SVD.vecOuter(tmp5, v);
let C = SVD.matScalarMult(2.0, tmp6);
SVD.matSubtractRows(QQ, k, C); // QQ[:, k:] -= C
} // for k
// return parts of QQ and RR
let result = [];
result[0] = SVD.matExtractFirstCols(QQ, n); // Q[:, :n]
result[1] = SVD.matExtractFirstRows(RR, n); // R[:n, :]
return result;
} // matDecompQR()
// --------------------------------------------------------
static matDecompEigen(A, tol = 1.0e-12, maxIter = 5000)
{
// eigen-decomp of symmetric matrix using QR
let m = A.length; let n = A[0].length;
let AA = SVD.matCopy(A); // don't destroy A
let V = SVD.matIdentity(n);
let iter = 0;
while (iter "lt" maxIter) {
let QR = SVD.matDecompQR(AA);
let Q = QR[0]; let R = QR[1];
let Anew = SVD.matProduct(R, Q);
V = SVD.matProduct(V, Q);
// check convergence (off-diagonal norm)
let tmp1 = SVD.matDiag(Anew);
let tmp2 = SVD.matCreateDiag(tmp1);
let tmp3 = SVD.matSubtract(Anew, tmp2);
let norm = SVD.matNorm(tmp3);
if (norm "lt" tol) {
AA = Anew; break;
}
AA = Anew;
++iter;
}
if (iter == maxIter)
console.log("Warn: exceeded maxIter ");
let result = [];
result[0] = SVD.matDiag(AA); // eigen vals
result[1] = V; // eigen vecs
return result;
} // matDecompEigen()
// --------------------------------------------------------
static argsort(vec)
{
let n = vec.length;
let augmented = [];
for (let i = 0; i "lt" n; ++i) {
augmented[i] = [vec[i], i]; // add indices to values
}
augmented.sort(function(a, b){return a[0] - b[0]});
let result = [];
for (let i = 0; i "lt" n; ++i) {
result[i] = augmented[i][1]; // fetch the indexes
}
return result; // 2, 4, 3, 0, 1
}
// --------------------------------------------------------
static matMake(nRows, nCols, val)
{
let result = [];
for (let i = 0; i "lt" nRows; ++i) {
result[i] = [];
for (let j = 0; j "lt" nCols; ++j) {
result[i][j] = val;
}
}
return result;
}
// --------------------------------------------------------
static matIdentity(n)
{
let result = SVD.matMake(n, n, 0.0);
for (let i = 0; i "lt" n; ++i)
result[i][i] = 1.0;
return result;
}
// --------------------------------------------------------
static matCopy(A)
{
let nRows = A.length;
let nCols = A[0].length;
let result = SVD.matMake(nRows, nCols, 0.0);
for (let i = 0; i "lt" nRows; ++i)
for (let j = 0; j "lt" nCols; ++j)
result[i][j] = A[i][j];
return result
}
// --------------------------------------------------------
static matDiag(A)
{
// return diag elements of A as a vector
let m = A.length; let n = A[0].length; // m == n
let result = SVD.vecMake(m, 0.0);
for (let i = 0; i "lt" m; ++i)
result[i] = A[i][i];
return result;
}
// --------------------------------------------------------
static matCreateDiag(v)
{
// create a matrix using v as diagonal elements
let n = v.length;
let result = SVD.matMake(n, n, 0.0);
for (let i = 0; i "lt" n; ++i)
result[i][i] = v[i];
return result;
}
// --------------------------------------------------------
static matTranspose(A)
{
let nRows = A.length; let nCols = A[0].length;
let result = SVD.matMake(nCols, nRows, 0.0);
for (let i = 0; i "lt" nRows; ++i)
for (let j = 0; j "lt" nCols; ++j)
result[j][i] = A[i][j];
return result;
}
// --------------------------------------------------------
static matSubtract(A, B)
{
let m = A.length; let n = A[0].length;
let result = SVD.matMake(m, n, 0.0);
for (let i = 0; i "lt" m; ++i)
for (let j = 0; j "lt" n; ++j)
result[i][j] = A[i][j] - B[i][j];
return result;
}
// --------------------------------------------------------
static matColToEnd(A, k)
{
// last part of col k, from [k,k] to end
let m = A.length; let n = A[0].length;
let result = SVD.vecMake(m - k, 0.0);
for (let i = 0; i "lt" m - k; ++i)
result[i] = A[i + k][k];
return result;
}
// --------------------------------------------------------
static matRowColToEnd(A, k)
{
// block of A, row k to end, col k to end, A[k:, k:]
let m = A.length; let n = A[0].length;
let result = SVD.matMake(m - k, n - k, 0.0);
for (let i = 0; i "lt" m - k; ++i)
for (let j = 0; j "lt" n - k; ++j)
result[i][j] = A[i + k][j + k];
return result;
}
// --------------------------------------------------------
static matVecProd(A, v)
{
// A * v returns a regular 1D vector
let nRows = A.length;
let nCols = A[0].length;
let n = v.length; // assume nCols = n
let result = SVD.vecMake(nRows, 0.0);
for (let i = 0; i "lt" nRows; ++i)
for (let k = 0; k "lt" nCols; ++k)
result[i] += A[i][k] * v[k];
return result;
}
// --------------------------------------------------------
static vecMatProd(v, A)
{
// v * A
let m = A.length; let n = A[0].length;
let result = SVD.vecMake(n, 0.0);
for (let j = 0; j "lt" n; ++j)
for (let i = 0; i "lt" m; ++i)
result[j] += v[i] * A[i][j];
return result;
}
// --------------------------------------------------------
static matScalarMult(x, A)
{
// x * A element-wise
let m = A.length; let n = A[0].length;
let result = SVD.matMake(m, n, 0.0);
for (let i = 0; i "lt" m; ++i)
for (let j = 0; j "lt" n; ++j)
result[i][j] = x * A[i][j];
return result;
}
// --------------------------------------------------------
static matSubtractCorner(A, k, C)
{
// subtract elements in C from lower right A
// A[k:, k:] -= C
let m = A.length; let n = A[0].length;
for (let i = 0; i "lt" m - k; ++i)
for (let j = 0; j "lt" n - k; ++j)
A[i + k][j + k] -= C[i][j];
return; // modify A in-place
}
// --------------------------------------------------------
static matSubtractRows(A, k, C)
{
// A[:, k:] -= C
let m = A.length; let n = A[0].length;
for (let i = 0; i "lt" m; ++i)
for (let j = 0; j "lt" n - k; ++j)
A[i][j + k] -= C[i][j];
return; // modify A in-place
}
// ------------------------------------------------------
static matExtractAllRowsColToEnd(A, k)
{
// A[:, k:]
let m = A.length; let n = A[0].length;
let result = SVD.matMake(m, n - k, 0.0);
for (let i = 0; i "lt" m; ++i)
for (let j = 0; j "lt" n - k; ++j)
result[i][j] = A[i][j + k];
return result;
}
// -----------------------------------------------------
static matExtractFirstCols(A, k)
{
// A[:, :n]
let m = A.length; let n = A[0].length;
let result = SVD.matMake(m, k, 0.0);
for (let i = 0; i "lt" m; ++i)
for (let j = 0; j "lt" k; ++j)
result[i][j] = A[i][j];
return result;
}
// --------------------------------------------------------
static matExtractFirstRows(A, k)
{
// A[:n, :]
let m = A.length; let n = A[0].length;
let result = SVD.matMake(k, n, 0.0);
for (let i = 0; i "lt" k; ++i)
for (let j = 0; j "lt" n; ++j)
result[i][j] = A[i][j];
return result;
}
// --------------------------------------------------------
static matProduct(A, B)
{
let aRows = A.length; let aCols = A[0].length;
let bRows = B.length; let bCols = B[0].length;
let result = SVD.matMake(aRows, bCols, 0.0);
for (let i = 0; i "lt" aRows; ++i)
for (let j = 0; j "lt" bCols; ++j)
for (let k = 0; k "lt" aCols; ++k)
result[i][j] += A[i][k] * B[k][j];
return result;
}
// --------------------------------------------------------
static matRemoveCols(A, idxs)
{
// remove the columns specified in idxs
let m = A.length; let n = A[0].length;
let countNonzeos = 0;
for (let i = 0; i "lt" idxs.length; ++i)
if (idxs[i] == 1) ++countNonzeos;
let result = SVD.matMake(m, countNonzeos, 0.0);
let k = 0; // destination col
for (let j = 0; j "lt" n; ++j) {
if (idxs[j] == 0) continue; // skip this col
for (let i = 0; i "lt" m; ++i) {
result[i][k] = A[i][j];
}
++k;
}
return result;
}
// --------------------------------------------------------
static matDivVector(A, v)
{
// v must be all non-zero
// divide each row of A by v, element-wise
let m = A.length; let n = A[0].length;
let result = SVD.matMake(m, n, 0.0);
for (let i = 0; i "lt" m; ++i)
for (let j = 0; j "lt" n; ++j)
result[i][j] = A[i][j] / v[j];
return result;
}
// --------------------------------------------------------
static matNorm(A)
{
// treat A as one big vector
let m = A.length; let n = A[0].length;
let sum = 0.0;
for (let i = 0; i "lt" m; ++i)
for (let j = 0; j "lt" n; ++j)
sum += A[i][j] * A[i][j];
return Math.sqrt(sum);
}
// --------------------------------------------------------
static matSortedUsingIdxs(A, idxs)
{
// arrange columns of A using order in idxs
let m = A.length; let n = A[0].length;
let result = SVD.matMake(m, n, 0.0);
for (let i = 0; i "lt" m; ++i) {
for (let j = 0; j "lt" n; ++j) {
let srcCol = idxs[j];
result[i][j] = A[i][srcCol];
}
}
return result;
}
// --------------------------------------------------------
static vecSortedUsingIdxs(v, idxs)
{
// arrange values in v using order specified in idxs
let n = v.length;
let result = SVD.vecMake(n, 0.0);
for (let i = 0; i "lt" n; ++i)
result[i] = v[idxs[i]];
return result;
}
// --------------------------------------------------------
static vecMake(n, val)
{
let result = [];
for (let i = 0; i "lt" n; ++i) {
result[i] = val;
}
return result;
}
// --------------------------------------------------------
static vecNorm(vec)
{
let sum = 0.0;
let n = vec.length;
for (let i = 0; i "lt" n; ++i)
sum += vec[i] * vec[i];
return Math.sqrt(sum);
}
// --------------------------------------------------------
static vecDot(v1, v2)
{
let n = v1.length; // assume len(v1) == len(v2)
let sum = 0.0;
for (let i = 0; i "lt" n; ++i)
sum += v1[i] * v2[i];
return sum;
}
// --------------------------------------------------------
static vecCopy(v)
{
let result = SVD.vecMake(v.length, 0.0);
for (let i = 0; i "lt" v.length; ++i)
result[i] = v[i];
return result;
}
// --------------------------------------------------------
static vecOuter(v1, v2)
{
// vector outer product
let result = SVD.matMake(v1.length, v2.length);
for (let i = 0; i "lt" v1.length; ++i)
for (let j = 0; j "lt" v2.length; ++j)
result[i][j] = v1[i] * v2[j];
return result;
}
// --------------------------------------------------------
static vecReversed(v)
{
// return vector of v elements in reversed order
let n = v.length;
let result = SVD.vecMake(n, 0.0);
for (let i = 0; i "lt" n; ++i)
result[i] = v[n - 1 - i];
return result;
}
// ------------------------------------------------------
static vecClip(v, minVal)
{
// used to make sure no negative values
let n = v.length;
let result = SVD.vecMake(n, 0.0);
for (let i = 0; i "lt" n; ++i) {
if (v[i] "lt" minVal)
result[i] = minVal;
else
result[i] = v[i];
}
return result;
}
// ------------------------------------------------------
static vecSqrt(v)
{
// sqrt each element in v
// assumes v has no negative values
let n = v.length;
let result = SVD.vecMake(n, 0.0);
for (let i = 0; i "lt" n; ++i)
result[i] = Math.sqrt(v[i]);
return result;
}
// ------------------------------------------------------
static vecIndicesWhereNonNegative(v, tol)
{
// indices of non-negative values (within tolerance)
let n = v.length;
let result = SVD.vecMake(n, 0.0);
for (let i = 0; i "lt" n; ++i) {
if (Math.abs(v[i]) "gt" tol)
result[i] = 1;
}
return result;
}
// ------------------------------------------------------
static vecRemoveZeros(v, idxs)
{
// if idx = 1, val is non-zero, if 0 val is near zero
let n = v.length;
let result = [];
for (let i = 0; i "lt" n; ++i) {
if (idxs[i] == 1)
result.push(v[i]);
}
return result;
}
// ------------------------------------------------------
} // class SVD
// ----------------------------------------------------------
// helper functions for main()
// matShow, vecShow, matLoad (not used)
// ----------------------------------------------------------
function vecShow(vec, dec, wid, nl)
{
let small = 1.0 / Math.pow(10, dec);
for (let i = 0; i "lt" vec.length; ++i) {
let x = vec[i];
if (Math.abs(x) "lt" small) x = 0.0 // avoid -0.00
let xx = x.toFixed(dec);
let s = xx.toString().padStart(wid, ' ');
process.stdout.write(s);
process.stdout.write(" ");
}
if (nl == true)
process.stdout.write("\n");
}
// ----------------------------------------------------------
function matShow(A, dec, wid)
{
let small = 1.0 / Math.pow(10, dec);
let nr = A.length;
let nc = A[0].length;
for (let i = 0; i "lt" nr; ++i) {
for (let j = 0; j "lt" nc; ++j) {
let x = A[i][j];
if (Math.abs(x) "lt" small) x = 0.0;
let xx = x.toFixed(dec);
let s = xx.toString().padStart(wid, ' ');
process.stdout.write(s);
process.stdout.write(" ");
}
process.stdout.write("\n");
}
}
// ----------------------------------------------------------
function matLoad(fn, delimit, usecols, comment)
{
// efficient but mildly complicated
let all = FS.readFileSync(fn, "utf8"); // giant string
all = all.trim(); // strip final crlf in file
let lines = all.split("\n"); // array of lines
// count number non-comment lines
let nRows = 0;
for (let i = 0; i "lt" lines.length; ++i) {
if (!lines[i].startsWith(comment))
++nRows;
}
let nCols = usecols.length;
// let result = matMake(nRows, nCols, 0.0);
let result = [];
for (let i = 0; i "lt" nRows; ++i) {
result[i] = [];
for (let j = 0; j "lt" nCols; ++j) {
result[i][j] = 0.0;
}
}
let r = 0; // into lines
let i = 0; // into result[][]
while (r "lt" lines.length) {
if (lines[r].startsWith(comment)) {
++r; // next row
}
else {
let tokens = lines[r].split(delimit);
for (let j = 0; j "lt" nCols; ++j) {
result[i][j] = parseFloat(tokens[usecols[j]]);
}
++r;
++i;
}
}
return result;
}
// ----------------------------------------------------------
function main()
{
console.log("\nSingular value decomposition " +
"(QR Householder) using JavaScript ");
let A = [];
A[0] = [ 1, 2, 3, 4, 5 ];
A[1] = [ 0, -3, 5, -7, 9 ];
A[2] = [ 2, 0, -2, 0, -2 ];
A[3] = [ 4, -1, 5, 6, 1 ];
A[4] = [ 3, 6, 8, 2, 2 ];
A[5] = [ 5, -2, 4, -4, 3 ];
console.log("\nSource (tall) matrix: ");
matShow(A, 1, 6);
console.log("\nComputing SVD ");
let USV = SVD.matDecomp(A);
let U = USV[0];
let s = USV[1];
let Vh = USV[2];
console.log("Done ");
console.log("\nU = ");
matShow(U, 6, 11);
console.log("\ns = ");
vecShow(s, 6, 11, true)
console.log("\nVh = ");
matShow(Vh, 6, 11);
console.log("\n============================ ");
let B = [];
B[0] = [ 3, 1, 1, 0, 5 ];
B[1] = [ -1, 3, 1, -2, 4 ];
B[2] = [ 0, 2, 2, 1, -3 ];
console.log("\nSource (wide) matrix: ");
matShow(B, 1, 6);
console.log("\nComputing SVD ");
USV = SVD.matDecomp(B);
U = USV[0];
s = USV[1];
Vh = USV[2];
console.log("Done ");
console.log("\nU = ");
matShow(U, 6, 11);
console.log("\ns = ");
vecShow(s, 6, 11, true)
console.log("\nVh = ");
matShow(Vh, 6, 11);
console.log("\nEnd demo");
}
main();
.NET Test Automation Recipes
Software Testing
SciPy Programming Succinctly
Keras Succinctly
R Programming
2026 Visual Studio Live
2025 Summer MLADS Conference
2025 DevIntersection Conference
2025 Machine Learning Week
2025 Ai4 Conference
2025 G2E Conference
2025 iSC West Conference
You must be logged in to post a comment.