nykergoto’s blog

機械学習とpythonをメインに

文章の埋め込みモデル: Sparse Composite Document Vectors を読んで実装してみた

自然言語処理である単語の意味情報を数値化したいという場合に単語を特定のベクトルに埋め込む(分散表現)手法として word 2 vec があります。 この word2vec と同じような発想で文章自体をベクトル化するという発想があり Doc2Vec やそのたもろもろも方法が存在しています。 今回はその中の一つである SCDV (Sparse Composite Document Vector) を実装したのでその記録です。

そもそも何者か

文章を表現するベクトルを取得する手法です。

どうやってやるか

SCDV はいくつかのフェーズに分かれています。以下では5つのフェーズに分けて説明します。 若干論文の notation と違う所があるのでそこだけ注意していただければと思います。

1. 単語の分散表現を取得する

はじめに文章全体をつかって単語の分散表現を学習します。 Word2Vec や fasttext などが有名なところですね。

以下での話をわかりやすくするため、文章と単語に関してすこし数式を定義します. まず、文章中に現れるすべての単語の数を $N$ とおきます。 そして $i \in \{1, 2, ..., N\}$ 番目の単語に対して分散表現 $w_i \in {\mathbb R}^M$ を得たとします。ここで $M$ は埋め込みベクトルの次元数を表します.

2. 単語を更に複数のクラスタに分類する

1 で取得した単語ベクトルを更に $K$ 個のクラスタに分類します。 このときクラスタリングのモデルとして混合ガウス分布をつかいます。 これによって単語ごとに、先の埋め込みベクトルとは別のクラス分だけのベクトル $p_i \in \mathbb{R}^K$ をえることが出来ます。

混合ガウス分布の学習方法については論文中では特に言及はありませんでしたが、著者の実装では「各クラスタの共分散行列が同じである」という仮定のもとで推定を行うようになっていました。これはおそらく、各クラスタの領域を同じぐらいになるような制限をかけることで各単語の負担率 $p_i$ に極端な偏りがないようにするためなのかなーと考えています。単に計算量を減らすためかもなのであくまで推定ですが。

3. 埋め込み表現とクラス表現を掛けあわせた後に idf をかける

1 で得られた単語の分散表現に 2 で得られたクラスタへの割当を表すベクトルをおのおのかけて各単語ごとに $MK$ 次元のベクトルを作ります。 そして出来上がったベクトルに単語 $i$ の idf 特徴 ${\rm idf}_i \in \mathbb{R}$ を掛けあわせます ここで得られる新しい単語-クラスタ ベクトルを $u_i \in \mathbb{R}^{MK}$ とすると以下のようになります

$$ u_i = {\rm idf}_i \left( \begin{array}{c} p_{i 1} w_i \\ p_{i 2} w_i \\ \vdots \\ p_{i k} w_i \end{array} \right) $$

ここで idf を掛け算しているのは出現回数の多い単語の影響を低くしたいから、です。 とくに 2 でクラスタに分割して次元を多くしたので、出現回数の多い単語が大きい影響度を持つ次元の割合は増えることが予想されるためこの処理は必要なんだろうな、と僕は理解しています。*1

4. 文章ごとに単語クラスタベクトルを足しあわせて正規化する

ここではじめて文章という単位が登場します(ここまでの演算には文章が関係しません)。 $j$ 番目の文章に現れる単語を $L_j$ とし、全部で $J$ 個の文章があるとします。 先ほど得られた単語クラスタベクトルを文章ごとに足しあわせて文章ベクトル $v_j \in \mathbb{R}^{MK}$ を得ます

$$ v_j = \sum_{i \in L_j} u_i $$

さらにこれを各文章ごとに正規化します。論文には After normalizing the vector としか書かれていませんでしたが実装を見るとユークリッドノルムの意味で1になるようにしていました。 それに習うと, 正規化後のベクトルを $\hat{v_j} \in \mathbb{R}^{MK}$ とすると

$$ \hat{v_j} = \frac{v_j}{\| v_j \|_2} $$

となります。

5. ゼロに近いものを 0 に押しつぶす

4 の正規化を終えた時点でじつはほとんどの文章ベクトルの要素が 0 になります。 容量の圧縮や意味の鮮明化の為、ゼロに近いもののうち特定の条件を満たすものをゼロとみなします。具体的には

$$ \hat{v_i} = \begin{cases} \hat{v_i} & {\rm if} | \hat{v_i} | \ge t \\ 0 & {\rm otherwise} \end{cases} $$

とします。ここで $t$ は

$$ \begin{align} t = \frac{p}{100} \times \frac{ | a_{{\rm min}} | + | a_{{\rm max}} | } {2} \\ a_{{\rm min}} = \frac{1}{J} \sum_{j = 1}^{J} {\rm min}\ \hat{v_j} \\ a_{{\rm max}}= \frac{1}{J}\sum_{j = 1}^{J} {\rm max}\ \hat{v_j} \end{align} $$

