『Kernel Random Projection Depth for Outlier Detection』をPythonで実装した話

機械学習 python

はじめに

最近，下記の論文がアップロードされていた．

『Kernel Random Projection Depth for Outlier Detection』

arxiv.org

Pythonによる簡単なデモンストレーションを実装したので，公開する．

以降，Kernel Random Projection DepthをKRPDと呼ぶ．

KRPDの中身

ざっくり言えば，Kernel PCAによる非線形な変換をかけた後に，random projection depthを計算するアルゴリズムである． そのdepth値に-1を掛けた値を外れ値スコアとして利用する．

実験および実装

データセット

いくつかのtoy dataset上で，外れ値スコアの等高線を表示するデモンストレーションを用意した． 具体的には 以下の4つである．

  • Blobs: 3つのガウス分布から発生させたデータセット
  • Cross: 十字の形に分布するデータセット
  • Circles: scikit-learnに付属のmake_circles利用
  • Moons: scikit-learnに付属のmake_moons利用

それぞれのデータセットのサンプル点は総数が400個，うちinlierが300個，outlierが100個である．

実装

KRPDの実装およびデモンストレーションは以下のNotebookである．要PyODライブラリ．

その他の代表的な検知器（LOFやOCSVM，Isolation Forestなど; 下記参照）によってどのような等高線が得られるかを比較するNotebookたちも用意した． 検知器の実装はPyODライブラリによった．

ちなみに検知器たちは以下の8つである．

  • Angle-based Outlier Detector (ABOD)
  • Histogram-base Outlier Detection (HBOS)
  • $k$-Nearest Neighbor (KNN)
  • Local Outlier Factor (LOF)
  • Isolation Forest (IForest)
  • Cluster-based LOF (CBLOF)
  • Minimum Covariance Determinant (MCD)
  • One-Class Support Vector Machine (OCSVM)

実験結果

KRPDの等高線プロットを図1に示す．オレンジ色で示された領域は外れ値の比率（25%）から決められた，外れ値スコア上位の領域を示す．いずれのデータセットにおいても，この領域はinlierの分布をうまく囲うことができている．

図1：KRPDの等高線プロット

以下，図2から図5まで比較のための等高線プロットをそれぞれ示す．

図2：Blobs上での各検知器による等高線プロット

図3：Cross上での各検知器による等高線プロット

図4：Circles上での各検知器による等高線プロット

図5：Moons上での各検知器による等高線プロット

おわりに

Python （具体的にはPyOD）でKRPDを実装した．諸課題は論文に書かれている．

おまけ

InlierのみでKRPDをfitするColab notebookも作った．

Deep Divergence Learning (ICML 2020) の論文に掲載された実験結果を検証する試み −分布クラスタリング 前編−

深層学習 機械学習 距離学習

はじめに

距離学習および深層距離学習について述べた記事の最後に，

では距離関数自体を学習可能にする枠組みは存在するのか？→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.

proceedings.mlr.press

著者による実験のソースコードが長らく公開されておらず，空のリポジトリのみが存在していたが，2022年の2月にjupyter notebookが公開された．

github.com

残念ながら，この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となる（右下）．

図1: 提案手法（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を参照しながらデータの作り方を説明しよう．

図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を見て分かる通り，分散が大きいため各クラスタの点群は互いに大きく重なっている．

比較手法

の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へのページをリンクしておく．

  • Rand Index (RI) URL

  • Adjusted Rand Index (ARI) URL

これらはいずれも，同じデータセットに対する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の経験分布単位ではクラスタごとに分かれた埋め込みが得られており，やや広がりを持って散布している．このような埋め込みが得られた理由は，損失関数が分布単位で計算されておらず，あくまでも単独のサンプル点単位で計算されているから，つまり「分布単位で点がまとまる」ように変換は学習されないから，と一応の説明がつけられる．

図3：埋め込みの可視化（baseline，すべてのサンプル点）

図4：埋め込みの可視化（baseline，経験分布単位の平均ベクトル）

最後にproposedに関して図5に埋め込みの可視化を示す．スコアの高さが示す通り，クラスタがきれいにまとまり，かつきれいに分かれた埋め込みが得られている．

図5：埋め込みの可視化（proposed）

実装

以下のリポジトリに置いた．深層距離学習の実装にはPyTorch Metric Learningが大いに役立った．動かし方の詳細はREADMEに書いてある．本記事では実装の詳細をいくつか抜粋して紹介する．

