I was sitting on a six-hour flight from Boston to Seattle. I was bored. I figured I’d pass the time by refactoring some matrix pseudo-inverse code from C# to JavaScript.

In ordinary arithmetic, the inverse of 4 is 0.25 because 4 * 0.25 = 1. In matrix arithmetic, the inverse of a square n-by-n matrix A is an n-by-n matrix Ainv where A * Ainv = I. (Also Ainv * A = I). The * is matrix multiplication, the I is the n-by-n identity matrix which has 1s on the upper-left to lower-right diagonal, 0s elsewhere.

For example, if A =

 3.0   2.0
 0.0   4.0

the inverse of A is Ainv =

 0.3333  -0.1667
 0.0000   0.2500

because A * Ainv = Ainv * A =

 1.0   0.0
 0.0   1.0

But for matrices that aren’t square, such as a 3-by-5 or a 4-by-2 matrix, there’s no inverse. But it’s sometimes useful to define a pseudo-inverse. There are several different types of pseudo-inverse, but the most common is called the Moore-Penrose pseudo-inverse.

For a matrix A, the MP pseudo-inverse is a matrix Apinv where A * Apinv * A = A. Computing a pseudo-inverse is tricky. The technique I usually prefer is the singular value decomposition (SVD) method. SVD breaks a matrix A down into three components: a matrix U, a matrix Vh, and a vector s. If you compute U * S * Vh you get the original matrix A. Here S is a square matrix with the elements of vector s on the diagonal, 0s elsewhere. (The values in the s vector are called the singular values of the matrix.

Computing an SVD decomposition is very difficult. Based on my experience (but no hard data), the two most common techniques are the Jacobi method and the QR decomposition method. If you new to matrix arithmetic, the details can be overwhelming at first.

If you can compute the U, Vh, and S from a source matrix A, then the pseudo-inverse of A is V * Sinv * Uh. Here V is the inverse of Vh, which fortunately is just the transpose (rows and columns switched) of Vh. The Uh is the transpose of U, which again, fortunately is just the transpose of U. The Sinv is the S matrix but each value on the diagonal is 1/s instead of s.

To summarize, in order to find the MP pseudo-inverse of a matrix A, perform SVD decomposition to get U, Vh, and s. Compute V = transpose(Vh), Uh = transpose(U), Sinv (as described above), and then compute Apinv = V * Sinv * U).

The code shows how easy it is:

function matPseudoInverse(M)
{
  let decomp = svdJacobi(M);
  let U = decomp[0];
  let s = decomp[1];
  let Vh = decomp[2];

  let Sinv = matDiag(s); 
  for (let i = 0; i "lt" Sinv.length; ++i)
    Sinv[i][i] = 1.0 / Sinv[i][i];

  let V = matTranspose(Vh);
  let Uh = matTranspose(U);
  let result = matMult(V, matMult(Sinv, Uh));
  return result;
}

Unfortunately, computing the SVD decomposition is very, very difficult. Luckily, I had previously implemented SVD using the Jacobi technique in JavaScript (based on code in the GNU scientific library) and so all the hard work was done. All I had to do was reuse the SVD code. Code reuse is a good thing.

Here’s the output of my demo program:

Begin matrix pseudo-inverse (via Jacobi SVD) using JavaScript

Source matrix A:
 1.00   2.00   3.00
 5.00   0.00   2.00
 8.00   5.00   4.00
 1.00   0.00   9.00

Pseudo-inverse A =
  -0.0784     0.1706     0.0314    -0.0257
   0.1568    -0.2390     0.1373    -0.0602
   0.0287    -0.0091    -0.0115     0.1087

Verifying A * Apinv * A = 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

Verified equal

source matrix B:
 1.00   5.00   8.00   1.00
 2.00   0.00   5.00   0.00
 3.00   2.00   4.00   9.00

Pseudo-inverse B =
  -0.0784     0.1568     0.0287
   0.1706    -0.2390    -0.0091
   0.0314     0.1373    -0.0115
  -0.0257    -0.0602     0.1087

Verifying B * Bpinv * B = B
   1.0000     5.0000     8.0000     1.0000
   2.0000     0.0000     5.0000     0.0000
   3.0000     2.0000     4.0000     9.0000

Verified equal

End demo

An interesting way to spend time on a flight from Boston to Seattle.

Code reuse is a key technique in software development. Prop reuse in science fiction movies is a key technique too.