のように定義される値です。全体の最大最小の $p$ %よりも小さかったら 0 にしましょうということみたいですね。(何故これで良いのかは特に記述がありませんでした。流石にノリで決めては無いと思うんですがよくわかりません)。

どのぐらいいいの?

定量評価のため論文中では SDCV と他の特徴量を SVM で訓練した結果が乗っています。

f:id:dette:20190224031410p:plain
Table1

これを見る限りわかりやすく一番つよいですね。

実装してみる

ここまで読んで楽しそうだったので自分でも実装をやってみました。

github.com

使うのはライブドアニュースのデータセットです。 docker-compose を使うと以下のように環境が作れます。

cp project.env .env
docker-compose build

docker-compose up -d

docker exec -it scdv-jupyter bash

SCDV の作成

作成は src/create.py で行います。著者実装のコードでは単語の分散表現もスクラッチで学習させていましたが今回は学習済みの fasttext 特徴量を用いています。 学習済みモデルは https://qiita.com/Hironsan/items/513b9f93752ecee9e670 よりお借りしました。感謝!!

def main():
    args = vars(get_arguments())

    word_vec = ja_word_vector()

    output_dir = os.path.join(setting.PROCESSED_ROOT)
    n_wv_embed = word_vec.vector_size
    n_components = args['components']

    docs, _ = livedoor_news()
    parsed_docs = create_parsed_document(docs)

    # w2v model と corpus の語彙集合を作成
    vocab_model = set(k for k in word_vec.vocab.keys())
    vocab_docs = set([w for doc in parsed_docs for w in doc])
    out_of_vocabs = len(vocab_docs) - len(vocab_docs & vocab_model)
    print('out of vocabs: {out_of_vocabs}'.format(**locals()))

    # 使う文章に入っているものだけ学習させるため共通集合を取得してその word vector を GMM の入力にする
    use_words = list(vocab_docs & vocab_model)

    # 使う単語分だけ word vector を取得. よって shape = (n_vocabs, n_wv_embed,)
    use_word_vectors = np.array([word_vec[w] for w in use_words])

    # 公式実装: https://github.com/dheeraj7596/SCDV/blob/master/20news/SCDV.py#L32 により tied で学習
    # 共分散行列全部推定する必要が有るほど低次元ではないという判断?
    # -> 多分各クラスの分散を共通化することで各クラスに所属するデータ数を揃えたいとうのがお気持ちっぽい
    clf = GaussianMixture(n_components=n_components, covariance_type='tied', verbose=2)
    clf.fit(use_word_vectors)

    # word probs は各単語のクラスタへの割当確率なので shape = (n_vocabs, n_components,)
    word_probs = clf.predict_proba(use_word_vectors)

    # 単語ごとにクラスタへの割当確率を wv に対して掛け算する
    # shape = (n_vocabs, n_components, n_wv_embed) になる
    word_cluster_vector = use_word_vectors[:, None, :] * word_probs[:, :, None]

    # はじめに文章全体の idf を作成した後, use_word だけの df と left join して
    # 使用している単語の idf を取得
    df_use = pd.DataFrame()
    df_use['word'] = use_words
    df_idf = create_idf_dataframe(parsed_docs)
    df_use = pd.merge(df_use, df_idf, on='word', how='left')
    idf = df_use['idf'].values

    # topic vector を計算するときに concatenation するとあるが
    # 単に 二次元のベクトルに変形して各 vocab に対して idf をかければ OK
    topic_vector = word_cluster_vector.reshape(-1, n_components * n_wv_embed) * idf[:, None]
    # nanで影響が出ないように 0 で埋める
    topic_vector[np.isnan(topic_vector)] = 0
    word_to_topic = dict(zip(use_words, topic_vector))

    np.save(os.path.join(output_dir, 'word_topic_vector.npy'), topic_vector)

    topic_vector = np.load(os.path.join(output_dir, 'word_topic_vector.npy'))
    n_embedding = topic_vector.shape[1]

    cdv_vector = create_document_vector(parsed_docs, word_to_topic, n_embedding)
    np.save(os.path.join(output_dir, 'raw_document_vector.npy'), cdv_vector)

    compressed = compress_document_vector(cdv_vector)
    np.save(os.path.join(output_dir, 'compressed_document_vector.npy'), compressed)

実験

著者の実験では SVM を使っていました。 普段モデリングするときには勾配ブースティング or Neural Network を使うことがおおいので、今回は lightGBM をつかってベンチマークを実装しました。 (実行するスクリプトsrc/SCDV_vs_SWEM.py にあります。)

タスクはライブドアニュースの文章から、そのカテゴリを予測する問題です。 カテゴリは全部で8個あり文章数はおおよそ7000程度です。

戦わせる特徴量は個人的に押しの論文 Baseline Needs More Love: On Simple Word-Embedding-Based Models and Associated Pooling Mechanisms で提案されている Simple Word Embedding Model (SWEM) をつかいます。

SWEM の文章 $k$ の特徴量 $z^k \in \mathbb{R}^{M}$ は以下で表されます