github.com

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は自前の実装である．

    • Adjusted Mutual Information URL
    • Normalized Mutual Information URL
    • Homogenity score URL
    • Purity score
    • クラスタリング後の分類性能という観点で，$\mathit{k}$-Nearest Neighbor (kNN) 分類器に基づくaccuracyとROC AUCも計算できるようにしてある．

  • "Davis & Dhillon"の分布クラスタリングを実行する entropic_clustering.py:
    以前の記事の実装をそのまま持ってきて，今回の実験で用いたデータセット上で訓練・評価する．この記事はまさに本記事の布石であった．Python実装が転がっているはずもないので，自前で実装する必要があったからだ．

tam5917.hatenablog.com

おわりに

今回の記事では分布クラスタリングに関して，実験の一部分について再現に成功した（と思う）．続く記事にて，分布クラスタリングに関して残された実験（"Human Activity Recognition Task"）についても報告する．

参考

著者らによるスライド＆プレゼン動画

slideslive.com

*1:ちなみに指数型分布族と（正則な）凸関数の族が対応することが知られており，またBregman divergenceを与えたとき，それをKL-divergenceに持つ指数型分布族が存在することも知られている．

*2:Davis, J. and Dhillon, I. S. Differential entropic clustering of multivariate Gaussians. In Advances in Neural Information Processing Systems (NIPS), 2006.

"Differential Entropic Clustering of Multivariate Gaussians"をNumbaを使って高速化してみた話

機械学習 python

はじめに

前回記事で実装した Differential Entropic Clustering をもう少し高速化したいなぁ，という話．

tam5917.hatenablog.com

実装

やり方は簡単で，numbaをインストールして，@jit デコレータをBurg matrix divergence およびMahalanobis距離を計算する関数につけるだけ．

@jit(nopython=True)
def comp_burg_div(mat_x, mat_y):
    """Compute Burg matrix divergence."""
    dim = mat_x.shape[0]
    mat_y_inv = LA.inv(mat_y)
    mat = mat_x @ mat_y_inv
    burg_div = np.trace(mat) - np.log(LA.det(mat)) - dim
    return burg_div

@jit(nopython=True)
def comp_maha_dist(vec_x, vec_y, cov):
    """Compute Mahalanobis distance."""
    delta = vec_x - vec_y
    m = np.dot(np.dot(delta, LA.inv(cov)), delta)
    return np.sqrt(m)

前回記事のdemoスクリプトを上記の関数たちに置き換えて，time moduleで簡単に計測してみたところ，

  • 高速化なし ... 58.50秒
  • Burg matrix divergence を計算する関数のみ高速化 ... 33.59秒
  • Mahalanobis 距離を計算する関数のみ高速化 ... 45.33秒
  • 両方の関数を高速化 ... 12.31秒

となった．

高速化に成功した．行列計算関連で繰り返し呼び出される関数には，jitデコレータをつけておくのが良さそうである．

以下，参考まで，今回の実装．main関数の前後で時間計測する（面倒なので）．

A demo script of "Differential Entropic Clustering of Multivariate Gaussians" in Python. This is a faster version thanks to numba and jit. · GitHub

おわりに

Numbaによる高速化の効果は大いにあったので，薦めたい．

Differential Entropic Clustering of Multivariate Gaussians (NIPS 2006) をPythonで実装した

機械学習 python

論文はこれ．

proceedings.neurips.cc

多変量ガウス分布の平均と共分散行列の集合が与えられたときに，分布を単位として（＝対応する平均と共分散をペアにして）クラスタリングするアルゴリズムが提案されている． 行列の距離を測るための"Burg matrix divergence"が利用されているのが興味深い．これは行列を対象にしたBregman divergenceに相当する．

実装してみたのがこれ．要progressbar2．

Section 4.1 "Synthetic Data"を参考に，simplex上から平均ベクトルをクラスタ数だけランダムにサンプリングする．共分散行列は対角共分散にして，かつ対角成分の大きさも論文を参考にした．論文と設定が異なるのは，クラスタリングの対象とするのが標本平均と標本共分散という点である．つまり平均ベクトルと共分散行列で指定されるガウス分布から点を指定個数だけサンプルし，それらサンプル点に基づいて標本平均と標本共分散を計算する．論文よりも厳しい条件で分布クラスタリングを試していると言える．

