Cholesky decomposition computes a matrix L from a source matrix A such that L * trans(L) = A. The L matrix is a lower triangular (zeros above the diagonal).

Cholesky decomposition applies only to special matrices that are square, symmetric and positive definite (PD). In machine learning, such special matrices occur during Gaussian process regression and Gaussian mixture model clustering. You can find the inverse of a PD matrix relatively easily, as opposed to finding the inverse of an arbitrary matrix, typically using Crout’s decomposition, which is significantly more difficult.

I spend 30 minutes each morning before work to explore a topic and write some code. So, one bright summer morning, I implemented Cholesky matrix inverse using JavaScript. Why JavaScript? Mostly to keep my JavaScript programming skills sharp. The key calling code is:

  let mat = matMake(4, 4, 0.0);  // set up a PD matrix
  mat[0] = [ 1.0, 0.2, 0.3, 0.4 ];
  mat[1] = [ 0.2, 1.0, 0.5, 0.6 ];
  mat[2] = [ 0.3, 0.5, 1.0, 0.7 ];
  mat[3] = [ 0.4, 0.6, 0.7, 1.0 ];

  let lower = matCholeskyDecomp(mat);  // decompose
  let inverse = matCholeskyInverse(lower);  // compute inverse
  let product = matProduct(mat, inverse);  // verify
  matShow(product2, 4, 11);  // should be Identity

It would be possible to combine the matCholeskyDecomp() and matCholeskyInverse() into a single function, but Cholesky decomposition can be used for other purposes, notably solving a system of linear equations.

There’s no big moral to this story. The minor moral is that all skills, including coding, must be practiced.

In a sense, a matrix inverse is a kind of an opposite matrix. When I was a young man, to a large extent I first learned to read from comic books. A fairly common theme was a hero who appeared to be evil instead of good — but of course it always turned out that it was a ploy to trick the villain. Left: Action Comics #311 from April 1964. Center: Batman #123 from April 1959. Right: The Flash #139 from September 1963.

Demo code. Replace “lt”, “gt”, “lte”, “gte” with Boolean operator symbols (my lame blog editor often chokes on symbols).

// test_cholesky_inverse.js

// let U = require("./Utils/utilities_lib.js");
// could put functions in a library

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

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 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 matProduct(ma, mb)
{
  let aRows = ma.length;
  let aCols = ma[0].length;
  let bRows = mb.length;
  let bCols = mb[0].length;
  if (aCols != bRows) {
    throw "Non-conformable matrices";
  }

  let result = matMake(aRows, bCols, 0.0);

  for (let i = 0; i "lt" aRows; ++i) { // each row of A
    for (let j = 0; j "lt" bCols; ++j) { // each col of B
      for (let k = 0; k "lt" aCols; ++k) { // could use bRows
        result[i][j] += ma[i][k] * mb[k][j];
      }
    }
  }

  return result;
}

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

function matTranspose(m)
{
  let r = m.length;
  let c = m[0].length;
  
  let result = matMake(c, r, 0.0);
  for (let i = 0; i "lt" r; ++i)
    for (let j = 0; j "lt" c; ++j)
      result[j][i] = 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 matCholeskyDecomp(m)
{
  // Cholesky decomposition
  // m is square, symmetric, positive definite
  let n = m.length;
  let result = matMake(n, n, 0.0);
  for (let i = 0; i "lt" n; ++i) {
    for (let j = 0; j "lte" i; ++j) {
      let sum = 0.0;
      for (let k = 0; k "lt" j; ++k)
        sum += result[i][k] * result[j][k];
      if (i == j) {
        let tmp = m[i][i] - sum;
        if (tmp "lt" 0.0)
          throw "matCholeskyDecomp fatal error ";
        result[i][j] = Math.sqrt(tmp);
      }
      else {
        if (result[j][j] == 0.0)
          throw "matCholeskyDecomp fatal error ";
        result[i][j] =
          (1.0 / result[j][j] * (m[i][j] - sum));
      }
    } // j
  } // i

  return result;
}

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

function matCholeskyInverse(L)
{
  // L is the Cholesky decomp of source matrix
  let n = L.length;
  let result = matIdentity(n);

  for (let k = 0; k "lt" n; ++k) {
    for (let j = 0; j "lt" n; j++) {
      for (let i = 0; i "lt" k; i++) {
        result[k][j] -= result[i][j] * L[k][i];
      }
      result[k][j] /= L[k][k];
    }
  }

  for (let k = n - 1; k "gte" 0; --k) {
    for (let j = 0; j "lt" n; j++) {
      for (let i = k + 1; i "lt" n; i++) {
        result[k][j] -= result[i][j] * L[i][k];
      }
      result[k][j] /= L[k][k];
    }
  }
  return result;
}

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

function main()
{
  console.log("\nBegin Cholesky using JavaScript ");

  let mat = matMake(4, 4, 0.0);
  mat[0] = [ 1.0, 0.2, 0.3, 0.4 ];
  mat[1] = [ 0.2, 1.0, 0.5, 0.6 ];
  mat[2] = [ 0.3, 0.5, 1.0, 0.7 ];
  mat[3] = [ 0.4, 0.6, 0.7, 1.0 ];

  console.log("\nSource matrix: ");
  matShow(mat, 1, 4);

  let lower = matCholeskyDecomp(mat);
  console.log("\nCholesky decomp L: ");
  matShow(lower, 4, 11);

  let lowerT = matTranspose(lower);
  console.log("\nL transpose: ");
  matShow(lowerT, 4, 11);

  let product1 = matProduct(lower, lowerT);
  console.log("\nL * Lt should be source: ");
  matShow(product1, 4, 11);

  let inverse = matCholeskyInverse(lower);
  console.log("\nCholesky inverse: ");
  matShow(inverse, 4, 11); 

  let product2 = matProduct(mat, inverse);
  console.log("\nsource * inverse should be Identity: ");
  matShow(product2, 4, 11); 

  console.log("\nEnd Cholesky using JavaScript ");
}

main();