$$ z_j^{k} = \max_{i \in L_k} \ w_{ij}. $$

ここで $w_{ij}$ は単語 $i$ の $j$ 番目の要素であり $L_k$ は $k$ 番目の文章に含まれる単語の添字集合です. 要するに文章中の単語に対して埋め込み次元方向に max をとったものです。 得られる文章ベクトルは単語のものと同じ $M$ 次元になります。

これの n-gram version である SWEM-hier なども提案されていて文章分類などの単純なタスクにおいては CNN や LSTM をつかったリッチなモデルと匹敵、場合によっては勝つ場合もある、単純ですがあまり侮れない特徴量であったりします。

モデルはLightGBM で k=5 の Stratified Fold を行い accuracy により評価します。

結果と感想

Out of Fold の Accuracy を特徴量ごとにプロットしたのが以下のグラフです。

f:id:dette:20190224035914p:plain
学習結果

SCDV の圧縮を行わないバージョンがもっともよいスコアになりました。CV全てで二番目となった SWEM の max-pooling version よりも良い値となりこの特徴量の強さが伺えます。 一方 PCA で圧縮したものはあまりワークしませんでした。これはせっかく混合ガウス分布まで持ちだして拡張したことによって得られる微妙な差分が PCA によって押しつぶされてしまったことが原因と考えられます。

また SCDV は学習がおおげさになるというのも浮き彫りになりました。というのも特徴量の次元が $MK$ なので fasttext の埋め込み次元が 300, 混合ガウス分布の次元が 60 なので各文章ごとに 1800 次元もあるのです。

ですので今回のタスクのような小さい文章セット (1万以下) であっても numpy でなにも考えずに保存すると約 7GB 程度になります。また学習も SWEM に比べて体感 10 ~ 20 倍の時間がかかりました。 以上のことからメモリや計算資源に余裕が有る場合に SCDV を使いさくっとそれなりのベンチマークがほしい時は SWEM という使い分けも良いかも知れません。

まとめ

  • SCDV では単語情報を $M$ 次元ベクトルから更に $K$ 次元の情報を引き出して $MK$ 次元に拡張する。そこから先は普通(たして正規化)
  • livedoor ニュースのデータ・セットを用いた実験では SWEM よりも精度がよかった。が時間とメモリは要る。

*1:SWEMでは特に出現回数に関する考慮はなく精度が出ているため次元拡張の影響を打ち消すという意味合いが強いのかも

Pandas で Index を dictionary で更新したい

pandas のデータフレームで以下のようなものが有るとします。

In [1]: import pandas as pd
In [2]: import numpy as np