クラスタリングの性能指標は以下を計算できるようにした．Purity Score以外は全てscikit-learnから利用できる． 参考までにAccuracyも計算している．

  • Rand Index
  • Adjusted Rand Index
  • Adjusted Mutual Information
  • Normalized Mutual Information
  • Purity Score
  • Homogenity Score

試しに動かしたらこんな感じ：

  • Rand Index = 0.776869 ± 0.071684
  • Adjusted Rand Index = 0.420437 ± 0.191200
  • Adjusted Mutual Information = 0.485865 ± 0.204090
  • Normalized Mutual Information = 0.494591 ± 0.200726
  • Purity Score = 0.630800 ± 0.135231
  • Homogenity Score = 0.485421 ± 0.196300

乱数の種を50回変えたときの結果（平均と標準偏差）である．データの次元数は4，クラスタ数も4である．今回はスコアのバラつきが大きくなった． スクリプト内のget_dataset関数をカスタマイズし，各クラスタの平均ベクトルの取り方や共分散行列の取り方を変えてみると，バラつきの傾向も変わるので色々と楽しめる．例えば，今回のように平均ベクトルをsimplex上に限定する必要はない．

クラスタリング結果の可視化はそのうち．

距離学習と深層距離学習の違い

距離学習 深層学習 機械学習

距離学習と深層距離学習の違いについて備忘録を残しておく．

共通していること

データセット中の2点 $\mathbf{x}, \mathbf{y}$間の距離関数$d(\mathbf{x}, \mathbf{y})$が与えられる． 距離関数はクラス分類やクラスタリング・異常検知などのタスクで活用される道具である． それらタスクに有利になる形で特徴量の変換を考えたい． その変換をパラメタ（行列など）で表現し，種々の損失関数 with 正則化項 に基づいてパラメタを最適化する． 距離に関する「望ましさ」は（不等式などの）制約によって表現される．例えば同一クラスのサンプルペアは近い距離，異なるクラスのサンプルペアは遠い距離にする制約である． したがって各サンプルには所属クラスなどのsupervisedなラベルが与えられることが典型である．

違い

代表的なユークリッド距離で考える．

  • 距離学習：$d(\mathbf{x}, \mathbf{y}) = ||G\mathbf{x} - G\mathbf{y}||_{2}$となる線形変換$G$（半正定値行列）を獲得する．

  • 深層距離学習：$d(\mathbf{x}, \mathbf{y}) = ||f_{\theta}(\mathbf{x}) - f_{\theta}(\mathbf{y})||_{2}$となる非線形変換$f_{\theta}$を獲得する． この非線形変換がニューラルネットにより実現され，$\theta$はそのパラメタである．

つまり距離学習は行列の最適化問題に帰着される．深層距離学習は当たり前だがニューラルネットのパラメタ最適化問題に帰着される． しばしば混同されるが，距離関数そのものはpre-definedであり最適化の対象ではない．この意味で距離関数の選択はハイパーパラメータといえる．

線形であれ非線形であれ，変換された特徴量はしばしば「埋め込み」（embedding）とも呼ばれる． 距離学習では特徴量エンジニアリングが依然として重要であるが，深層距離学習では埋め込みの獲得が特徴量抽出を兼ねている面がある．

では距離関数自体を学習可能にする枠組みは存在するのか？→Yes.

日本語x-vectorから感情成分を分離するニューラルネットワークを構築してみた −感情分類に敵対的な損失関数の導入−

音声 深層学習 機械学習

はじめに

本記事は前回記事の続編に相当する．

前回記事では声優統計コーパスの3話者・3感情の音声データに対してx-vector抽出器を適用し，UMAPで可視化を試みた． この可視化の実験を通じて，感情成分が分離できていない傾向が見られた．すなわち，本来は話者3クラスにも関わらず，疑似的な9クラス（= 3話者 × 3感情） が存在するように見える，というものである（x-vector抽出器の学習データを考えてみれば，それはそうなのだが）．せっかくx-vectorが手元にあるのだから，感情成分を分離/除去するフィルタの役割を果たす手法を実装してみたいと考えた．本記事はその実装の詳細と簡単な検証実験に関する報告である．

感情成分を分離するニューラルネットワーク

先行研究と論文

今回の実装にあたり下記の論文を参考にした．本論文では，音響特徴量（ベクトル系列）に含まれる話者成分とテキスト情報を表す成分を分離するよう，損失関数に制約が加えられたニューラルネットワークによってnon-parallelな音声変換を実現している．この「損失関数とネットワーク構造を工夫することで特徴量に含まれる特定の成分を分離する」というコンセプトが， 今回のx-vectorを用いた感情成分の分離にも使えるはずだ，という個人的な直感があった．

