コレスキー分解をJavaScriptで実装する

1. 背景

※ この記事では実行列のみを対象としている。

行列 $A$ が対称かつ正定値のとき、下三角行列 $L$ とその転置によって $A$ を以下のように分解できることが知られており、そのような分解をコレスキー分解（Cholesky Decomposition, Cholesky Factorization）と呼ぶ。

$$
A = L L^T
$$

Pythonならnumpy.linalg.choleskyによってコレスキー分解をすることが可能だが、JavaScriptのメジャーなライブラリには実装されていない。個人的にJavaScriptで行列計算を行いたい場合はmath.jsやTensorFlow.jsをまず検討するのだが、前者にはLU分解やQR分解、後者にはQR分解が実装されているだけである。
コレスキー分解を提供しているJavaScriptのライブラリも無くはないのだが、今回は自分で実装することにした。

2. 実装

以下のページを参考に実装した。コード上のindex等はこの記事に合わせている。
https://rosettacode.org/wiki/Cholesky_decomposition

function choleskyDecomposition(matrix) {
  // Argument "matrix" can be either math.matrix or standard 2D array
  const A = math.matrix(matrix);
  // Matrix A must be symmetric
  console.assert(math.deepEqual(A, math.transpose(A)));

  const n = A.size()[0];
  // Prepare 2D array with 0
  const L = new Array(n).fill(0).map(_ => new Array(n).fill(0));

  d3.range(n).forEach(i => {
    d3.range(i+1).forEach(k => {
      const s = d3.sum(d3.range(k).map(j => L[i][j]*L[k][j]));
      L[i][k] = i === k ? math.sqrt(A.get([k, k]) - s) : 1/L[k][k] * (A.get([i, k]) - s);
    })
  });
  return L;
}

d3.jsとmath.jsに依存している。d3の依存性を除きたければfor文等で容易に代替できるはず。

3. Usage

二次元配列またはmath.jsのmatrixを引数に与えると、二次元配列として行列 $L$ を返す。
以下のようにするとコレスキー分解の正しさを確かめることができる。

const m1 = [[25, 15, -5],
            [15, 18,  0],
            [-5,  0, 11]];
const L1 = choleskyDecomposition(m1);

math.deepEqual(m1, math.multiply(L1, math.transpose(L1)))
// => true

4. Demo

JavaScriptでインタラクティブなノートブックを作成できるObservableというサービスがあり、そこに公開している。
https://observablehq.com/@sw1227/cholesky-decomposition