階層的クラスタリングの計算速度を改善（高速化）する方法

# はじめに

[階層的クラスタリング（Hierarchical clustering）](https://en.wikipedia.org/wiki/Hierarchical_clustering)という手法を知った。
詳しい使い方や解説は、以下を参照してほしい。
- [凝集型(階層型)クラスタリングを理解する](https://qiita.com/g-k/items/8f0d9905d3e106caed59)
- [階層的クラスタリング](https://chokkan.github.io/mlnote/unsupervised/02hac.html)

具体的な使用法は`scipy`を使用した場合は以下のようになる。
```python
from scipy.cluster.hierarchy import linkage, fcluster

X = [[i] for i in [2, 8, 0, 4, 1, 9, 9, 0]]

Z = linkage(X, method='average', metric='euclidean')
clusters = fcluster(Z, t=0.15, criterion='distance')
```


# 問題
しかし、入力の数（`X`のサイズ）が大きくなってくると、非常に時間がかかるようになる。
3000次元のベクトルを10000個なんて入れた日には目も当てられない。

階層的クラスタリングの場合、計算量は$O(n^2\log n)$や$O(n^2)$[^1]となるらしい。
各要素間の距離情報が必要であり、要素の組み合わせの数だけ計算が必要になるため、非常に計算コストがかかることになる。

# 解決方法

処理に時間がかかっているのは距離行列を求めるところであるため、計算アルゴリズムの速い関数を使って、事前に距離行列を求めておき`linkage`に入力する。

```python
import numpy as np
from sklearn.metrics import pairwise_distances_chunked
from scipy.cluster.hierarchy import linkage, fcluster

gen = pairwise_distances_chunked(X, metric='euclidean')
s_tmp = np.concatenate(list(gen), axis=0)
S = s_tmp[np.triu_indices(s_tmp.shape[0], k=1)]
Z = linkage(S, method='average')
clusters = fcluster(Z, t=0.15, criterion='distance')
```
全然計算速度が違う。感動もの！
この方法は[ここ](https://stackoverflow.com/questions/11251148/parallel-construction-of-a-distance-matrix)で見つけた（metricをmethodと書き間違えていたり、カッコ（]）を書き忘れたりしてるので注意）。

## 解説

[`linkage`関数](https://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.linkage.html)は、入力したデータに対し、距離行列を求めてクラスタリングを行っているが、距離行列を入力することもできる。
`linkage`の処理は以下の様に書き直すことが出来る。

```python
Z = linkage(X, method='average', metric='euclidean')
```
↓
```python
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage

S = pdist(X, metric='euclidean')
Z = linkage(S, method='average')
```
[`pdist`関数](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.pdist.html)は、距離行列を計算する関数である。
この`pdist`があまり早くないようで、これを[`pairwise_distances`関数](https://scikit-learn.org/dev/modules/generated/sklearn.metrics.pairwise_distances.html)に変更する。
↓
```python
import numpy as np
from sklearn.metrics import pairwise_distances
from scipy.cluster.hierarchy import linkage

gen = pairwise_distances(X, metric='euclidean')
S = gen[np.triu_indices(gen.shape[0], k=1)]
Z = linkage(S, method='average')
```
`pairwise_distances`関数が早いのは、おそらく並列処理をしていると思われる。
[`np.triu_indices`関数](https://numpy.org/doc/stable/reference/generated/numpy.triu_indices.html)を使って配列を操作しているのは凝集距離行列にするためである。
`pairwise_distances`関数は、各要素間の距離を要素数×要素数の行列で出力する。これは対称行列であり、階層的クラスタリングでは、対角要素と片方は不要となる。
必要な距離のみを再配置した行列を凝集距離行列（a condensed distance matrix）と呼ぶらしい。
`np.triu_indices`関数は上三角のインデックスを取得する関数であり、距離行列の必要な距離情報を抽出するのに使われている。

ただ大規模なデータの場合、`pairwise_distances`関数がエラーを返すことがあるらしい。その問題を回避するために、チャンクごとに距離行列を計算する[`pairwise_distances_chunked`関数](https://scikit-learn.org/dev/modules/generated/sklearn.metrics.pairwise_distances_chunked.html)を用いると良いとのこと。（ただし、これは古い情報かもしれない）
↓
```python
import numpy as np
from sklearn.metrics import pairwise_distances_chunked
from scipy.cluster.hierarchy import linkage, fcluster

gen = pairwise_distances_chunked(X, metric='euclidean')
s_tmp = np.concatenate(list(gen), axis=0)
S = s_tmp[np.triu_indices(s_tmp.shape[0], k=1)]
Z = linkage(S, method='average')
```

# おわりに

扱うデータの大きさも増えてきて、「数の暴力」に圧倒される日々が続きますが、先人が切り開いたイケてるアルゴリズムとマシンパワーで乗り切りたいと思います。

[^1]: [クラスタリング (クラスター分析)](https://www.kamishima.net/jp/clustering/)