Jing-Xuan Zhang, Zhen-Hua Ling, and Li-Rong Dai. 2020. Non-Parallel Sequence-to-Sequence Voice Conversion With Disentangled Linguistic and Speaker Representations. IEEE/ACM Trans. Audio, Speech and Lang. Proc. 28 (2020), 540–552. https://doi.org/10.1109/TASLP.2019.2960721

ネットワーク構造，損失関数およびモデル訓練アルゴリズム

ネットワーク構造の概略を図1に示す．

図1: ネットワーク構造

ねらい

仮にx-vectorに話者分類器をくっつけてsoftmax層の直前のベクトル（ logitではない）を取り出せば，感情成分は話者ラベルに丸めこまれる形となって，見かけ上で感情成分は消えてしまう．しかし，それでは面白くない．話者分類器を手元のx-vector群に関してファインチューニングすると，x-vectorに埋め込まれていることが期待される多様な発話スタイルを含む話者特徴が消えてしまうおそれがある．逆に感情分類器のみをくっつけると，今度はもとのx-vectorを忘れてしまう．つまり感情分類器は入力がx-vectorいかんに関わらず構築できるため，話者特徴が消えてしまい，本来の目的から遠ざかる． そこでまず，話者成分と感情成分をそれぞれ独立に抽出するencoderを用意する（speaker encoderとemotion encoder）．各encoderで抽出された各成分をdecoderに通し，もとのx-vectorを復元する制約を付与する．これによりencoder出力とx-vectorとの間の暗黙的な相互情報量の最大化を期待している．emotion encoder側には感情分類器（線形層 + softmax）をくっつけることで，感情成分の抽出を促進するねらいがある． 一方，speaker encoder側にも，補助の感情分類器（auxiliary classifier）をくっつける．ただしこのclassifierが強くなりすぎるとspeaker encoderの出力に感情成分が混在してしまうので，auxiliary classifierの感情分類に敵対的な損失関数（adversarial loss）を導入し，感情成分の分離・除去を促進する．つまり，もし仮に話者成分と感情成分が混在したx-vectorから感情成分がすっかり除去されれば，auxiliary classifierの予測はすべて等確率（3感情の分類ならばすべて1/3）になるはずだという直感のもと，敵対的損失関数により感情の予測分布が一様分布へと近づくようにclassifierおよびencoderを成長させる．

損失関数たち

図1より，以下の損失関数を考える：

  • $\mathcal{L}_{rec}$ : encoder-decoderによる再構成損失（平均二乗誤差）
  • $\mathcal{L}_{enc}$ : emotion encoder側の感情分類に関する損失（交差エントロピー）
  • $\mathcal{L}_{cls}$ : speaker encoder側の感情分類に関する損失（auxiliary classifierの交差エントロピー）
  • $\mathcal{L}_{adv}$ : speaker encoder側の感情分類に関する敵対的損失

$\mathcal{L}_{adv}$の計算には選択の余地がある：

(1) 一つは auxiliary classifierに関する感情クラスの予測分布を$\hat{\mathbf{p}}$とし，一様分布を表すベクトルを$\mathbf{u}$とするとき，$\hat{\mathbf{p}}$と$\mathbf{u}$の平均二乗誤差を損失に設定する方法である．その場合の損失は次式で書くことができる：

$$ \mathcal{L}_{adv} = \frac{1}{N} \sum_{i = 1}^{N} ||\hat{\mathbf{p}}^{(i)} - \mathbf{u}||_{2}^{2} $$

ただし$N$はミニバッチのサイズであり，一様分布を$\mathbf{u} = [1/C, 1/C, \ldots, 1/C]$というベクトルで表現している．$C$は感情クラスの数を表し，今回は$C = 3$である．したがって$\mathbf{u} = [1/3, 1/3, 1/3]$である．$\hat{\mathbf{p}}$の肩に書いてある添字はミニバッチ内の$i$番目のサンプルに関して計算した予測であることを明示する．平均二乗誤差の最小化により，予測分布を一様分布に近づける，すなわちどの感情クラスも均等に分類を誤るようauxiliary classifierに制約を課す．

