はじめに
距離学習および深層距離学習について述べた記事の最後に，
では距離関数自体を学習可能にする枠組みは存在するのか？→Yes.
とあった．その枠組みに関する論文が，"Deep Divergence Learning" と題して ICML2020にて発表された．
Cilingir, H.K., Manzelli, R. and Kulis, B., "Deep Divergence Learning," Proceedings of the 37th International Conference on Machine Learning, 2020.
著者による実験のソースコードが長らく公開されておらず，空のリポジトリのみが存在していたが，2022年の2月にjupyter notebookが公開された．
残念ながら，このnotebookは論文に掲載された実験を再現するには余りにも不十分だった．それならば，最初から全て PyTorch & scikit-learn で書いて実験を再現してやろうと決めた（再現できるとは言ってない）．本記事はそれら実験結果のうち特に分布クラスタリングに関して再現実験を試みたときの報告である． 当初は分布クラスタリングで単独の記事のつもりとしていたが，実験単位で記事を2つに分けることにした．本記事はそのうちの最初の記事である．
論文の概要
本記事に関連する範囲で論文の概要（という名のメモ/抜粋そして私的な解釈）を述べる．
図1に提案手法と既存の距離学習手法との関係図を示す．図は上記 "Deep Divergence Learning"から引用する． 従来の距離学習における距離関数はユークリッド距離に留まり，また特徴量変換は線形なものに留まっていた（左上）．特徴量変換をニューラルネットによる非線形なものに拡張すると deep metric learning となる（左下）． 一方，距離関数自体をBregman divergenceに（学習可能な形で）拡張すると Bregman divergence learning となる（右上）． 提案手法は素のBregman divergenceではなく，functional Bregman divergenceを活用したdeep divergence learningとなる（右下）．
Bregman divergenceは2点（ベクトル）の距離を測る道具であり，凸関数$\phi$で特徴づけられる．定義を書いておこう． $\phi$をある凸集合$\Omega$上で定義された連続的微分可能な狭義凸関数とするとき，2点$\mathbf{x}, \mathbf{y} \in \Omega$に関するBregman divergenceは次式で与えられる．
$$ D_\phi(\mathbf{x}, \mathbf{y}) := \phi(\mathbf{x})-\phi(\mathbf{y})-\langle \mathbf{x} - \mathbf{y}, \nabla \phi (\mathbf{y}) \rangle. $$
機械学習や信号処理の分野でもよく出会うユークリッド距離，KL-divergence，板倉−斎藤距離などは凸関数のとり方をそれぞれ変えて得られる*1 ．
functional Bregman divergenceは関数ペアの間のdivergenceを測る道具である．こちらも定義を書いておこう． 2つの関数ペア$p$および$q$，そして狭義凸汎関数$\phi$を考える．凸汎関数は関数の凸集合で定義されており，実数に値を取るものとする．functional Bregman divergenceは次式で定義される．
$$ D_\phi(p, q) := \phi(p) - \phi(q) - \int [p(x) - q(x)] \delta \phi (q) (x) dx $$
ここで$\delta \phi (q)$は$\phi$の汎関数微分（Fréchet 微分）である．本論文では特に分布関数のペアを考える．
論文ではfunctional Bregman divergenceが対称となる必要十分条件を与えている．ここでfunctional Bregman divergenceが対称とは，（その定義域に入る）全ての関数$p, q$に対して$D_{\phi}(p, q) = D_{\phi}(q, p)$が成り立つものとして定義する．
定理 functional Bregman divergence $D_{\phi}(p, q) $が対称となるためには，対称な半正定値関数$\psi (x, y)$を用いて以下のように書けることが必要十分である． $$ D_{\phi}(p, q) = \iint (p(x) - q(x))(p(y) - q(y)) \psi(x, y) dx dy $$
そこから従来の（深層）距離学習に現れる距離関数たちが特殊な具体例として導出される．Linear Metric Learning，Kernel Metric Learningの距離関数，そして深層距離学習との関連でmoment-matching 型の距離関数が挙げられている．
$$ D_{\phi}(p, q) = \left\lVert \mathbb{E}_{p} [f_{W}] - \mathbb{E}_{q} [f_{W}] \right\rVert ^{2} $$
ただし$p$と$q$はそれぞれ分布関数，$\phi$は凸汎関数である．$f_{W}$は$W$をパラメタに持つニューラルネットによる埋め込みを表す．これは上記定理で$\psi (x, y) = f_{W}(x) ^{T} f_{W}(y)$と置いたケースである．特に$p$と$q$がそれぞれ観測データ集合$P$，$Q$に関する経験分布だとすると，次式で書ける．
$$ D_{\phi}(p, q) = \left\lVert \frac{1}{|P|} \sum_{x \in P} f_{W} (x) - \frac{1}{|Q|} \sum_{y \in Q} f_{W} (y) \right\rVert ^{2} $$
ただし$|P|$はその集合$P$の濃度を表す． $|P|=|Q|=1$の場合，経験分布はデルタ関数となり，従来のユークリッド距離に基づく深層距離学習に帰着する．上式は一次モーメントの意味でマッチする分布の距離を測っている．
さらに論文は学習可能な凸汎関数を導入したdeep Bregman divergenceの定式化へと繋がっていくが，本記事は上記moment-matchingまでの範囲で実験を行う． deep Bregman divergenceの実装と再現実験は次回以降の記事に回す．というのも，deep Bregman divergenceはクラスタリングの実験には使われていないからである．
応用可能性 −分布クラスタリング−
通常のクラスタリングはデータ点の集合に対して行われるが，ここでは経験分布（関数）の集合に対するクラスタリングを考える．経験分布はまたデータ点集合を伴う．例えばオンラインショップにおける各商品の評価値（レーティング）の分布が挙げられる; Amazonの各商品はユーザ投票による星1つ〜星5つの5段階評価値の集合と，そのヒストグラムすなわち経験分布を有している．2つの経験分布の距離を測る道具を我々は既に手にしているわけなので，その距離に基づくクラスタリングアルゴリズム，例えば$\mathit{k}$-meansクラスタリングを駆動できるだろう，と考えるのは自然な発想である．
Davis and Dhillon (2006) *2 は教師なし分布クラスタリングアルゴリズムを提案している．ただし各分布が多次元ガウス分布という仮定と，分布の平均ベクトル群が各クラスタごとに線形分離可能という暗黙の仮定を置いていると著者らは指摘する．後者の仮定は，分布群の配置があまり複雑に入り組んでおらず，「クラスタリングしやすい」状況を指摘している．その一方でdeep divergence learningの本論文において著者らは，分布クラスタリングに関して「半教師あり」（semi-supervised）な枠組みを考察している．同じクラスタに属するべき分布ペア，もしくは別のクラスタに分けるべき分布ペアが訓練データ上で与えられるという想定のもとに，contrastive loss あるいは triplet loss に基づく深層距離学習の適用を見据えている※． たとえ元の分布群の配置が線形分離不可能であっても，埋め込みを実現するニューラルネット$f_{W}$によって，埋め込み空間上で線形分離可能な配置へと分布群を移動できると期待する．あとはfunctional Bregman divergenceに基づくクラスタリングアルゴリズムを駆動すればよい．加えて各分布が多次元ガウスである仮定を置く必要もない．"Davis and Dhillon" の手法とは，これらの点で区別される．
※ 著者らの「半教師あり」はむしろ「弱教師あり」（weakly supervised）の意味に解釈すべきだろう．contrastive loss あるいは triplet lossを効率よく計算しようとするならば，ラベルの所与を前提に置かねばならない．件の分布ペアを訓練時のミニバッチ内で都度作り出す（on-the-fly）場合は特にそうである．ただしラベルとは分布相互の近隣・遠隔に関する何らかの情報を指す．本記事の実験ではクラスタラベルをそのまま利用している．いずれにせよ，各分布をあたかも個々の観測データと見なすという点を除けば，従来の深層距離学習と同じ枠組みに帰着できる．
実験： Clustering Multivariate Gaussian Distributions
実験条件
データセット生成
分布クラスタリングのデモンストレーション実験のために，著者らは多変量ガウス分布の群を取り上げた．以下，図2を参照しながらデータの作り方を説明しよう．
半径を変えた3つの同心円の円周付近にガウス分布の平均ベクトルを500個発生させる（訓練用; テスト用はこれとは別に200個）．同心円の半径はそれぞれ0.2, 0.6, 1.0であり，3つの同心円は3つのクラスタにそれぞれ対応する．1つの平均ベクトルは円周上の点をランダムにサンプリングし，2次元のガウスノイズの重畳により得られる．このときのノイズの大きさは論文に明記されていなかったが，本実験では各次元で標準偏差0.05とした．以上が図2の一番左に対応する．色付けにより，対応するクラスタを見やすくしている．
続いて各ガウス分布から複数個の点をランダムサンプリングする; 各ガウス分布の共分散の大きさについては幸いにも論文に記述があり，対角共分散行列を仮定し，対角成分はそれぞれ0.1とした．$\sqrt{0.1}$でスケーリングした2次元ガウスノイズを独立に複数個発生させ，各平均ベクトルに加えているわけである．論文には具体的にいくつ点をサンプリングするのかについて明確な言及がなかったが，論文のFigure 1には各ガウス分布から50個の点をサンプリングしたときの図が掲載されていた点を鑑みて，本実験でも1つのガウス分布あたり50個の点をサンプリングした． 図2の真ん中3つの散布図はそれぞれ，各クラスタのガウス分布からランダムサンプリングした点群を示している．一番右はこれら3つの散布図を全て重ねた図である．平均ベクトルのサンプリングと，対応するガウス分布からのサンプリングという2段階の構造になっているので気をつけて欲しい．
最終的に，手元には500 × 50 = 25,000個からなる訓練用の点群と200 × 50 = 10,000個からなるテスト用の点群があり，各点は3つのクラスタのいずれかに属する（＝ラベルが付与されている），という設定である．図2を見て分かる通り，分散が大きいため各クラスタの点群は互いに大きく重なっている．
比較手法
Davis and Dhillon (2006)：分布クラスタリングに関する従来手法
moment-matching型のfunctional Bregman divergenceに基づく深層距離学習 & $\mathit{k}$-means アルゴリズム
の3つを比較する．従来の深層距離学習では，個々のサンプル点は独立な観測データとして扱われ，クラスタリングもまた観測データ単位で行われる．"Davis and Dhillon"および提案手法の枠組みでは経験分布をクラスタリングの対象にするために，1つのガウス分布から発生させた50個のサンプル点が1つの経験分布を構成するよう，各点と経験分布とを一意に紐づけて管理する（下記補足参照）．各経験分布のクラスタラベルは所与とする．
"Davis and Dhillon"では経験分布から計算される標本平均ベクトルと標本共分散行列に基づいてアルゴリズムを駆動する．教師なしのクラスタリングアルゴリズムなので，訓練時にクラスタラベルは使用されない．評価用のreference（正解ラベル）として利用する．つまり分布クラスタリング後に，経験分布に対してクラスタラベルが予測されるので，referenceとの比較により性能を評価する．
深層距離学習に関する2手法では所与のクラスタラベルを活用してサンプル点ペアないし分布ペアを作成する．訓練が完了した後に訓練済みのニューラルネットによって各データの埋め込みを計算する（非線形変換）．そして変換後の埋め込み空間上で$\mathit{k}$-meansアルゴリズムを駆動する．学習済みの「距離関数」に基づいた教師なしクラスタリングというわけである．予測されるクラスタラベルは，従来の深層距離学習では各点（観測データ）単位で算出されるが，提案手法では経験分布単位で算出されるので，後述する評価指標の計算時には注意が必要である．
ネットワークアーキテクチャ および 損失関数
埋め込みネットワークは入力層に続いて隠れ層が2層続く構成である．ユニット数は入力層から最終層に向かって順に(2)-(1000)-(500)-(2)とした．最終層のユニット数が2となっているのは埋め込み空間の可視化を容易にするためである．層の間の活性化関数はReLUを用いた．
損失関数はcontrastive lossおよびtriplet lossである．
評価指標
クラスタリングの性能評価のために以下の2指標を用いた．参考まで，scikit-learnへのページをリンクしておく．
これらはいずれも，同じデータセットに対する2つのクラスタリング結果の間の類似度を計算するものである．ARIは特にクラスタラベルのランダム割り当てにかかるバイアスを補正した（adjusted）RIである．実験では予測クラスタラベル（predicted）と真のクラスタラベル（ground truth, reference）の間の類似度を計算する．
その他の実験条件
| 項目 | 値 |
|---|---|
| エポック数 | 50 |
| 最適化アルゴリズム | Adam |
| 学習率 | 1e-4 |
| バッチサイズ | 32 |
| triplet loss のマージン | 0.2 |
| contrastive loss のマージン | 0.2 |
補足：経験分布の構成
経験分布の構成を簡単に式に書いておこう（本記事オリジナル）．経験分布の集合$E$を
$$ E = \left\lbrace e_{1}, e_{2}, \ldots, e_{N} \right\rbrace $$
と書いてみる．$N$は経験分布の総数であり，訓練用は$N=500$である．テスト用には$N=200$として別に作成しておく．各経験分布が50個のサンプル点を持っている． つまり$e_{i}\; (1 \leq i \leq N)$ には
$$ \left\lbrace \mathbf{x}^{(i)}_{1}, \mathbf{x}^{(i)}_{2}, \ldots, \mathbf{x}^{(i)}_{50} \right\rbrace , \;\;\; \mathbf{x}^{(i)}_{j} \sim \mathcal{N} \left( \mathbf{m}^{(i)}; \begin{bmatrix} 0.1 & 0 \\ 0 & 0.1 \end{bmatrix} \right) \;\;\; (1\leq j \leq 50) $$
という点集合が付随している．ただし$\mathbf{m}^{(i)} \in \mathbb{R}^{2}$は平均ベクトルである．これらを用いて
$$ e_{i} (\mathbf{x}) = \frac{1}{50} \sum_{j=1}^{50} \delta \left(\mathbf{x} - \mathbf{x}^{(i)}_{j} \right),\; \; \; \mathbf{x} \in \mathbb{R}^{2} $$
書けるものとする．$\delta$はDiracのデルタ関数である．各経験分布はクラスタのラベルを有している．便宜上，3クラスタのラベルを$0, 1, 2$として表せば，$E$からラベル集合$C = \left\lbrace 0, 1, 2 \right\rbrace$への写像$g$が定義されている．
$$ g: E \to C; \;\; e_{i} \mapsto c \in C $$
それゆえにまた，各サンプル点にも，$e_{i}$を介してラベルが付与される．
実験結果および考察
実験結果
評価指標ごとに表を分けて示す．baselineは従来の深層距離学習，proposedは提案手法である． 平均および標準偏差は乱数の種（seed）を10回変えて得られたスコアから計算した． なおseedの値が関わるのはデータセット生成の部分のみである（経験分布を生成する際のサンプリング）． もしseedを固定すればたとえGPUを使用しても常に同じ結果が得られる．
- RI
| 手法（損失関数） | 平均 | 標準偏差 |
|---|---|---|
| Davis & Dhillon | 0.649826 | 0.008840 |
| baseline (constrative loss) | 0.650168 | 0.013206 |
| proposed (constrative loss) | 0.990864 | 0.007971 |
| baseline (triplet loss) | 0.657529 | 0.012002 |
| proposed (triplet loss) | 0.992834 | 0.006162 |
- ARI
| 手法（損失関数） | 平均 | 標準偏差 |
|---|---|---|
| Davis & Dhillon | 0.241574 | 0.009874 |
| baseline (constrative loss) | 0.238745 | 0.030552 |
| proposed (constrative loss) | 0.979434 | 0.017977 |
| baseline (triplet loss) | 0.238126 | 0.028162 |
| proposed (triplet loss) | 0.983881 | 0.013887 |
どちらの評価指標においても提案手法は"Davis & Dhillon"とbaselineを大きく上回った結果を確認できる．以上より，論文の実験結果をほぼ再現できたが，"Davis & Dhillon"の手法は今回の実験でより高いスコアを出した．
考察
すでに述べた通り，経験分布に属する各サンプル点の埋め込み平均によって，moment-matching型のfunctional Bregman divergenceを計算しながら深層距離学習を遂行するというのが提案手法の枠組みだった．今回のデータセットに関していえば，あたかも図2の一番左にプロットした各分布の平均ベクトルたちを対象に深層距離学習を行っているようなものである．平均ベクトルたちはクラスタ間で重なるようには生成されていないわけで，このデータセット上では提案手法が上手く動いたと説明できる．また経験分布単位のクラスタリングという観点で"Davis and Dhillon"と比較した場合，今回のような線形分離不可能な経験分布の初期配置であっても，クラスタリングがうまく働くように非線形変換が学習できており，提案手法は有利であるといえる．
baselineは25,000個全ての点を対象にクラスタリングしていたわけだが，もし，baselineについても埋め込み後に経験分布単位で平均ベクトルを計算し，得られた500個の点を対象にクラスタリングしたらどうなるだろうか？ 独自に実装したからこそ，このような実験も新たにできるわけである． その結果は次の表に示す通りとなった．
| 評価指標（損失関数） | 平均 | 標準偏差 |
|---|---|---|
| RI (constrative loss) | 0.972352 | 0.017127 |
| RI (triplet loss) | 0.986995 | 0.005775 |
| ARI (constrative loss) | 0.937757 | 0.038647 |
| ARI (triplet loss) | 0.970741 | 0.013020 |
baselineのスコアは大幅に向上しproposedとほぼ同等になったが，「平均」ではわずかながらproposedが勝る結果となった．以上より，baselineにおいても深層距離学習は問題なく動いていることが確認できた．
おまけ 埋め込みの可視化
まずbaselineに関して，図3と図4に埋め込みの可視化をそれぞれ示す．ただし訓練データに対する埋め込みの可視化であり，埋め込みネットワークはtriplet lossで学習したときのものを使っている． 図3は25000点すべてをプロットした場合，図4は埋め込み後に経験分布単位で平均ベクトルを計算した場合にそれぞれ対応している．
図3ではクラスタ間の重なりが大きいように見える．したがってクラスタリングのスコアも大幅に低下したと考えられる．図4の経験分布単位ではクラスタごとに分かれた埋め込みが得られており，やや広がりを持って散布している．このような埋め込みが得られた理由は，損失関数が分布単位で計算されておらず，あくまでも単独のサンプル点単位で計算されているから，つまり「分布単位で点がまとまる」ように変換は学習されないから，と一応の説明がつけられる．
最後にproposedに関して図5に埋め込みの可視化を示す．スコアの高さが示す通り，クラスタがきれいにまとまり，かつきれいに分かれた埋め込みが得られている．
実装
以下のリポジトリに置いた．深層距離学習の実装にはPyTorch Metric Learningが大いに役立った．動かし方の詳細はREADMEに書いてある．本記事では実装の詳細をいくつか抜粋して紹介する．
clustering/gaussianに置いてある．
moment-matching型 functional Bregman divergence のクラス
MomentMatching:
model.pyにある．PyTorch Metric LearningのBaseDistanceを継承し，compute_matとpairwise_distanceを独自に実装する．平均ベクトルを計算する処理が特別に入る以外はユークリッド距離の計算とほぼ同じであり，実装自体は容易である．データセットのクラス
DistributionDataset:
dataset.pyにある． moment-matching型のdivergenceの計算には経験分布に関する点集合を伴う． それら点集合たちをリストで管理する．リストの各要素が2次元の行ベクトル × 50個つまり[50, 2]のarrayである．実験条件の箇所で述べた通り，リストの全長は訓練用が500, テスト用が200である．なお従来法（深層距離学習のbaseline）はこのような特別のデータセットは使わずとも，全ての点をフラットに扱う一つの巨大なarray（例えば訓練用は[25000, 2]のサイズ）を作り，TensorDatasetでラップするだけで済む． ただし推論時に経験分布単位のフェアな比較をする場合はbaselineと提案法ともにDistributionDatasetを使う．評価指標を計算する
util.py
元論文ではRIとARIを計算していたが，scikit-learnには幸い"Clustering metrics"のカテゴリにて種々の指標を計算するAPIが提供されているので，併せて計算できるようにしておいた．ただしPurity scoreは自前の実装である．"Davis & Dhillon"の分布クラスタリングを実行する
entropic_clustering.py:
以前の記事の実装をそのまま持ってきて，今回の実験で用いたデータセット上で訓練・評価する．この記事はまさに本記事の布石であった．Python実装が転がっているはずもないので，自前で実装する必要があったからだ．
おわりに
今回の記事では分布クラスタリングに関して，実験の一部分について再現に成功した（と思う）．続く記事にて，分布クラスタリングに関して残された実験（"Human Activity Recognition Task"）についても報告する．
参考
著者らによるスライド＆プレゼン動画
slideslive.com