t-SNE-CUDAで大規模データの超速次元圧縮&可視化

# 概要
以前[t-SNEの説明記事を書きました](https://qiita.com/stfate/items/8988d01aad9596f9d586)が，[t-SNE-CUDA](https://github.com/CannyLab/tsne-cuda)なるパッケージが出ていました．
名前の通りt-SNEをCUDAで高速化するものですが，[HPML2018(High Performance Machine Learning 2018 Workshop)](https://hpml2018.github.io/)で発表されました．

- [論文](https://arxiv.org/abs/1807.11824)
- [スライド](https://hpml2018.github.io/HPML2018_6.pdf)

単にGPU実装がんばりました，だけではなくアルゴリズム上の工夫もしていて，
SGDの更新式を物理の知見を用いてAttractive ForcesとRepulsive Forcesの項に分けた上で，
以下の各ステップの計算をCUDAで高速化しているようです．

- $P_{ij}$の計算
- $P_{ij}$と$Q_{ij}$の積
- Attractive forcesの計算
- Repulsive forcesにおけるBarnes-Hut treeの構築
- Replusive forcesのtree走査
- 低次元空間上の点におけるAttractive/Repulsive forcesの適用

t-SNEはサンプル数が増えると計算時間ものすごく増えるので，高速化されるのはわりと嬉しいですが，
やはり気になるのが「**はたしてどのくらい速くなるのか？**」
ということで試してみましたのでレポートします．

# インストール
詳細は[公式リポジトリのInstallation](https://github.com/CannyLab/tsne-cuda/wiki/Installation)を参照いただければと思いますが，
Condaでのインストールが楽です．とりあえずUbuntu16.04+Anaconda3では成功しました．
ソースからのビルドは試してないですが，依存ライブラリが多く結構めんどくさそうです．

# 試してみた
とりあえずMNISTでデータ数変えて試してみました．
ただt-SNEと比較するだけだと論文と同じことしてるだけなので，[UMAP](https://github.com/lmcinnes/umap)との比較もしてみました．

```python
import scipy as sp
from sklearn.datasets import fetch_mldata
import json
from tqdm import tqdm
from tsnecuda import TSNE
import adelheid.plot as pp

import time


if __name__ == "__main__":
    COLORS = ["#000000", "#ffff33", "#66ff66", "#00ccff", "#6633ff", "#cc3399", "#ff3300", "#0066cc", "#9933cc", "#006633"]
    mnist = fetch_mldata("MNIST original", data_home="./")
    n_obs_all = len(mnist["data"])
    n_obs_list = [100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 10000, 20000, 30000, 40000, 50000, 60000, n_obs_all]

    result = {}
    for n_obs in tqdm(n_obs_list):
        mnist_data = mnist["data"][:n_obs]
        mnist_labels = mnist["target"][:n_obs].astype(int)

        st = time.clock()

        model = TSNE(n_components=2, perplexity=30.0, theta=0.5, n_iter=1000)
        predicted = model.fit_transform(mnist_data)

        ed = time.clock()
        result[n_obs] = ed-st

        xmin = predicted[:,0].min()
        xmax = predicted[:,0].max()
        ymin = predicted[:,1].min()
        ymax = predicted[:,1].max()

        output_fn = f"mnist_tsnecuda_obs{n_obs}.png"
        fig,ax = pp.subplots(figsize=(16,12))
        unique_labels = sp.unique(mnist_labels)
        for _label in unique_labels:
            inds = sp.where(mnist_labels == _label)[0]
            ax.scatter(predicted[inds,0], predicted[inds,1], c=COLORS[_label], label=_label)

        ax.set_xlim(xmin, xmax)
        ax.set_ylim(ymin, ymax)
        ax.set_xlabel("component 0")
        ax.set_ylabel("component 1")
        ax.set_title("MNIST t-SNE-CUDA visualization")
        ax.legend()
        pp.savefig(output_fn)

    json.dump(result, open("tsnecuda_results.json", "w"), ensure_ascii=False, indent=2, sort_keys=True)
```

t-SNEやUMAPも含めたテスト用コードは[こちら](https://github.com/stfate/mnist-tsnecuda-test)に置いてあります．
(上記コードを動作させるには[こちらのパッケージ](https://github.com/stfate/adelheid)が必要です)

![tsne-cuda-proctime.png](https://qiita-image-store.s3.amazonaws.com/0/36710/0bb1f85c-6676-6ebb-bb22-cfbc5b65b48b.png)


データ数が10000を超えたあたりからt-SNE-CUDAが有意に速くなりだし，MNISTフルサイズの70000点だと`174.46`倍ほど高速化されますね．
論文で報告されている数値とほぼ同一です．
UMAPと比較しても`4.47`倍ほど速いですが，UMAPはt-SNEと比べて各クラスタの中心が離れてクラスタが凝集しやすくなる利点があるので，
このくらいの速度差であれば結局UMAP使ったほうがいいかも，という気もします．

次にt-SNE/t-SNE-CUDA/UMAPで得られた2次元マップを見てみます．

## t-SNE
![mnist_tsne_obs70000.png](https://qiita-image-store.s3.amazonaws.com/0/36710/128f7ac9-939c-73c8-8ad2-23e7696a70c0.png)

9のクラスタが2つできてしまっていますが，概ね妥当なクラスタですね．

## t-SNE-CUDA
![mnist_tsnecuda_obs70000.png](https://qiita-image-store.s3.amazonaws.com/0/36710/ea1c1b9b-bfdd-b5fb-9e55-f074522f08d4.png)

t-SNEとほぼ同一の傾向が得られていますが，複数クラスタに分かれてしまっているクラスが多くなっています．

## UMAP
![mnist_umap_obs70000.png](https://qiita-image-store.s3.amazonaws.com/0/36710/0405f3b1-f3a2-1efc-45bd-55e9a1e3dfee.png)

クラスタがより凝集し，各クラスタが見やすくなっています．
やはり見やすさの観点ではt-SNEよりも上ですね．

# まとめ
t-SNE-CUDAを動かしてMNISTにてBarnes Hut t-SNE/t-SNE-CUDA/UMAPの比較を行いました．

- データ数10000を超えたあたりから有意に速くなる
- 次元圧縮の特性としてはt-SNEとほぼ同一

純粋な速度比較だけであればUMAPよりも高速ですが，低次元空間の特性やクラスタの見やすさを考えるとUMAPを使う方が総合的にはよいかもしれません．
数万点以上の規模のデータに対してt-SNEを適用する必要がある場合には有用だと思います．

おしまい

# 参考文献
- [t-SNE-CUDA: GPU-Accelerated t-SNE and its Applications to Modern Data](https://arxiv.org/abs/1807.11824)