(2) もう一つの選択は$\mathcal{L}_{cls}$の符号を反転したものを$\mathcal{L}_{adv}$に設定する方法である．すなわち$\mathcal{L}_{adv} = -\mathcal{L}_{cls}$とする．これはauxiliary classifierが相反する損失関数を持つという意味で明解である．後述する実装ではこれら敵対的損失をオプションで切り替えられるようにした．

損失関数に基づく訓練アルゴリズム

ネットワーク全体のうち，auxiliary classifier以外の部分を"main"と称する．

  1. x-vectorをネットワークに流すことで$\mathcal{L}_{rec}$ , $\mathcal{L}_{enc}$, $\mathcal{L}_{cls}$, $\mathcal{L}_{adv}$を計算する．

  2. "main"の損失$\mathcal{L}_{\mathrm{main}} = \mathcal{L}_{rec} + \mathcal{L}_{enc} + w_{adv} \mathcal{L}_{adv}$を計算し，この損失関数に基づいて勾配を流しネットワーク重みを更新する．その間，auxiliary classifierのネットワークパラメタは固定（freeze）する．ここで$w_{adv}$は$\mathcal{L}_{adv}$の影響を調整する重み係数である．

  3. auxiliary classifierのパラメタの固定を解除し，逆に"main"のパラメタを固定する．$w_{cls} \mathcal{L}_{cls}$に基づいて勾配を流し，auxiliary classifier のパラメタを更新する．その後，"main"のパラメタの固定を解除する．ここで$w_{cls}$は$\mathcal{L}_{cls}$の影響を調整する重み係数である．

  4. 指定回数だけ1. から 3. を繰り返す．

実装

以下のリポジトリに置いた．Enjoy!

github.com

動かし方

  1. config.yaml を各自の環境に合わせて修正する．root_dirの修正のみで十分と思われる．

  2. $ python preprocess.py: 前処理を実行する． 具体的には，声優統計コーパスのダウンロード（および解凍），x-vector抽出のための事前学習済モデルのダウンロード，および声優統計コーパスの音声群からx-vectorを抽出して保存する処理からなる．

  3. $ python training.py: モデルの訓練を実行する．

  4. $ python inference.py: 訓練済みのモデルにx-vectorを通すことで，感情成分が除去されたx-vectorを手に入れ，保存する．

  5. $ python plot_umap_all.py: UMAPによりx-vectorを可視化する．

実験

実験条件

前回記事と同様に，声優統計コーパスの音声サンプルを利用している．敵対的損失関数は感情分類の予測分布と一様分布の間の平均二乗誤差とした． そのほか詳細は折りたたみ形式で示す。ネットワークの次元数など、どうしても気になるひとだけクリック。

実験条件を表示する

計算機環境を示す．

計算機環境 バージョンなど
OS Ubuntu 22.04
CPU Intel i9-9900K
GPU RTX2070
Python 3.10.6
PyTorch 2.0.1

訓練の基本的な設定は以下の通りである．

項目 設定
ミニバッチサイズ 32
エポック数 1500
学習率 0.001からスタートし，500エポックおきに0.6倍

ニューラルネットワーク全体を通じて，すべて全結合層からなる．speaker encoderにおける各layerの次元数は 入力側から順に(512)-(256)-(256)とした．emotion encoderも同様の次元数である．emotion encoderに直接くっつけるclassifierはlogit計算のための線形層が1つ挟まっている．auxiliary classifierはencoder側からsoftmax層に向かって(256)-(256)-(256)-(256)-(3) である．decoderの次元数はencoder側から出力方向に向かって(256+256)-(512)-(512)-(512)-(512)である．次元数を256+256と表現しているのはspeakerとemotionのencoder出力（それぞれ256次元）を結合しているためである． これら全結合層にはバッチ正規化が組み込まれており，また活性化関数にはReLUを採用した．

以上の設定はconfig.yamlに書いてある．

実験結果

図2に可視化結果を示す．話者tsuchiyaのごく一部のサンプルが話者fujitouのクラスタに紛れ込んでいるが，概ね話者ごとに3クラスタに分かれており，感情成分の分離はある程度成功したのではないかと考えられる．話者uemuraや話者fujitouは各感情のサンプルがよく重なっており美しいが，話者tsuchiyaはangryクラスが依然として独立する傾向にある．実際の音声サンプルをさらに詳細に調べることで，このような傾向の違いが生じた原因が見えてくる可能性はある．その際，発話長や抑揚， ボリュームなど，個々の発話から計算できる特徴を統計的な分析にかけることが有用であると考えられる（例：平均発話長，基本周波数の平均および分散）．

