Affinity PropagationをJavaScriptで実装した

# はじめに

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

今回はAffinity Propagationの実装について解説します。

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

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

# 概説

更新式について、いくつかパターンが存在するようです。

1.  2006年提案手法

    [別の方の投稿](https://qiita.com/daiki_yosky/items/98ce56e37623c369cc60)に紹介された更新式です。

    この元となった論文は「Frey et al. "Mixture Modeling by Affinity Propagation"(2006)」になるかと思います。
    最初はこちらを実装した（`__fit`）のですが、上手く動きませんでした。

2.  2007年提案手法

    調べると、同じ著者による別の手法が存在しました。
    論文は「Frey et al. "Clustering by Passing Messages Between Data Points"(2007)」になります。

    「Fujiwara et al. "Fast Algorithm for Affinity Propagation"(2011)」に同様のアルゴリズムの詳細が載っていましたので、この数式を参考に実装しました（`fit`）。

    なお、「松下ら "メッセージ集約に基づくAffinity Propagationの高速化"(2019)」にもほぼ同様のアルゴリズムが載っていました。英語が苦手な方はこちらも参考にしてください。

## アルゴリズム

「Fujiwara et al. "Fast Algorithm for Affinity Propagation"(2011)」からの引用となります。
計算アルゴリズム部分のみですので、理論が知りたい方は論文を参照してください。

### 初期化

データ$x_i (i = 1, 2, ..., N)$に対して、零行列$\bf\it R, A \in \mathbb{R}^{N \times N}$を用意します。

また行列$\bf\it S \in \mathbb{R}^{N \times N}$は、適当な類似度を計算する関数$\it{D} \left(x_i, x_j \right)$を用いて

```math
s_{i, j} = \it{D} \left(x_i, x_j \right)
```
と設定します。
実装では、$\it D$はユークリッド距離の二乗に$-1$を掛けた値を使用しています。ただし、対角成分に関しては非対角成分の最小値としています。

また$0 < \lambda < 1$となる変数を用意します。

### 更新

```math
\begin{eqnarray}
\rho_{i, j} &=& \left\{ \begin{array}{ll}
  s_{i, j} - \max_{k \not= j} \left( a_{i, k} + s_{i, k} \right) & (i \not= j) \\
  s_{i, j} - \max_{k \not= j} s_{i, k} & (i = j)
\end{array} \right.
\\
\alpha_{i, j} &=& \left\{ \begin{array}{ll}
  \min \left(0, r_{j, j} + \sum_{k \not= i, j} \max \left(0, r_{k, j} \right) \right) & (i \not= j) \\
  \sum_{k \not= i} \max \left(0, r_{k, j} \right) & (i = j)
\end{array} \right.
\end{eqnarray}
```

を計算し、これより

```math
\begin{eqnarray}
r_{i, j} &=& (1 - \lambda) \rho_{i, j} + \lambda r_{i, j}
\\
a_{i, j} &=& (1 - \lambda) \alpha_{i, j} + \lambda a_{i, j}
\end{eqnarray}
```

と更新します。

### 推定

データ$x_i$のexemplar（模範）となるデータは、次の通り求めます

```math
\newcommand{\argmax}{\mathop{\rm arg~max}\limits}

\argmax_{j} \left( r_{i, j} + a_{i, j} \right)
```

# コード

## Affinity Propagation

```javascript
class AffinityPropagation {
	// https://qiita.com/daiki_yosky/items/98ce56e37623c369cc60
	// https://tjo.hatenablog.com/entry/2014/07/31/190218
	constructor() {
		this._epoch = 0;
		this._l = 0.8;

		this._x = [];
		this._y = null;
	}

	get centroidCategories() {
		const y = this.predict();
		return [...new Set(y)]
	}

	get centroids() {
		return this.centroidCategories.map(i => this._x[i]);
	}

	get size() {
		const y = this.predict();
		return new Set(y).size
	}

	get epoch() {
		return this._epoch
	}

	init(datas) {
		this._x = datas;
		const n = datas.length;
		this._r = Array(n);
		this._a = Array(n);
		this._ar = Array(n);
		this._s = Array(n);
		this._as = Array(n);
		for (let i = 0; i < n; i++) {
			this._r[i] = Array(n).fill(0)
			this._a[i] = Array(n).fill(0)
			this._ar[i] = Array(n).fill(0)
			this._s[i] = Array(n)
			this._as[i] = Array(n)
		}
		this._y = null;
		this._epoch = 0;

		let min = Infinity
		for (let i = 0; i < n; i++) {
			for (let j = 0; j < i; j++) {
				if (i === j) continue
				const d = -this._x[i].reduce((s, v, k) => s + (v - this._x[j][k]) ** 2, 0)
				this._s[i][j] = this._s[j][i] = d;
				this._as[i][j] = this._as[j][i] = d;
				min = Math.min(min, d)
			}
		}
		for (let i = 0; i < n; i++) {
			this._s[i][i] = this._as[i][i] = min;
		}
	}

	fit() {
		// Frey. et al. "Clustering by Passing Messages Between Data Points" (2007)
		// "Fast Algorithm for Affinity Propagation"
		const x = this._x;
		const n = x.length;
		const l = this._l

		for (let i = 0; i < n; i++) {
			for (let k = 0; k < n; k++) {
				let m = -Infinity
				const ss = (i === k) ? this._s[i] : this._as[i];
				for (let kd = 0; kd < n; kd++) {
					if (k === kd) continue
					m = Math.max(m, ss[kd])
				}
				this._r[i][k] = l * this._r[i][k] + (1 - l) * (this._s[i][k] - m)
			}
		}

		for (let i = 0; i < n; i++) {
			for (let k = 0; k < n; k++) {
				let s = (i === k) ? 0 : this._r[k][k];
				for (let id = 0; id < n; id++) {
					if (id !== i && id !== k) {
						s += Math.max(0, this._r[id][k])
					}
				}
				if (i !== k) s = Math.min(0, s)
				const aik = l * this._a[i][k] + (1 - l) * s;
				this._a[i][k] = aik;
				this._ar[i][k] = aik + this._r[i][k];
				this._as[i][k] = aik + this._s[i][k];
			}
		}

		this._y = null;
		this._epoch++;
	}

	__fit() {
		// Frey. et al. "Mixture Modeling by Affinity Propagation" (2006)
		const x = this._x;
		const n = x.length;

		for (let i = 0; i < n; i++) {
			for (let j = 0; j < n; j++) {
				let s = 0;
				for (let k = 0; k < n; k++) {
					if (k === j) continue
					s += this._a[i][k] * this._s[i][k];
				}
				this._r[i][j] = this._s[i][j] / s
			}
		}

		for (let i = 0; i < n; i++) {
			let p = 1;
			for (let k = 0; k < n; k++) {
				if (k === i) continue
				p *= (1 + this._r[k][i]);
			}
			this._a[i][i] = p - 1
			this._ar[i][i] = this._a[i][i] + this._r[i][i]

			for (let j = 0; j < n; j++) {
				if (i === j) continue;
				p = 1 / this._r[i][i] - 1;
				for (let k = 0; k < n; k++) {
					if (k === i || k === j) continue
					p *= 1 / (1 + this._r[k][i])
				}
				this._a[i][j] = 1 / (1 + p)
				this._ar[i][j] = this._a[i][j] + this._r[i][j]
			}
		}
		this._y = null;
		this._epoch++;
	}

	predict() {
		if (!this._y) {
			this._y = [];
			const n = this._x.length;
			for (let i = 0; i < n; i++) {
				let max_v = -Infinity;
				let max_i = -1;
				for (let j = 0; j < n; j++) {
					const v = this._ar[i][j];
					if (max_v < v) {
						max_v = v;
						max_i = j;
					}
				}
				this._y.push(max_i);
			}
		}
		return this._y;
	}
}
```