In [3]: df_train = pd.DataFrame(data=np.random.uniform(size=(3, 2)), index=['one
   ...: ', 'two', 'three'])

In [4]: df_train
Out[4]: 
              0         1
one    0.289881  0.649603
two    0.229427  0.811377
three  0.498204  0.779105

でこの index を例えば日本語の "いち", "に", "さん" になおしたいなという時。愚かな方法だと index.str.replace につらつらと書いていく方法があります.

df_train.index = df_train.index.str.replace('one', 'いち').str.replace('two', 'に')

でもなんかださいですね。dict で対応関係を記述して, それを元に変換するというふうにコードの責任分離をやりたいところです。

調べた所 pandas.DataFrame.rename を index にたいして使えとありました。

In [5]: en_ja_map = dict(one='いち', two='に', three='さん')

In [6]: df_train = df_train.rename(en_ja_map, axis='index')

In [7]: df_train
Out[7]: 
           0         1
いち  0.289881  0.6496030.229427  0.811377
さん  0.498204  0.779105

おしゃれでいいですね;)

さんこう

敵対的サンプリング検出のための基準としての相互情報量 - Understanding Measures of Uncertainty for Adversarial Example Detection

Understanding Measures of Uncertainty for Adversarial Example Detection
https://arxiv.org/pdf/1803.08533.pdf

概要 (200文字程度)

敵対的サンプルを判別する基準として相互情報量 (Mutual Information) が優れていることを主張する論文. MI の推定に Dropout を用いた MC を採用している.

Mnist と隠れ層2次元 Variational Autoencoder を用いて不確実性の分布を可視化し既存の基準である softmax entropy と MI の違いをわかりやすく説明している.

この論文の新規性

予測の不確実性を意味によって2つに分離し、モデルの知識不足部分の不確実性を取り出せれば敵対的サンプルを判別できるという主張を情報量の観点から議論している点.

不確実性の基準

そもそも不確実な入力には2つの種類がある

  1. 知識不足による認知の不確実性
  2. データの偶然性による不確実性

機械学習の文脈で言うと、前者はデータが欠乏しているために、その入力データに対する事後分布が平らになってしまうような現象を指す。

後者はそもそも入力データに対して出力のばらつきが激しくて予測が無理な場合を指す。 例えば表も裏も確率 1/2 で出るコイン投げの予測モデルなどは、知識ではなくデータの生成過程自体が完全にランダム性に支配されている為に起る不確実である。

一般に用いられる不確実性を測る基準として、予測値と入力のエントロピーがある。

$$ H[p(y | x)] = - \sum_{y \in Y} p(y| x) \ln p(y | x) $$

ここで $Y$ はモデルの予測値が取る離散的な空間を表している。 この基準は先の2つの不確実性を分離していない為、不確実性がモデルの知識不足から来るのか、それともデータの分散が激しいことから来るのかに関して判断することはできない。

そこでデータ $D$ と入力 $x$ が与えられた時のモデルのパラメータ $w$ とその予測値 $y$ の相互情報量 (Mutual Information) を考える。これは以下の式で表される。

$$ I(w; y | x, D) = H[p(y|x, D)] - {\mathbb E}_{p(w|D)}[H[p(y|w,x)]] $$

A, B の相互情報量とは B の情報を得た時に減少するAの不確かさのこと、と解釈することができる(A, B は逆でも成り立つ)。

今回でいうと、ラベル $y$ が与えられた時に減る $w$ に関する不確かさである。これは即ち $w$ に関して新しく得られた情報量とみなすことができる。

相互情報量が大きいような $x$ はどういう値かということを考えてみる。 相互情報量が大きいとは上記の議論より、得られる情報量が大きいということであるから、そのラベルの値 $y$ の価値が大きい(ひとつの $y$ を得るだけでとても大きな情報(知識)を得ることができる)。一方で相互情報量が小さいというのはラベルの価値が小さい、すなわち既にその入力 $x$ に対する $w$ の不確かさが十分に小さく、ラベルを与えられてもあまり情報が増えないということを表している。

よってこの相互情報量は、先の「不確実性の種類」で言うところの知識不足による不確実性(1) に相当する量である、とみなすことができる。

この情報量は直接的に計算できないが Dropout を用いた MC をすることで近似的に求めることが可能である。具体的にはラベルの事後分布を

$$ p(y|D,x) \simeq \frac{1}{T} \sum_{i=1}^T p(y|w_i, x) := p_{MC} (y|x) $$

の用に近似することで以下のように計算することができる

$$ I(w; y | x, D) \simeq H[p_{MC}(y | D, x)] - \frac{1}{T}\sum_{i=1}^T H[ p(y|w_i, x)] $$

不確実性の基準による違い

相互情報量であれ予測値のエントロピーであれ、入力画像 $x$ の分布(多様体)から著しく離れた $\hat{x}$ が与えられると大きな値を取る。 一方で予測値のエントロピーは入力画像 $x$ の分布(多様体) の近くであっても大きな値を取る場合がある。これはそういう画像が本質的に曖昧なもの, すなわちラベルがつけにくいものであるからである。

例えば人間が見ても 1 とも 7 とも見えるような画像があるとする。これは本質的に 1 or 7 が決定的に決められない曖昧な画像であるから予測値の softmax の値は複数のラベル(今回であれば 1 と 7) の確率が高い予測となる。エントロピーはひとつのラベルにピークを持つものほど小さくなる性質を持つので1、このような一つに決められないものに対してはエントロピーは大きくなる。 一方でMIは低い値を保つ(MIは (2) を測るものであるからである)。

さらに、敵対的サンプルとは元画像の持っているクラスとは別のクラスを持つ画像を作る行為であり、普通元のクラスの予測値を小さくし、高くしたいクラスの確率を大きくするように画像を修正する。

これは予測値のエントロピーを小さくするという効果もあるが MI によって測ることができる不確実性にも影響を与える、というのもエントロピーの下界が MI に相当するので、エントロピー最小化を行ってもMIより小さい値になることは無いからである。

MI と softmax の分散

softmax の分散と MI のテイラー展開した第一項は一致するので softmax の分散は MI の近似とみなすことができる。(詳細は論文参照のこと)

実験の対象と結果

MNIST に対して latend dim = 2 の VAE を用いて潜在空間の不確実性を可視化する実験と、犬と猫の分類問題に対する敵対的サンプルを entropy と MI とで予測しその不確実性の違いの特徴を調べる実験を行っている

MI と entropy の曖昧さの捉え方の違い

潜在空間で内挿したあと decoder を通して画像にしたものと、入力画像自体を内挿した画像とを用意してそれぞれ entropy と MI の値がどうなるかを表したのが以下の Figure2. である。

f:id:dette:20181128025644p:plain

上段の上の画像が潜在空間の内挿の画像で下が実空間上での内挿の画像を表している。
中段がそれらの entroy の値であり、下段が MI の値である。これを見ると entropy では潜在空間 (Latent Sapce) の内挿画像 (図の赤線) にたいしてピークを持っているのに対して MI ではピークがないことがわかる。これは潜在空間で内挿された画像は「不確実であるが本質的に曖昧で予測が曖昧なのは当然」な画像であるためデータの持つ偶然性も図ってしまう entropy は大きくなるが、モデルの知識に関する不確実性のみを測る MI ではそういう画像は曖昧で当然であるという判断になり小さいままなためである。

一方入力画像自体を内挿した場合 (Image Space) は entroy, MI ともにピークを持つ。これはそういう画像は訓練データに存在しておらず「この画像に対する知識は無い(言い換えるとモデルパラメータの分布の裾が広い)」ため、ラベルが与えられた時に得られる情報が大きいからである。

MNIST with VAE

隠れ層が2次元の VAE を用いて潜在空間上の点とその点での不確実性を可視化して、基準による違いを見る。 やり方は

  1. 潜在空間上で適当に点を取る
  2. decoder に通して画像空間上に射影して対応する画像を作成
  3. 作成した画像の不確実性を計算
  4. 以上を空間すべての点で計算して画像化する

という流れになっている。 fig1 が MI, fig4 が entropy の値を可視化したもの。どちらも白いところが不確実性が高い場所で黒くなるに連れて低い値となるようにプロットしている。色の付いている部分は MNIST の各クラスの画像の潜在空間上での値を表す。

f:id:dette:20181128030028p:plain f:id:dette:20181128030025p:plain

共通する傾向として中心から遠い(画像の多様体から遠い)場所では白い(=不確実性が高い)。

一方でよーく見てみると fig4 の方は色のついたクラスの間のところが白くなっているのがわかる。これは2つのクラスの間の画像で fig2 の上段の上の画像に相当するもので「数字っぽいけど本質的に曖昧でなんとも言えない画像」に対応している。 エントロピーだと不確実性を区別しないので、クラスの間は不確実性が高くなるというわけである。

もうひとつ気になるのは Dropout で不確実性をどの程度禁じできているのか、という点である。 これに関しては Dropout による変分推論は局所的なモードしか捉えられない為に事後分布全体の近似としては上手く働いていない。2 このため decoder によって生成された画像が意味をなしていないような場合でも高い信頼度を出してしまうような領域(holes と作者は呼んでいる)が出てしまう (figure5 がその例。 figure4 でも2時の方向に黒い holes が確認できる。)。3

このような穴があるがゆえ、敵対的サンプルはこの脆弱性をついて攻撃することが可能である。

そこで事後分布をちゃんと表現できればこの弱点を緩和できるのでは?という観点に基づいて修正を行う。 具体的には Dropout モデルをいくつも作ってアンサンブルする。こうすることで事後分布の mode の一つしか表現できないという弱点を解消できる。4

それをやったのが以下の figure6 で figure4 と比べるとだいぶ穴が減っているのがわかる

f:id:dette:20181128030125p:plain

犬猫分類

ASSIR cat and dogs データセットを用いて敵対的サンプルがどうかの判別を行い基準ごとの性質を確認する

使用するモデル

Imagenet で訓練済みの ResNet の Convolution 層を固定し Full Connection レイヤのみを fine-tuning したモデルをつかう。 Dropout は FC 層のみに適用し、 MC は 20 サンプルを抽出する。

確定的な予測は Dropout を使わないだけで重みの値は全く同じネットワークを使用する。

敵対的サンプルのアルゴリズム

  • Basic Iterative Method
  • Fast Gradient Method
  • Momentum Iterative Method

の3つを使う。

結果と考察

各種の敵対的サンプリング作成法によって作られた画像の分類精度を AUC で測ったのが以下の表である。

f:id:dette:20181128030140p:plain

DETERMINISTIC と MC の違いは予測時にも dropout を使って重みの事後分布の期待値として予測をしているかどうか、である(MCがテスト時にもdropoutを行う)。 相互情報量は dropout による近似を行わないと計算不能なため、DETERMINISTIC では N.A と記述されている。

結果を見ると MI が良い性能を出しておりまたエントロピーでは 0.50 を下回っている。 0.5 を下回る AUC は ランダムに判定するよりも性能が悪い ことになる。これは敵対的サンプルに対して過度に信頼度をおいていることを示している。

一方犬猫判定の精度においては、MCを用いたモデルの精度よりも確定的な値を用いたモデルの精度のほうが良いという結果になった。 これは Dropout の確率が高すぎたために畳み込み層の特徴の一部しか使えずに転移学習の足を引っ張る形になってしまうのが原因と考えられるので、 dropout は精度に関しては良い影響は与えられない。

議論と結論

  • 敵対的サンプルに対して制約を与えずにただ不確実性だけを見れば自動的に良いロバスト性を得られることを示したのは進展である。
  • Dropout だけが敵対的サンプルに対抗する手段であるとは思っていない。ただ確定的なモデルよりも攻撃するのが難しくなっただけ。
  • 隠れ層の可視化でわかったことは敵対的サンプルに対抗するためには第一にモデルのロバスト性を高めて不確実性を扱えるようにすれば、結果として騙されにくいモデルが出来上がる、ということ。

読んでみた感想

VAE を用いて潜在空間を可視化し、既存の基準と提案基準の違いをわかりやすく表現していて読みやすかった。 数式が出てくる論文は気分が良い。

次に読むべき論文は

敵対的サンプル系の論文読んだことないので読みたいね、ということで実験で使っていた手法一覧です。


  1. 最もエントロピーが小さいのはディラックデルタ関数となる場合であり、最もエントロピーが大きくなるのはすべてのクラスの確率が同じとなる場合である。

  2. これは KL距離が非対称であるがゆえ、VI が事後分布のピークの一つを近似するようなものになるという性質が原因

  3. (本当は訓練データに似たものの無い画像に対しては不確実だよーと言ってほしい)

  4. 当然計算量はモデル数が $n$ ならば $n$ 倍になってしまうという欠点は持っている

seaborn の clustermap をちゃんと理解する

python の可視化ツールは matplotlib が有名ですがそのプラグイン的役割を果たすモジュールとして seaborn があります。seaborn でプロットすると勝手に回帰直線をプロットしてくれたりするのですが今回は seaborn.clustermap を取り上げます。

seaborn.clustermap とは

seaborn.clustermap は入力された二次元データフレームを行と列それぞれでクラスタリングし、近いクラスタ同士を近い場所に配置するよう並び替えを行った後 heatmap + デンドログラム(樹形図)として可視化することができる関数です。

API: https://seaborn.pydata.org/generated/seaborn.clustermap.html

適当にあやめデータセットでやると以下のようになります。それっぽい図が出てきて楽しい気分になります。

import seaborn as sns
df_iris = sns.load_dataset('iris')
species = df_iris.pop("species")

g = sns.clustermap(df_iris)
g.savefig('clustermap.png', dpi=150)

f:id:dette:20181119022900p:plain

一方でこのクラスタリングがどういうロジックに基づいて動いているのか、カスタマイズするときはどうやったらええのか、ベストプラクティス的なことはないのか、など気になる点があるのでそれについてまとめて行きます。

ロジックについて

clustermap ではクラスタリングをすべて scipy.cluster.hierarchy.linkage にまかしているという記述があります。実際 clustermap 関数を呼び出すと

  1. Cluster Class インスタンスが新しく生成され、インスタンスの plot メソッドが呼ばれる
  2. plot 内部で self.plot_dendgram が呼ばれる
  3. self.dendgram では dendgram method が呼ばれ内部で _DendrogramPlotter インスタンスが新しく生成され、インスタンスの plot が呼び出される
  4. Dendgram.plot で scipy.cluster.hierarchy.linkage によって linkage が生成される

という流れで処理が走ります。というわけでまずは scipy の linkage を理解する必要が有ることがわかりました。

Scipy.cluster.hierarchy.linkage

linkage は階層型クラスタリングを行う関数です。

Scipy.cluster.hierarchy.linkage API scipy.cluster.hierarchy.linkage — SciPy v1.1.0 Reference Guide

具体的には特定の指標のベクトル $y$ を入力としてその hierarchy linkage を返します。 階層型クラスタリングとはすべての点が集まった状態を頂点として各点がひとつのクラスタになった状態を枝とするような木構造で表現できるクラスタリングのことです。階層型ではないクラスタリングで有名なのは k-means です。 k-means はある点が特定のクラスタに属するという状態が出力なためクラスタ同士に階層がありません。

入力値

一次元の距離尺度ベクトル、若しくは二次元の shape = (クラスタにしたい要素数, クラスタの特徴ベクトル) という形状の配列を指定します。

一次元ベクトルを渡す場合は、クラスタにしたい要素同士 $n$ のおのおのの距離が入った配列を指定する必要があります。そのため長さは必ず $nC_2$ である必要があります。

返り値

返り値は shape = (n-1, 4) の形状の matrix $Z$ です。
$Z$ の第 $i$ 番目の行ベクトルは $i$ 番目のクラスタリングの結果が入っています。

仮に $i$ 番目の行ベクトルを $v$ とします。この時 $v[0]$ と $v[1]$ は新しい $n+i$ 番目のクラスタとして結合される要素の index を表しています。 また $v[2]$ には結合されたクラスタ $v[0], v[1]$ の距離が、 $v[3]$ には結合後のクラスタの要素数がそれぞれ格納されています。

クラスタリングの指定

クラスタリングをどのように実行するか、を制御することができます。これはクラスタリングをする際に何を決めなくてはならないか、を考えるとわかりやすいです。

第一にクラスタリングする際にはクラスタごとの距離を計算する必要があります。これはアルゴリズム

  1. 今存在しているクラスタ同士の距離をすべて計算し
  2. そのなかでもっとも近い物同士を結びつける (linkage)
  3. 上記を繰り返す

というふうに進んでいく為 step1 でクラスタ間の距離計算を行う必要があるためです。

このクラスタごとの距離の計算を行う関数は「クラスタ化を考えているふたつの集合 $Y, Y'$ の中に存在する2つの要素 $u, v \in Y, Y'$ を入力としてその近さを返す距離尺度 $d: Y \times Y' \to \mathbb{R}$ 」である必要があります。

scipy にはいくつかその関数 (引数としては method で指定します) が用意されていますが、初めは単純な single を使います.
single は以下の数式で表される指標で, クラスタ内の要素の組み合わせすべての中で距離の最小値を返すものです

$$ d(u, v) = \min_{a \in u, b \in v}\left( {\rm dist}(a, b) \right) $$

ここで ${\rm dist}$ は2つのクラスタの要素を与えられた時に実数値を返す関数です. もし入力 y として1次元配列が渡された時は単に絶対値を返す関数となります。 一方で二次元配列 (n, m) が渡された時は m 次元ベクトル二点を入力としてその距離を返す関数 ${\mathbb R}^m \times {\mathbb R}^m \to {\mathbb R}$ を指定する必要があります。(scipy では metrics で指定します)

from scipy.cluster.hierarchy import linkage

y = [2, 8, 0, 4, 1, 9, 1, 0]

# y を特徴ベクトルとしたいので shape = (8, 1) に変換する
x = [[d] for d in y]
z = linkage(x, method='single', metric='euclid')

この場合だと dist はユークリッド距離になります。各要素の特徴次元が一次元の整数値なのでまあ大体の metric を使っても結果は変わりません。例えばチェビシェフ距離とかでも出力結果は一緒です。

>>> z_abs = linkage(x, method='single', metric='chebyshev')
>>> z_abs == z

    array([[ True,  True,  True,  True],
           [ True,  True,  True,  True],
           [ True,  True,  True,  True],
           [ True,  True,  True,  True],
           [ True,  True,  True,  True],
           [ True,  True,  True,  True],
           [ True,  True,  True,  True]])

つぎにこの出力値を解釈していきます。

>>> z

    array([[ 4.,  6.,  0.,  2.],
           [ 2.,  7.,  0.,  2.],
           [ 0.,  8.,  1.,  3.],
           [ 9., 10.,  1.,  5.],
           [ 1.,  5.,  1.,  2.],
           [ 3., 11.,  2.,  6.],
           [12., 13.,  4.,  8.]])

まず初めの行について考えてみます。

>>> z[0]

array([4., 6., 0., 2.])

index 0, 1 の値を見ると 4 番目と 6 番目の要素を結合していることがわかります(値はともに 1 です). それらの距離は 0 となっています. 要素がともに同じ値であるので明らかです。そしてこれらの要素は $8 + 0 = 8$ 番目のクラスタに割り当てられました。

>>> z[2]

array([0., 8., 1., 3.])

つぎに第2行目をみてみます。ここではクラスタ 0, 8 を結合しています。 クラスタ0 は即ち第0番めの要素でこれは 2 です。一方クラスタ8は8番目の要素(すなわち第1行目で出来たクラスタ)を指しているのでこれらの値はどっちも1でした。

2 と (1, 1,) のユークリッド距離の最小値は 1 なので距離は1です。そしてまとめられたクラスタは要素数3になり, 8 + 2 = 10 番目のクラスタに割り当てられました。

こんな感じにもっとも距離の近いクラスタを集結させていくのが linkage で行われている hierarchical/agglomerative clusteringというものです。

method の種類

method は"single" 以外にもいくつかの種類が用意されています.

https://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.linkage.html

single や complete といった最小値、最大値を取るものは特徴の中に一つでも近い(遠い)ものがあるときに極端なクラスタを作りがちなのでノイズや変数間のスケールに左右されてしまいます。一般には ward 法が外れ値やスケールに堅牢であるとされているようです。1

Seaborn.clustermap

scipy の linkage の理解が出来たところで seaborn.clustermap の方に戻りましょう. 主な引数は以下のようになっています.

基本

  • data: pandas.DataFrame
  • method: default は "average" . linkage のクラスタ間の距離尺度を変えたい時はこの値を変更します。scipy.cluster.hierarchy.linkagemethod と役割は同じです。
  • metric: default は "euclidien". 要素間の距離尺度を変えたいときはこの値を変更します。 scipy.cluster.hierarchy.linkage と役割は同じです。 例えば予測値同士の相関に関してクラスタにしてほしいなーという時は metric="correlation" とかにすると相関ごとにクラスタにしてくれるので良い様な気がします。
  • z_score: 数値が代入されるとその軸方向で 平均 0 分散 1 に正規化します。ずなわち 0 の時が行の方向で 1 の時が列方向です。行方向に正規化するのは画像データとかが該当するでしょうか(各画像ごとに平均ゼロ分散1に揃えたい)。普通の集計データの時は列方向を使うことが多いと思います(列方向に正規化すると各列の値がどれだけ平均からずれているかという無次元情報に変換されるため、特徴ごとのスケールなどを考慮しなくて済むからです)。
  • standard_scale: 数値が代入されるとその軸方向を $[0, 1]$ の範囲に圧縮(拡大)します。分散を1にしないスケーリングです。

ちょっとテクいことをする時

  • {row,col}_linkage: 行(列)の linkage を指定できます。デフォルトの挙動では行と列に同じ method 及び metric での linkage を行います。違うロジックを使いたいーという時には行の linkage, 列の linkage を予め scipy.cluster.hierarchy.linkage で計算しその結果を渡すとそのクラスタで linkage グラフを描写してくれます。

では実際にいろいろ試してみます。まずは average method でクラスタリングしてみます。

g = sns.clustermap(df_iris, method='average', cmap='viridis')
g.ax_col_dendrogram.set_title('Average Euclidean')
g.savefig('average_euclid.png', dpi=150)

f:id:dette:20181119030411p:plain

これでもいいんですが少し列間の類似度に差があまり無いように表示されています。これは列に対してユークリッド距離で距離を測っていてかつ正規化などをしていない為に列の値の平均値のズレに引っ張られている為です。 ちょっといたずらで petal_witdh の行に 100足してみましょう

df_fix = df_iris.copy()
df_fix['petal_width'] += 100
g = sns.clustermap(df_fix, method='average', cmap='viridis')
g.ax_col_dendrogram.set_title('Average Euclidean')
g.savefig('fixed_average_euclid.png', dpi=150)

f:id:dette:20181119031140p:plain

petal_width だけ仲間はずれになってしまいました。これはユークリッド距離を使っているために元データと比べ列比較の際に各要素の差分が 100 増えたことが原因です。

では次は距離をユークリッド距離ではなく相関係数にしましょう。相関係数の式はその系列の平均で引いた後に要素ごとに積を行うので平均を 100 増やしても問題ないはずです。

g = sns.clustermap(df_fix, method='average', metric='correlation', cmap='viridis')
g.ax_col_dendrogram.set_title('Average Correlation')
g.savefig('fixed_average_correlation.png', dpi=150)

f:id:dette:20181119031759p:plain

列のクラスタはいい感じになりました! でもちょっと色が意味をなしていないですね。というわけで今度は列方向に正規化します。

g = sns.clustermap(df_fix, method='average', metric='correlation', z_score=1, cmap='viridis')
g.ax_col_dendrogram.set_title('Average Correlation')
g.savefig('fixed_and_normalized_average_correlation.png', dpi=150)

f:id:dette:20181119031820p:plain

だいぶ見栄えも良くなりました。ちょっとあれな点があるとすると行方向にも相関を使っている所が微妙かも知れません。 というわけで行にはコサイン類似度を使ってクラスタリングをやるようにしてみます。

# 列方向は相関で
col_link = linkage(df_iris.values.T, method='average', metric='correlation')

# 行方向はユークリッド距離でクラスタ化
z = (df_fix.values - df_fix.values.mean(axis=0)) / df_fix.values.std()
row_link = linkage(z, method='ward', metric='euclidean')
g = sns.clustermap(df_fix, col_linkage=col_link, row_linkage=row_link, z_score=1, cmap='viridis')

g.ax_col_dendrogram.set_title('Average-Corr  Ward-Euclidiean')
g.savefig('average-corr__ward-euclidean.png', dpi=150)

f:id:dette:20181119032923p:plain

先よりは sepal_width の右上あたりに3あたりを取るクラスタっぽいものができるようになりました。

気をつけたほうが良さそうなこと

以上色々と試してきましたが、実際に可視化をする際に気をつけたほうが良さそうなポイントについて考えてみます。

scaling

今回の iris のような基本的なテーブルデータでは列ごとに意味を持っている場合が多いです。

そのためスケーリングを行うとするならば列方向の指定, 即ち z_score=1standard_scale=1 を指定するのが良さそうです。
また初めからスケーリングしてしまうと列ごとの値の大小がわからなくなってしまうため、列ごとの大きさの情報がわかっていない場合はスケーリングせずプロットしてみる、というのもありかも知れません。

アルゴリズムの選定

metric については列に対してはユークリッド距離ではなく相関係数を使ったほうが良いでしょう。ユークリッド距離を使うのであれば、列ごとのスケーリングもしくは平均値を均すような前処理を行うべきです。

method については色々と試してみましたが metric さえちゃんと選んでいれば average や ward など選んでおけば大体事足りそうです。2 一方で "single""double" のような最小値/最大値しか見ないロジックを指定すると出力結果が極端になりがちな印象を受けました。のでちょっといつ使ったらいいのかなーという気分です。

またテーブルデータでは列方向には相関 correlation を使いたい場面が多いと思いますが行方向 (レコードごとの距離測定) にも相関でいいのか?という純粋な疑問があります。3 なのでちょっと面倒ですが時間が有る場合は自分で linkage を計算してやるのがよいように思えます。


  1. http://ibisforest.org/index.php?Ward%E6%B3%95 がわかり易すぎるので自分の努力要らなかった説))

  2. ward を指定する際には metric は "euclidean" にする必要があります

  3. このあたりは SVM などデータ点どうしの距離を使うモデルでの前処理でも問題になるところだと思うのですが、未だに何がいいのかよくわかっていません。ノリで列方向に正規化したりはしていますが…