図2: 感情成分を分離したx-vectorのUMAPによる可視化

TIPS

敵対的損失関数にかかる重み係数の調整には少し気を使うが，ある程度の範囲で係数を変えても同様の結果は得られそうである（係数の値の変化にはあまり敏感でないかも？）．またモデル訓練のエポック数も多くしすぎるのは避けたほうがよい．実験結果はエポック数1500だったが，3000回では各話者のx-vectorが連続的につながるような可視化結果を得ている（図3）．ほどほどのエポック数で止めておくのが良さそうである．各ネットワークのユニット数（各レイヤーの次元数） や層数も調整の余地はあるだろうが，あまり次元数を絞りすぎるとencoderやdecoder，auxiliary classifierが弱くなり，その場合も期待通りの分離結果は得られない．

図3: エポック数3000のときの結果

おわりに

本記事では日本語x-vectorから感情成分を分離するニューラルネットの実装を試みた． 感情成分を分離するために，感情分類器に敵対的な損失関数を導入した． UMAPによる可視化実験の結果，前回記事におけるx-vectorの疑似9クラス（= 3話者 × 3感情）が概ね話者ごとに独立した3つのクラスタを作るようになり，本記事の試みはある程度成功したといえるだろう．

感情成分が分離されたx-vectorを何らかの後段のタスクに使ってみるのも一興だろう．今回は幸いにも話者ラベル・感情ラベルが利用可能なので，各種の距離学習の導入も考えられる．そのほか，教師なし（= ラベル情報なし）で分離する手法への発展も考えられるだろう．今後の課題としたい．

参考文献

冒頭で示した論文の公式実装．実装面で大いに参考になった． github.com

日本語x-vector抽出器により声優統計コーパスからx-vectorを抽出してUMAPで可視化した話

深層学習 音声 python

はじめに

最近，日本語 x-vector 抽出器がPyPIに登録された旨がツイートされた．

そこで早速，『声優統計コーパス』の音声データに対して適用し，抽出したx-vectorを可視化しようと試みたので，その報告が本記事の主旨である．

日本語 x-vector 抽出器

抽出器はxvector-jtubespeechという名前で，すでにgithubのリポジトリが公開されている．

github.com

x-vector (話者表現ベクトル) を抽出するための学習済みモデルを提供します．このモデルは，JTubeSpeechコーパスと呼ばれる，YouTubeから収集した日本語音声から学習されています．

ということなので，ぜひ使ってみたいと．

声優統計コーパス

下記のサイトにて公開されている．3名の話者が「通常」「怒り」「喜び」の3パターンの感情で音素バランス文を読み上げた際の音声データである．オリジナルの音声が48kHzサンプリングで収録されているので，x-vector抽出の際にはlibrosaなどにより16kHzにリサンプリングしておく必要がある．なぜならば，xvector-jtubespeech の事前学習済モデルは16kHzの音声で訓練されているため，入力音声を16kHzサンプリングに揃えなければならないからである．

voice-statistics.github.io

実装

以下のリポジトリにサンプルコードを置いた．

github.com

具体的には，

  • preprocess.py

    • 声優統計コーパスのダウンロード（および解凍）
    • 事前学習済モデルのダウンロード
    • x-vectorを抽出して保存

  • plot_umap.py

    • UMAPによりx-vectorを可視化（話者ごとに図を作成）

  • plot_umap.py

    • UMAPによりx-vectorを可視化（全話者をまとめて1つの図として作成）

  • config.yaml

    • 設定ファイル

からなる．事前学習済モデルは上記xvector-jtubespeechのリポジトリからダウンロードしている．config.yamlは各自の環境に応じて書き換えられたい（root_dirのみで十分）．

UMAPによる可視化スクリプトの例としてColab Notebookも用意した．

colab.research.google.com

可視化例

図1に声優統計コーパス3話者から抽出したx-vectorを可視化した例を示す．3話者とは「fujitou, tsuchiya, umemura」の3名である．話者ごと・感情ごとにx-vectorがきれいに分かれてプロットされていることが見て取れる．

図1. UMAPによるx-vectorの可視化（声優統計コーパス3話者）

おわりに

x-vector抽出器の公開はとても嬉しい．いろいろな音声データセットに対して抽出を試み，傾向の違いなどを観察するのは興味深いことだろう．