階層型クラスタリングをJavaScriptで実装した

# はじめに

「[色々な機械学習処理をブラウザ上で試せるサイトを作った](https://qiita.com/norimy/items/669c19ecf2256b3ecae5)」中で実装したモデルの解説の三回目です。

今回は階層型クラスタリングの実装について解説します。

デモは[こちら](https://ai-on-browser.github.io/)から。（TaskをClusteringにして、Modelのhierarchyを選択）
実際のコードは[hierarchy.js](https://github.com/ai-on-browser/ai-on-browser.github.io/blob/master/model/hierarchy.js)にあります。

なお、可視化部分については一切触れません。

# 概説

階層型クラスタリングのうち、以下の７手法を実装します。

- Complete Linkage（完全連結法）
- Single Linkage（単純連結法）
- Group Average（群平均法／UPGMA）
- Ward's（ウォード法）
- Centroid（重心法／UPGMC）
- Weighted Average（重み付き平均法／WPGMA）
- Median（メジアン法／WPGMC）

実装は[こちらの資料](http://www.kamishima.net/archive/clustering.pdf)を参考にしました。

手法ごとで処理が異なる部分は、クラスタ間の距離の計算とその更新式です。
[k-means法](https://qiita.com/norimy/items/c2f99d7700a0fb50f156)ではDependency Injectionの雰囲気で実装しましたが、こちらはクラスを継承して実装しました。

## 処理の流れ

### 階層構造の解析

`fit`関数で木構造を使って階層構造を一気に解析します。

処理の流れは次の通りです。

1. 木構造のルート直下に、葉としてデータが一つのみ存在するクラスタを設定する。また同時に、各データ間の距離を計算する
2. ルート直下の全てのノード間の距離を計算する
3. ルート直下のノードが一つの場合は、処理を終了する
4. ルート直下のノードの中から、最も距離の近い二つのノードを選択する。この二つのノードをまとめたものが、新しいクラスタとなる。
5. ルート直下から選択された二つのノードを削除し、その二つを子に持つ部分木をルート直下に入れる
6. ルート直下の全てのノード間の距離を再計算し、3.に戻る

データ間の距離はユークリッド距離、マンハッタン距離、チェビシェフ距離の三つを使用できるようにしています。

ノード間の距離を`distances`として二次元配列で保持し、再計算時はノードの削除・挿入と位置が合うように縮めていきます。

また、各ノードには子となる二つのクラスタ間距離を保持しておきます。

### Lance-Williamsの更新式

ノード間の距離の再計算は部分的な情報のみから計算することができ、計算量は$O(1)$になります。

結合前のクラスタを$C_a, C_b$、結合後のクラスタを$C_c$とし、それらとは別のクラスタを$C_k$とします。
また、クラスタ$C_x, C_y$間の距離を$d_{xy}$、クラスタ$C_x$に含まれるデータ数を$N_{C_x}$とします。

結合前のクラスタ間の距離$d_{ka}, d_{kb}$が分かっていたとき、結合後のクラスタ間の距離$d_{kc}$は、[こちらの資料](http://www.kamishima.net/archive/clustering.pdf)によると

```math
d_{kc} = \alpha_a d_{ka} + \alpha_b d_{kb} + \beta d_{ab} + \gamma \left| d_{ka} - d_{kb} \right|
```

となります。ただし式中の$\alpha_a, \alpha_b, \beta, \gamma$は、手法別に次の通りとなります。

| | $\alpha_a$ | $\alpha_b$ | $\beta$ | $\gamma$ |
| --- | :-: | :-: | :-: | :-: |
| Complete Linkage | $\frac{1}{2}$ | $\frac{1}{2}$ | $0$ | $\frac{1}{2}$ |
| Single Linkage | $\frac{1}{2}$ | $\frac{1}{2}$ | $0$ | $-\frac{1}{2}$ |
| Group Average | $\frac{N_{C_a}}{N_{C_c}}$ | $\frac{N_{C_b}}{N_{C_c}}$ | $0$ | $0$ |
| Ward's | $\frac{N_{C_k} + N_{C_a}}{N_{C_k} + N_{C_c}}$ | $\frac{N_{C_k} + N_{C_b}}{N_{C_k} + N_{C_c}}$ | $-\frac{N_{C_k}}{N_{C_k} + N_{C_c}}$ | $0$ |
| Centroid | $\frac{N_{C_a}}{N_{C_c}}$ | $\frac{N_{C_b}}{N_{C_c}}$ | $-\frac{N_{C_a}N_{C_b}}{N_{C_c}^2}$ | $0$ |
| Weighted Average | $\frac{1}{2}$ | $\frac{1}{2}$ | $0$ | $0$ |
| Median | $\frac{1}{2}$ | $\frac{1}{2}$ | $-\frac{1}{4}$ | $0$ |

### クラスタの取得

上記の流れで得られた木を使用すると、クラスタの取得はシンプルになります。
木の中のクラスタ間距離が大きいノードで分割していき、指定されたクラスタ数になるまで続けるだけです。

クラスタ間距離は親ノードよりも子ノードの方が小さいことが保証されますので、幅優先探索の要領でルートから順番にたどっていけば、不要な探索が無くなります。

処理の流れは次の通りです。

1. クラスタ一覧にすべてのデータが含まれたノードを設定する
2. クラスタ一覧に含まれるノード数が指定されたクラスタ数になった場合は、そのクラスタ一覧を返して処理を終了する
3. クラスタ一覧の中から、分割したのちのクラスタ間距離が一番大きいノードを選択する
4. 選択したノードをクラスタ一覧から取り除き、その子ノードを取り出してクラスタ一覧に追加する
5. 2.に戻る

なお木の配列を返しているのは、可視化時に使用するためです。
実際のデータは葉ノードにありますので、その部分を走査すれば取得できます。

## 木

木構造は、手作り数学ライブラリ[math.js](https://github.com/ai-on-browser/ai-on-browser.github.io/blob/master/js/math.js)に定義した`Tree`を使用します。

ここで使用している処理は次の通りです。

- プロパティ
    - length : 子ノードの数
    - value : ノードが持っている値
- 関数
    - at : 指定された位置の子ノードを返す
    - push : 子ノードの最後に渡された値を持った木を設定する。ただし、木が渡された場合はそのまま設定する
    - set : 指定した位置に渡された値を持った木を設定する。ただし、木が渡された場合はそのまま設定する
    - removeAt : 指定された位置の子ノードを削除する
    - leafCount : 葉の数を返す
    - isLeaf : 葉ノードの場合は`true`を返す
    - leafValues : 葉ノードの値を配列として返す

# コード

## 親クラス

子クラスで実装する関数は`distance`と`update`で、それぞれクラスタ間距離の計算とクラスタ間距離の更新計算を行います。

また、`_lanceWilliamsUpdater`によりLance-Williamsの更新式を計算する関数を返します。

```javascript
class HierarchyClustering {
	constructor(metric = 'euclid') {
		this._root = null;
		this._metric = metric;

		switch (this._metric) {
		case 'euclid':
			this._d = (a, b) => Math.sqrt(a.reduce((s, v, i) => s + (v - b[i]) ** 2, 0));
			break
		case 'manhattan':
			this._d = (a, b) => a.reduce((s, v, i) => s + Math.abs(v - b[i]), 0)
			break
		case 'chebyshev':
			this._d = (a, b) => a.reduce((s, v, i) => Math.max(s, Math.abs(v - b[i])), -Infinity)
			break;
		}
	}

	fit(points) {
		this._root = new Tree()
		points.forEach((v, i) => {
			this._root.push({
				point: v,
				index: i,
				distances: points.map(p => this._d(v, p))
			});
		});

		const distances = []
		for (let i = 0; i < this._root.length; i++) {
			if (!distances[i]) distances[i] = [];
			for (let j = 0; j < i; j++) {
				if (!distances[i][j]) distances[i][j] = distances[j][i] = this.distance(this._root.at(i), this._root.at(j));
			}
		}
		while (this._root.length > 1) {
			let n = this._root.length;

			let min_i = 0;
			let min_j = 1;
			let min_d = distances[0][1];
			for (let i = 1; i < n; i++) {
				distances[i].forEach((d, j) => {
					if (d < min_d) {
						min_i = i;
						min_j = j;
						min_d = d;
					}
				});
			}
			let min_i_leafs = this._root.at(min_i).leafCount();
			let min_j_leafs = this._root.at(min_j).leafCount();
			distances.forEach((dr, k) => {
				if (k != min_j && k != min_i) {
					dr[min_i] = this.update(min_i_leafs, min_j_leafs, this._root.at(k).leafCount(), dr[min_i], dr[min_j], distances[min_j][min_i]);
					distances[min_i][k] = dr[min_i];
					dr.splice(min_j, 1);
				}
			});
			distances[min_i].splice(min_j, 1);
			distances.splice(min_j, 1);
			this._root.set(min_i, new Tree({
				distance: min_d,
			}, [this._root.at(min_i), this._root.at(min_j)]));
			this._root.removeAt(min_j);
		}
		this._root = this._root.at(0);
	}

	getClusters(number) {
		const scanNodes = [this._root]
		while (scanNodes.length < number) {
			let max_distance = 0;
			let max_distance_idx = -1;
			for (let i = 0; i < scanNodes.length; i++) {
				const node = scanNodes[i];
				if (!node.isLeaf() && node.value.distance > max_distance) {
					max_distance_idx = i;
					max_distance = node.value.distance
				}
			}
			if (max_distance_idx === -1) {
				break
			}
			const max_distance_node = scanNodes[max_distance_idx];
			scanNodes.splice(max_distance_idx, 1, max_distance_node.at(0), max_distance_node.at(1))
		}
		return scanNodes;
	}

	distance(c1, c2) {
		throw new Error('Not Implemented');
	}

	_mean(d) {
		const m = Array(d[0].length).fill(0);
		for (let i = 0; i < d.length; i++) {
			for (let k = 0; k < d[i].length; k++) {
				m[k] += d[i][k]
			}
		}
		return m.map(v => v / d.length);
	}

	_lanceWilliamsUpdater(ala, alb, bt, gm) {
		return (ka, kb, ab) => ala * ka + alb * kb + bt * ab + gm * Math.abs(ka - kb);
	}

	update(ca, cb, ck, ka, kb, ab) {
		throw new Error('Not Implemented');
	}
}
```

## Complete Linkage（完全連結法）

クラスタ間距離は、それぞれのクラスタに含まれるデータの中で一番遠いデータ同士の距離になります。

```JavaScript
class CompleteLinkageHierarchyClustering extends HierarchyClustering {
	distance(c1, c2) {
		let f1 = c1.leafValues();
		let f2 = c2.leafValues();
		return Math.max.apply(null, f1.map(v1 => {
			return Math.max.apply(null, f2.map(v2 => v1.distances[v2.index]));
		}));
	}

	update(ca, cb, ck, ka, kb, ab) {
		return this._lanceWilliamsUpdater(0.5, 0.5, 0, 0.5)(ka, kb, ab)
	}
}
```

## Single Linkage（単純連結法）

クラスタ間距離は、それぞれのクラスタに含まれるデータの中で一番近いデータ同士の距離になります。

```JavaScript
class SingleLinkageHierarchyClustering extends HierarchyClustering {
	distance(c1, c2) {
		let f1 = c1.leafValues();
		let f2 = c2.leafValues();
		let minDistance = Math.min.apply(null, f1.map(v1 => {
			return Math.min.apply(null, f2.map(v2 => v1.distances[v2.index]));
		}));
		return minDistance;
	}

	update(ca, cb, ck, ka, kb, ab) {
		return this._lanceWilliamsUpdater(0.5, 0.5, 0, -0.5)(ka, kb, ab)
	}
}
```

## Group Average（群平均法／UPGMA）

クラスタ間距離は、それぞれのクラスタに含まれるデータの間の距離の平均になります。

```JavaScript
class GroupAverageHierarchyClustering extends HierarchyClustering {
	distance(c1, c2) {
		let f1 = c1.leafValues();
		let f2 = c2.leafValues();
		let totalDistance = f1.reduce((acc1, v1) => {
			return acc1 + f2.reduce((acc2, v2) => acc2 + v1.distances[v2.index], 0);
		}, 0);
		return totalDistance / (f1.length * f2.length);
	}

	update(ca, cb, ck, ka, kb, ab) {
		return this._lanceWilliamsUpdater(ca / (ca + cb), cb / (ca + cb), 0, 0)(ka, kb, ab);
	}
}
```

## Ward's（ウォード法）

クラスタ間距離は、結合後のクラスタの重心との距離の平均から、結合前それぞれのクラスタの重心との距離の平均を引いた値になります。

```JavaScript
class WardsHierarchyClustering extends HierarchyClustering {
	distance(c1, c2) {
		let f1 = c1.leafValues().map(f => f.point);
		let f2 = c2.leafValues().map(f => f.point);
		let fs = f1.concat(f2);
		let ave1 = this._mean(f1);
		let ave2 = this._mean(f2);
		let aves = this._mean(fs);
		let e1 = f1.reduce((acc, f) => acc + this._d(f, ave1) ** 2, 0);
		let e2 = f2.reduce((acc, f) => acc + this._d(f, ave2) ** 2, 0);
		let es = fs.reduce((acc, f) => acc + this._d(f, aves) ** 2, 0);
		return es - e1 - e2;
	}

	update(ca, cb, ck, ka, kb, ab) {
		return this._lanceWilliamsUpdater((ck + ca) / (ck + ca + cb), (ck + cb) / (ck + ca + cb), -ck / (ck + ca + cb), 0)(ka, kb, ab);
	}
}
```

## Centroid（重心法／UPGMC）

クラスタ間距離は、それぞれのクラスタの重心の距離になります。

```JavaScript
class CentroidHierarchyClustering extends HierarchyClustering {
	distance(c1, c2) {
		let f1 = c1.leafValues().map(f => f.point);
		let f2 = c2.leafValues().map(f => f.point);
		let d = this._d(this._mean(f1), this._mean(f2));
		return d * d;
	}

	update(ca, cb, ck, ka, kb, ab) {
		return this._lanceWilliamsUpdater(ca / (ca + cb), cb / (ca + cb), -ca * cb / ((ca + cb) * (ca + cb)), 0)(ka, kb, ab);
	}
}
```

## Weighted Average（重み付き平均法／WPGMA）

クラスタ間距離は、二つのクラスタ間の距離の平均を再帰的に計算していくことにより求めます。

```JavaScript
class WeightedAverageHierarchyClustering extends HierarchyClustering {
	distance(c1, c2) {
		let calcDistRec = function calcDistRec(h1, h2) {
			if (h1.leafCount() == 1 && h2.leafCount() == 1) {
				return h1.value.distances[h2.value.index];
			} else if (h2.leafCount() == 1) {
				return (calcDistRec(h2, h1.at(0)) + calcDistRec(h2, h1.at(1))) / 2;
			} else {
				return (calcDistRec(h1, h2.at(0)) + calcDistRec(h1, h2.at(1))) / 2;
			}
		}
		return calcDistRec(c1, c2);
	}

	update(ca, cb, ck, ka, kb, ab) {
		return this._lanceWilliamsUpdater(0.5, 0.5, 0, 0)(ka, kb, ab);
	}
}
```

## Median（メジアン法／WPGMC）

クラスタ間距離は、それぞれのクラスタの重心の中間点と、一方の重心との距離になります。

```JavaScript
class MedianHierarchyClustering extends HierarchyClustering {
	distance(c1, c2) {
		let m1 = this._mean(c1.leafValues().map(f => f.point));
		let m2 = this._mean(c2.leafValues().map(f => f.point));
		let m = m1.map((v, i) => (v + m2[i]) / 2);
		return this._d(m, m2) ** 2;
	}

	update(ca, cb, ck, ka, kb, ab) {
		return this._lanceWilliamsUpdater(0.5, 0.5, -0.25, 0)(ka, kb, ab);
	}
}
```

# さいごに

文末にセミコロンが付いたり付かなかったり、`const`でいいところを`let`を使用していたりと、気になる部分が多い。