Right: The Forbidden Planet uniforms were used again in an episode of the old “Twilight Zone” TV series titled “The Little People” (1962). A U.S. space explorer lands on a planet inhabited by tiny people. He goes insane and forces the little people to worship him as a god. The story ends with a pair of gigantic alien space explorers (wearing the Forbidden Planet uniforms) landing on the planet. It does not end well for the pseudo-god, now-tiny explorer. A classic episode.

Demo program. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols (my blog editor chokes on symbols).

// matrix_pseudoinverse_svd_jacobi.js
// pseudo-inverse using svd (Jacobi method)
// node.js ES6

let FS = require('fs');  // for matLoad()

// ----------------------------------------------------------

function matPseudoInverse(M)
{
  let decomp = svdJacobi(M);
  let U = decomp[0];
  let s = decomp[1];
  let Vh = decomp[2];

  let Sinv = matDiag(s); 
  for (let i = 0; i "lt" Sinv.length; ++i)
    Sinv[i][i] = 1.0 / Sinv[i][i];

  let V = matTranspose(Vh);
  let Uh = matTranspose(U);
  let tmp = matMult(Sinv, Uh);
  let result = matMult(V, tmp);
  return result;
}

// ----------------------------------------------------------

function svdJacobi(M)
{
  // A = U * S * Vh
  // returns an array where result[0] = U matrix,
  //  result[1] = s vector, result[2] = Vh matrix

  // see github.com/ampl/gsl/blob/master/linalg/svd.c
  // numpy svd() full_matrices = False version

  let DBL_EPSILON = 1.0e-15;

  let A = matCopy(M); // working U
  let m = A.length;
  let n = A[0].length;
  let Q = matIdentity(n); // working V
  let t = vecMake(n, 0.0);  // working s

  // initialize the rotation and sweep counters
  let count = 1;
  let sweep = 0;

  let tolerance = 10 * m * DBL_EPSILON; // heuristic

  // always do at least 12 sweeps
  sweepmax = Math.max(5 * n, 12); // heuristic

  // store the column error estimates in t, for use
  // during orthogonalization

  for (let j = 0; j "lt" n; ++j) {
    let cj = matGetColumn(A, j);
    let sj = vecNorm(cj);
    t[j] = DBL_EPSILON * sj;
  }

  // orthogonalize A by plane rotations
  while (count "gt" 0 && sweep "lte" sweepmax) {
    // initialize rotation counter
    count = n * (n - 1) / 2;

    for (let j = 0; j "lt" n - 1; ++j) {
      for (let k = j + 1; k "lt" n; ++k) {
        let cosine = 0.0; let sine = 0.0;

        let cj = matGetColumn(A, j);
        let ck = matGetColumn(A, k);

        let p = 2.0 * vecDot(cj, ck);
        let a = vecNorm(cj);
        let b = vecNorm(ck);

        let q = a * a - b * b;
        let v = hypot(p, q);

        // test for columns j,k orthogonal,
        // or dominant errors 
        let abserr_a = t[j];
        let abserr_b = t[k];

        let sorted = (a "gte" b);
        let orthog = (Math.abs(p) "lte" 
          tolerance * (a * b));
        let noisya = (a "lt" abserr_a);
        let noisyb = (b "lt" abserr_b);

        if (sorted && (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 (let i = 0; i "lt" m; ++i) {
          let Aik = A[i][k];
          let 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 (Vh)
        for (let i = 0; i "lt" n; ++i) {
          let Qij = Q[i][j];
          let 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
  let prevNorm = -1.0;

  for (let j = 0; j "lt" n; ++j) {
    let column = matGetColumn(A, j);
    let norm = vecNorm(column);

    // determine if singular value is zero
    if (norm == 0.0 || prevNorm == 0.0
      || (j "gt" 0 && norm "lte" tolerance * prevNorm)) {
      t[j] = 0.0;
      for (let i = 0; i "lt" m; ++i)
        A[i][j] = 0.0;
      prevNorm = 0.0;
    }
    else {
      t[j] = norm;
      for (let i = 0; i "lt" m; ++i)
        A[i][j] = A[i][j] * 1.0 / norm;
      prevNorm = norm;
    }
  }

  if (count "gt" 0) {
    console.log("Jacobi iterations did not" +
      " converge");
  }

  let U = A;
  let s = t;
  let Vh = matTranspose(Q);

  if (m "lt" n) {
    U = matExtractFirstColumns(U, m);
    s = vecExtractFirst(s, m);
    Vh = matExtractFirstRows(Vh, m);
  } 

  let result = [];
  result[0] = U;
  result[1] = s;
  result[2] = Vh;

  return result;
} // svdJacobi()

// ----------------------------------------------------------

// === helper functions =====================================
//
// matMake, matCopy, matIdentity, matGetColumn,
// matExtractFirstColumns, matExtractFirstRows,
// matTranspose, matDiag, matMult, vecMake, vecNorm,
// vecDot, hypot, matShow, vecShow, matAreEqual, matLoad
//
// ==========================================================

function matMake(rows, cols, val)
{
  let result = [];
  for (let i = 0; i "lt" rows; ++i) {
    result[i] = [];
    for (let j = 0; j "lt" cols; ++j) {
      result[i][j] = val;
    }
  }
  return result;
}

// ----------------------------------------------------------

function matCopy(m)
{
  let rows = m.length;
  let cols = m[0].length;
  let result = matMake(rows, cols, 0.0);
  for (let i = 0; i "lt" rows; ++i)
    for (let j = 0; j "lt" cols; ++j)
      result[i][j] = m[i][j];
  return result;
}

// ----------------------------------------------------------

function matIdentity(n)
{
  let result = matMake(n, n, 0.0);
  for (let i = 0; i "lt" n; ++i)
    result[i][i] = 1.0;
  return result;
}

// ----------------------------------------------------------

function matGetColumn(m, j)
{
  let rows = m.length;
  let result = vecMake(rows, 0.0);
  for (let i = 0; i "lt" rows; ++i)
    result[i] = m[i][j];
  return result;
}

// ----------------------------------------------------------

function matExtract(src, r, c)
{
  // extract upper left corner
  let result = matMake(r, c, 0.0);
  for (let i = 0; i "lt" r; ++i)
    for (let j = 0; j "lt" c; ++j)
      result[i][j] = src[i][j];
  return result;
}

// ----------------------------------------------------------

function matExtractFirstColumns(src, n)
{
  let nRows = src.length;

  let result = matMake(nRows, n);
  for (let i = 0; i "lt" nRows; ++i)
    for (let j = 0; j "lt" n; ++j)
      result[i][j] = src[i][j];
  return result;
}

// ----------------------------------------------------------

function matExtractFirstRows(src, n)
{
  let nCols = src[0].length;
  let result = matMake(n, nCols);
  for (let i = 0; i "lt" n; ++i)
    for (let j = 0; j "lt" nCols; ++j)
      result[i][j] = src[i][j];
  return result;
}

// ----------------------------------------------------------

function matTranspose(m)
{
  let nRows = m.length;
  let nCols = m[0].length;
  let result = matMake(nCols, nRows, 0.0);
  for (let i = 0; i "lt" nRows; ++i)
    for (let j = 0; j "lt" nCols; ++j)
      result[j][i] = m[i][j];
  return result;  
}

// ----------------------------------------------------------

function matDiag(vec)
{
  let n = vec.length;
  let result = matMake(n, n, 0.0);
  for (let i = 0; i "lt" n; ++i)
    result[i][i] = vec[i];
  return result;
}

// ----------------------------------------------------------

function matMult(matA, matB)
{
  let aRows = matA.length;
  let aCols = matA[0].length;
  let bRows = matB.length;
  let bCols = matB[0].length;

  let result = 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] += matA[i][k] * matB[k][j];

  return result;
}

// ----------------------------------------------------------

function matNormColumns(m)
{
  // norm by columns
  let rows = m.length;
  let cols = m[0].length;
  let result = matCopy(m);

  for (let j = 0; j "lt" cols; ++j) {
    let col = matGetColumn(m, j);
    let norm = vecNorm(col);
    for (let i = 0; i "lt" rows; ++i)
      result[i][j] /= norm;
  }
  return result;
}

// ----------------------------------------------------------

function vecMake(n, val)
{
  let result = [];
  for (let i = 0; i "lt" n; ++i) {
    result[i] = val;
  }
  return result;
}

// ----------------------------------------------------------

function 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);
}

// ----------------------------------------------------------

function vecDot(v1, v2)
{
  let n = v1.length;
  let sum = 0.0;
  for (let i = 0; i "lt" n; ++i)
    sum += v1[i] * v2[i];
  return sum;
}

// ----------------------------------------------------------

function hypot(a, b)
{
  // wacky sqrt(a^2 + b^2) used by most SVD implementations
  let xabs = Math.abs(a); let yabs = Math.abs(b);
  let min, max;

  if (xabs "lt" yabs) {
    min = xabs; max = yabs;
  }
  else {
    min = yabs; max = xabs;
  }

  if (min == 0)
    return max;
  else {
    let u = min / max;
    return max * Math.sqrt(1 + u * u);
  }
}

// ----------------------------------------------------------

function vecScale(vec, factor)
{
  let n = vec.length;
  for (let i = 0; i "lt" n; ++i)
    vec[i] *= factor;
}

// ----------------------------------------------------------

function vecSetZero(vec)
{
  let n = vec.length;
  for (let i = 0; i "lt" n; ++i)
    vec[i] = 0.0;
}

// ----------------------------------------------------------

function vecCopy(vec)
{
  let n = vec.length;
  let result = vecMake(n, 0.0);
  for (let i = 0; i "lt" n; ++i)
    result[i] = vec[i];
  return result;
}

// ----------------------------------------------------------

function vecExtractFirst(vec, n)
{
  let result = vecMake(n, 0.0);
  for (let i = 0; i "lt" n; ++i)
    result[i] = vec[i];
  return result;
}

// ----------------------------------------------------------

function matShow(m, dec, wid)
{
  let rows = m.length;
  let cols = m[0].length;
  for (let i = 0; i "lt" rows; ++i) {
    for (let j = 0; j "lt" cols; ++j) {
      let v = m[i][j];
      if (Math.abs(v) "lt" 0.000001) v = 0.0  // avoid -0.00
      let vv = v.toFixed(dec);
      let s = vv.toString().padStart(wid, ' ');
      process.stdout.write(s);
      process.stdout.write("  ");
    }
    process.stdout.write("\n");
  }
}

// ----------------------------------------------------------

function vecShow(vec, dec, wid)
{
  for (let i = 0; i "lt" vec.length; ++i) {
    let x = vec[i];
    if (Math.abs(x) "lt" 0.000001) 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 matAreEqual(A, B, tol)
{
  let nRows = A.length; let nCols = A[0].length;
  for (let i = 0; i "lt" nRows; ++i) {
    for (let j = 0; j "lt" nCols; ++j) {
      if (Math.abs(A[i][j] - B[i][j]) "gt" tol) {
        return false;
      }
    }
  }
  return true;
}

// ----------------------------------------------------------

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;
  }
  nCols = usecols.length;
  let result = matMake(nRows, nCols, 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("\nBegin matrix pseudo-inverse " +
    "(via Jacobi SVD) using JavaScript ");

  let A = [];
  A[0] = [ 1.0, 2.0, 3.0 ];
  A[1] = [ 5.0, 0.0, 2.0 ];
  A[2] = [ 8.0, 5.0, 4.0 ];
  A[3] = [ 1.0, 0.0, 9.0 ];

  console.log("\nSource matrix A: ");
  matShow(A, 2, 5);

  let Apinv = matPseudoInverse(A);
  console.log("\nPseudo-inverse A = ");
  matShow(Apinv, 4, 9);

  console.log("\nVerifying A * Apinv * A = A ");
  let AAp = matMult(A, Apinv);
  let AApA = matMult(AAp, A);
  matShow(AApA, 4, 9);

  let checkA = matAreEqual(AApA, A, 1.0e-8);
  if (checkA == true)
    console.log("\nVerified equal ");
  else
    console.log("\nVerification FAIL ");

  let B = matTranspose(A);
  console.log("\nsource matrix B: ");
  matShow(B, 2, 5);

  let Bpinv = matPseudoInverse(B);
  console.log("\nPseudo-inverse B = ");
  matShow(Bpinv, 4, 9);

  console.log("\nVerifying B * Bpinv * B = B ");
  let BBp = matMult(B, Bpinv);
  let BBpB = matMult(BBp, B);
  matShow(BBpB, 4, 9);

  let checkB = matAreEqual(BBpB, B, 1.0e-8);
  if (checkB == true)
    console.log("\nVerified equal ");
  else
    console.log("\nVerification FAIL ");

  console.log("\nEnd demo ");
}

main();