次元削減 - 静かなる名辞

【python】scikit-learnで大規模疎行列を扱うときのTips

2019-08-14T02:33:31+09:00

はじめに

　自然言語処理などで大規模疎行列を扱うことがあります。一昔前はNLPといえばこれでした（最近は低次元密行列で表現することのほうが多いですが）。

　疎行列はその特性をうまく生かして扱うとパフォーマンス上のメリットが得られる反面、うかつにdenseな表現に展開してしまうと効率が悪くなって激遅になったり、あっさりメモリから溢れたりします。

　scikit-learnでやる場合、うっかり使うと自動的にdenseな表現に展開されてしまう、という事故が起こりがちで、要するに使えるモデルに制約があり、注意が必要です。その辺の基本的なことをまとめておきます。

はじめに
疎行列ってなに？
特徴抽出する
特徴選択する
標準化する
次元削減する
その他の各種transformerを使う
分類や回帰など、予測に使う
実際にやってみる
まとめ

疎行列ってなに？

　まず、Pythonで疎行列といえばscipy.sparse.csr_matrixなどのことです。

scipy.sparse.csr_matrix — SciPy v1.4.1 Reference Guide

　内部の構造の詳細などは、こちらの記事が参考になります。

scipy.sparseの内部データ構造 – はむかず！

　重要なのは、まずこの方式にすると「メモリ効率がいい」ということです。これは単純に嬉しいですし、CPUキャッシュのことを考えても、パフォーマンス上大きなメリットがあります。また、0の要素は飛ばして探索できるので、うまく使うと効率も良くなります。

　大規模疎行列を相手にするときは、できるだけ疎行列のまま取り扱うことが重要になります。

特徴抽出する

　自然言語処理系のタスクだと、CountVectorizerやTfidfVectorizerが使えます。どちらもデフォルトでcsr_matrixを返してくれるので、素直に使えば疎行列が出てきます。

sklearn.feature_extraction.text.CountVectorizer — scikit-learn 0.22.1 documentation
sklearn.feature_extraction.text.TfidfVectorizer — scikit-learn 0.22.1 documentation

　もう少し幅広いタスクで使いたい場合は、DictVectrizerが便利でしょう。こちらもデフォルトではsparseな表現を返します（オプションでnumpy配列を返すようにすることも可能）。

sklearn.feature_extraction.DictVectorizer — scikit-learn 0.22.1 documentation

特徴選択する

　特徴抽出したあと素直に使うとだいたい変数が多すぎて使いづらいので、普通は変数選択をすると思います。sklearn.feature_selectionのものなら、これはだいたいデフォルトで疎行列のままの取り扱いに対応しています。

1.13. Feature selection — scikit-learn 0.22.1 documentation

　疎行列としてinputすれば疎行列で出てきます。速度もそうした方が速いです。

sklearnの変数選択は疎行列型（csr_matrix）でやると速いっぽいよ - 静かなる名辞

標準化する

　StandardScalerで標準化する場合は、with_mean=Falseを指定してください。これは平均0にしない標準化です。標準化の式の分子で平均を引かないものです。

　それだけで疎行列型のまま標準化することができます。

This scaler can also be applied to sparse CSR or CSC matrices by passing with_mean=False to avoid breaking the sparsity structure of the data.
sklearn.preprocessing.StandardScaler — scikit-learn 0.22.1 documentation

scikit-learnのStandardScalerで疎行列型のまま標準化する - 静かなる名辞

次元削減する

　Truncated SVDという素晴らしい手法があり、実装上も疎行列に対応しているので、こちらを使ってください。逆に、これ以外の選択肢は（おそらく）ありません。

sklearn.decomposition.TruncatedSVD — scikit-learn 0.22.1 documentation

　ただし、次元削減した方が良いのか、しない方が良いのかはなんとも言えません。次元削減は行わないで、疎行列型のまま後段のモデルに突っ込むという選択もあるからです。ぶっちゃけ性能は大して変わらないし、次元削減に時間がかかるのと大規模密行列になってしまうぶんだけ遅くなるかもしれない……という微妙な性質があります。

　それでも、次元削減が必要ならやればできます。

その他の各種transformerを使う

　transformの返り値がsparse matrixになるかどうかを確認してください。油断しているとnumpy配列に変換されます。

　リファレンスからある程度は読み取れますが、わからないことも多いので、一度動かして確かめた方が良いと思います。

分類や回帰など、予測に使う

　Truncated SVDで次元削減をした場合は勝手にnumpy配列になっているので、どんなモデルにも入力できます（実用的な速度と性能が両立できるかは別）。

　csr_matrixのまま突っ込む場合は、そのまま入力できるestimatorとできないestimatorがあるので、注意が必要です。これを確認するには、リファレンスのfitメソッドのパラメータを見ます。

　たとえばRandomForestClassifierだと、

X : array-like or sparse matrix of shape = [n_samples, n_features]
3.2.4.3.1. sklearn.ensemble.RandomForestClassifier — scikit-learn 0.22.1 documentation

　という記述があり、sparse matrixと書いてあるのが「疎行列型でも受け付けて、適切に取り扱ってあげますよ」という意味です。一方、たとえばLinearDiscriminantAnalysisだと（あまり使う人もいないと思いますが）、

X : array-like, shape (n_samples, n_features)
sklearn.discriminant_analysis.LinearDiscriminantAnalysis — scikit-learn 0.22.1 documentation

　と書いてあります。array-likeのときは、渡せば動くけど、内部的にはdenseな表現（numpy配列）に変換されてしまう公算が大きいです。でもけっきょくのところはよくわからないので、実際に入れて動くかどうか試してみた方が良いでしょう。

　他に実例は省略しますがnumpy arrayと書いてあるときも（たぶん）あり、この場合はたぶんsparse matrixだとエラーを吐きます。

　実際に動かしてみないと挙動がわからないこともままあるので、突っ込んでみてエラーが出ないか、メモリ消費が異常に膨れ上がらないかを確認しながら作業した方が良いと思います。

　以下は疎行列型でも行ける（と思う）代表的なestimatorのリストです。

　分類器

sklearn.ensemble.RandomForestClassifier
sklearn.svm.SVC
sklearn.svm.LinearSCV
sklearn.naive_bayes.MultinomialNB

　非負のみ。文書分類向き

sklearn.linear_model.LogisticRegression

　回帰モデル

sklearn.ensemble.RandomForestRegressor
sklearn.svm.SVR
sklearn.linear_model.ElasticNet

　代表的なものはだいたい対応しているけど、たまに使えないのもあるという感じです。

実際にやってみる

　20newsgroupsの文書ベクトルを返してくれるものがあるので、それでやります。

Classes 20
Samples total 18846
Dimensionality 130107
Features real
sklearn.datasets.fetch_20newsgroups_vectorized — scikit-learn 0.22.1 documentation

sklearnのfetch_20newsgroups_vectorizedでベクトル化された20 newsgroupsを試す - 静かなる名辞

　ご覧の通りでかいので、ナイーブにnumpy配列に変換して扱おうとすると苦労します。前にやったときは密行列に変換しようとしていろいろ苦労していましたが、疎行列型のままやった方がシンプルです。

　もちろん通例通り、Pipelineを使ってモデルを組み合わせます。

【python】sklearnのPipelineを使うとできること - 静かなる名辞

from sklearn.datasets import fetch_20newsgroups_vectorized
from sklearn.feature_selection import SelectKBest
from sklearn.decomposition import TruncatedSVD
from sklearn.svm import LinearSVC
from sklearn.pipeline import Pipeline
from sklearn.metrics import classification_report

def main():
    train = fetch_20newsgroups_vectorized(subset="train")
    test = fetch_20newsgroups_vectorized(subset="test")
    X_train, y_train = train.data, train.target
    X_test, y_test = test.data, test.target
    target_names = train.target_names

    skb = SelectKBest(k=5000)
    tsvd = TruncatedSVD(n_components=1000)
    svm = LinearSVC()
    clf = Pipeline([("skb", skb), ("tsvd", tsvd), ("svm", svm)])
    clf.fit(X_train, y_train)
    y_pred = clf.predict(X_test)
    print(classification_report(
        y_test, y_pred, target_names=target_names))

if __name__ == "__main__":
    main()

　実行したところ、1分くらいで処理が完了しました。パフォーマンスのよさが伺えます。

　結果。

                          precision    recall  f1-score   support

             alt.atheism       0.70      0.68      0.69       319
           comp.graphics       0.70      0.71      0.70       389
 comp.os.ms-windows.misc       0.71      0.71      0.71       394
comp.sys.ibm.pc.hardware       0.66      0.64      0.65       392
   comp.sys.mac.hardware       0.76      0.76      0.76       385
          comp.windows.x       0.79      0.72      0.75       395
            misc.forsale       0.81      0.87      0.84       390
               rec.autos       0.83      0.83      0.83       396
         rec.motorcycles       0.91      0.90      0.90       398
      rec.sport.baseball       0.84      0.89      0.87       397
        rec.sport.hockey       0.92      0.95      0.94       399
               sci.crypt       0.89      0.88      0.89       396
         sci.electronics       0.66      0.65      0.66       393
                 sci.med       0.81      0.78      0.79       396
               sci.space       0.86      0.87      0.86       394
  soc.religion.christian       0.76      0.91      0.83       398
      talk.politics.guns       0.68      0.88      0.77       364
   talk.politics.mideast       0.91      0.81      0.85       376
      talk.politics.misc       0.71      0.54      0.62       310
      talk.religion.misc       0.61      0.41      0.49       251

                accuracy                           0.78      7532
               macro avg       0.78      0.77      0.77      7532
            weighted avg       0.78      0.78      0.78      7532

　今回は性能を重視していないのでこの程度です。このタスクだとできるだけ次元を維持したまま（疎行列型のまま）ナイーブベイズに入れたほうが速くて性能が出るという知見を以前に得ています。その場合は0.83くらいまで出ています。

【python】sklearnのfetch_20newsgroupsで文書分類を試す(5) - 静かなる名辞

まとめ

　うまく疎行列型配列を使うと数桁くらい時間を節約できます。ぜひ活用してみてください。

　こちらの記事もおすすめです。
scikit-learnのモデルに疎行列（csr_matrix）を渡したときの速度 - 静かなる名辞
　

カーネルPCAで文字列の編集距離を可視化してみる

2019-07-08T01:20:53+09:00

はじめに

　以前に編集距離が計算された文字列間の位置関係をMDSを使ってまったく同じことをしましたが、今度はカーネルPCAでやってみます。

　違いとしては、MDSは距離行列から計算を行うのに対してカーネルPCAは類似度行列から計算を行えるということがあると思います。逆に言えば、それくらいしか違わないし、だいたい同じようなことができるということです。

　なお、この記事を読む前に、以下の2記事を先にご覧いただくとやっていることがわかりやすいと思います。

多次元尺度構成法（MDS）で文字列の編集距離を可視化してみる - 静かなる名辞
 scikit-learnのSVMを自分で計算したカーネルで使う - 静かなる名辞

実験

　カーネルPCAのドキュメントはこちらです。

sklearn.decomposition.KernelPCA — scikit-learn 0.21.3 documentation

　kernel="precomputed"として、あとは淡々と実装します。

import Levenshtein
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import KernelPCA

def main():
    # データ（自作）
    X = ["hogehoge", "hogghheg", "fogefoge", "hagbhcgd", "hogeratt",
         "hohohoho", "fugefuge", "hokehoke", "hogehope", "kogekoge",
         "fugafuga", "fugafuge", "fufufufu", "faggufaa", "fuuuuuuu",
         "fhunfhun", "ufagaguf", "agufaguf", "fogafoga", "fafafaoa"]
    label = np.array(["hoge"]*10 + ["fuga"]*10, dtype=object)

    # 類似度行列の作成（levenshtein距離）
    A = np.zeros((len(X), len(X)))
    for i in range(len(X)):
        for j in range(i + 1, len(X)):
            d = Levenshtein.distance(X[i], X[j])
            A[i,j] = A[j,i] = d
    A = -A  # 距離->類似度変換のためマイナスをかける

    # kernel pcaによる変換
    kpca = KernelPCA(n_components=2, kernel="precomputed")
    X_2d = kpca.fit_transform(A)

    # plot
    for (x, y), l in zip(X_2d, X):
        plt.text(x, y, l)
    
    for l in set(label):
        plt.scatter(X_2d[label==l, 0], X_2d[label==l, 1], label=l)

    plt.xlim(X_2d[:,0].min() - 1, X_2d[:,0].max() + 1)
    plt.ylim(X_2d[:,1].min() - 1, X_2d[:,1].max() + 1)
    plt.legend()

    plt.savefig("result.png")
   
if __name__ == "__main__":
    main()

result.png

　んー、こんなもんですかね。以前の結果も再掲します。

以前の結果（再掲）

　なぜか色の関係が逆転している気もしますがそれはともかく、基本的な構造はどちらでもつかめていると思います。

違いについて

　カーネルPCAだとtransformを呼べるので、データを追加して元の空間でどの辺に来るのかを予想する、といった用途には便利です。

　たとえば、を追加してみます。

import Levenshtein
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import KernelPCA

def main():
    # データ（自作）
    X = ["hogehoge", "hogghheg", "fogefoge", "hagbhcgd", "hogeratt",
         "hohohoho", "fugefuge", "hokehoke", "hogehope", "kogekoge",
         "fugafuga", "fugafuge", "fufufufu", "faggufaa", "fuuuuuuu",
         "fhunfhun", "ufagaguf", "agufaguf", "fogafoga", "fafafaoa"]
    label = np.array(["hoge"]*10 + ["fuga"]*10, dtype=object)

    # 類似度行列の作成（levenshtein距離）
    A = np.zeros((len(X), len(X)))
    for i in range(len(X)):
        for j in range(i + 1, len(X)):
            d = Levenshtein.distance(X[i], X[j])
            A[i,j] = A[j,i] = d
    A = -A  # 距離->類似度変換のためマイナスをかける

    # kernel pcaによる変換
    kpca = KernelPCA(n_components=2, kernel="precomputed")
    X_2d = kpca.fit_transform(A)

    # plot
    for (x, y), l in zip(X_2d, X):
        plt.text(x, y, l)
    
    for l in set(label):
        plt.scatter(X_2d[label==l, 0], X_2d[label==l, 1], label=l)

    # hogefugaを写像（追加）
    data = "hogefuga"
    A = np.zeros((1, len(X)))
    for i in range(len(X)):
        d = Levenshtein.distance(data, X[i])
        A[0,i] = d
    A = -A  # 距離->類似度変換のためマイナスをかける
    hg_X = kpca.transform(A)
    plt.text(hg_X[0,0], hg_X[0,1], data, color="r")

    # 見た目調整
    plt.xlim(X_2d[:,0].min() - 1, X_2d[:,0].max() + 1)
    plt.ylim(X_2d[:,1].min() - 1, X_2d[:,1].max() + 1)
    plt.legend()

    plt.savefig("result2.png")
   
if __name__ == "__main__":
    main()

result2.png

　このようにhogefugaがどのあたりに位置するのかを、既存の点を再計算せずに写像して見てみることができます。MDSだと全体を再計算しないとできないはずなので、これができることがカーネルPCAの明確なメリットの一つだと思います。

　transformに渡す配列のshapeに関しては若干自信がなかったのですが、逆にしたら（つまり.Tをつけたら）

ValueError: Precomputed metric requires shape (n_queries, n_indexed). Got (20, 1) for 20 indexed.

　なる見たことも聞いたこともないようなエラーが出たので、たぶんこれで合っています。

まとめ

　まあ同じようなことができるんだなぁと思いました。

多次元尺度構成法（MDS）で文字列の編集距離を可視化してみる

2019-06-27T23:14:42+09:00

はじめに

　ベクトルとして表現するのが難しいけど、個体間の距離（非類似度）は定義できる……というデータがたまにあります。こういうとき、多次元尺度構成法を使うと可視化がうまくいきます。

　ということで、編集距離を可視化してみようと思います。

データ

　hogehogeとfugafugaを適当に変えた文字列を10個ずつ考えます。こんなのです。

X = ["hogehoge", "hogghheg", "fogefoge", "hagbhcgd", "hogeratt",
     "hohohoho", "fugefuge", "hokehoke", "hogehope", "kogekoge",
     "fugafuga", "fugafuge", "fufufufu", "faggufaa", "fuuuuuuu",
     "fhunfhun", "ufagaguf", "agufaguf", "fogafoga", "fafafaoa"]
label = np.array(["hoge"]*10 + ["fuga"]*10, dtype=object)

　私のセンスが疑われそうな気もしますが、ご容赦ください。

距離行列の計算

　距離（非類似度）の行列は自分で計算する必要があります。今回は編集距離なので、Levenshteinパッケージを使います。

A = np.ones((len(X), len(X)))
for i in range(len(X)):
    for j in range(i + 1, len(X)):
        d = Levenshtein.distance(X[i], X[j])
        A[i,j] = A[j,i] = d

　これも特に難しいことはありません。

コード全体

　importから可視化のプロットなども含めたコードの全体像がこちらです。

import Levenshtein
import numpy as np
import matplotlib.pyplot as plt
from sklearn.manifold import MDS

def main():
    X = ["hogehoge", "hogghheg", "fogefoge", "hagbhcgd", "hogeratt",
         "hohohoho", "fugefuge", "hokehoke", "hogehope", "kogekoge",
         "fugafuga", "fugafuge", "fufufufu", "faggufaa", "fuuuuuuu",
         "fhunfhun", "ufagaguf", "agufaguf", "fogafoga", "fafafaoa"]
    label = np.array(["hoge"]*10 + ["fuga"]*10, dtype=object)

    A = np.zeros((len(X), len(X)))
    for i in range(len(X)):
        for j in range(i + 1, len(X)):
            d = Levenshtein.distance(X[i], X[j])
            A[i,j] = A[j,i] = d

    mds = MDS(n_components=2, dissimilarity="precomputed")
    X_2d = mds.fit_transform(A)
    for (x, y), l in zip(X_2d, X):
        plt.text(x, y, l)
    
    for l in set(label):
        plt.scatter(X_2d[label==l, 0], X_2d[label==l, 1], label=l)

    plt.xlim(X_2d[:,0].min() - 1, X_2d[:,0].max() + 1)
    plt.ylim(X_2d[:,1].min() - 1, X_2d[:,1].max() + 1)
    plt.legend()

    plt.savefig("result.png")
   
if __name__ == "__main__":
    main()

　可視化では文字列のプロットと、散布図をやっています。これでわかりやすく見えるはずです。

実行結果

　結果を以下に示します。

result.png

　似ているものが近くに来るように配置されることがわかります。だいたいきれいに2グループに分かれています。

まとめ

　このように多次元尺度構成法（MDS）が使えるときがあるので、名前だけは覚えておいても損はしません。こういうものもある、ということを知っておきましょう。

追記

　書き上げてから「同じことをやっている人がいるかなぁ」と思って「多次元尺度構成法編集距離」で検索したら、いくつか出てきました。

https://www.jstage.jst.go.jp/article/iieej/38/5/38_5_634/_pdf/-char/ja
　これは画像処理かな。専門外なのでよくわかってはいません。

主座標分析について簡単に紹介するよ！ - ほくそ笑む
　編集距離も使える、という説明だけ。

http://db-event.jpn.org/deim2012/proceedings/final-pdf/c10-1.pdf
　学会の抄録ですかね。

　そんなに目ぼしいものはないといえばない状況です。今回は単純なデータで実験したのでいい感じの結果になりましたが、実用的にはなかなか厳しいのかもしれません。

　あと、MDSでは距離の公理という面倒くさいものが絡んでくることがあります。これを満たさないものは普通のMDSでは扱えないので、「非計量多次元尺度構成法（nMDS）」で取り扱え、ということが言われていたりします。レーベンシュタイン距離は一応大丈夫らしいですが、（研究やビジネスなどで）ちゃんと使いたい場合はこの辺りにも気を配るべきでしょう。

追記2

　同じデータに対して、SVMで分類もやってみました。

scikit-learnのSVMを自分で計算したカーネルで使う - 静かなる名辞

　カーネルPCAでも見てみました。

カーネルPCAで文字列の編集距離を可視化してみる - 静かなる名辞

sklearnとmatplotlibでiris（3クラス）の予測確率を可視化した話

2019-06-22T23:51:33+09:00

はじめに

　よく分類器の性質などを把握するために、2次元で可視化している図があります。

　特に予測確率なんかを平面的に出せるとかっこいいですよね。つまり、こういうのです。

Classifier comparison — scikit-learn 0.21.3 documentation

以前の記事より

君はKNN（k nearest neighbor）の本当のすごさを知らない - 静かなる名辞

　ただ、これが素直にできるのは2クラス分類までで、3クラス分類だと下のような図にしかなりません。

以前の記事より

【python】高次元の分離境界をなんとか2次元で見る - 静かなる名辞

　ということでずっと諦めていたのですが、ふと思いました。

「RGBに各クラスの予測確率あてればできるじゃん」

　簡単にできると思ったら思いの外手間取ったので、備忘録として書いておきます。

まずやる

　とりあえずirisを二次元でプロットします。この辺は定石どおりにやるだけです。

import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    X = pca.fit_transform(iris.data)
    ax = plt.subplot()
    ax.scatter(X[:,0], X[:,1], c=iris.target, cmap="brg")
    plt.savefig("fig1.png")

if __name__ == "__main__":
    main()

fig1.png

　kNNを学習させて、まずは普通に分離境界を描きます。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.neighbors import KNeighborsClassifier

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    X = pca.fit_transform(iris.data)
    ax = plt.subplot()
    ax.scatter(X[:,0], X[:,1], c=iris.target, cmap="brg")
    
    clf = KNeighborsClassifier()
    clf.fit(X, iris.target)
    XX, YY = np.meshgrid(np.arange(-5, 5, 0.025),
                         np.arange(-2, 2, 0.025))
    Z = clf.predict(np.stack([XX.ravel(), YY.ravel()], axis=1))
    ZZ = Z.reshape(XX.shape)
    ax.pcolormesh(XX, YY, ZZ, alpha=0.05, cmap="brg", shading="gouraud")

    plt.savefig("fig2.png")

if __name__ == "__main__":
    main()

　参考：matplotlibのpcolormeshでalphaを小さくすると網目が出てくる対策 - 静かなる名辞

fig2.png

　さ、次はpredict_probaを呼ぶ訳ですが……pcolormeshとかこの辺の関数にはRGBのデータは渡せません。

matplotlib.pyplot.pcolormesh — Matplotlib 3.1.0 documentation

　しばし思案したあと、imshowならできると思いました。

　なにも考えずに書くと下のようなコードになります。

（これは動かないので注意してください）

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.neighbors import KNeighborsClassifier

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    X = pca.fit_transform(iris.data)
    ax = plt.subplot()
    ax.scatter(X[:,0], X[:,1], c=iris.target, cmap="brg")
    
    clf = KNeighborsClassifier()
    clf.fit(X, iris.target)
    XX, YY = np.meshgrid(np.arange(-5, 5, 0.025),
                         np.arange(-2, 2, 0.025))
    Z = clf.predict_proba(np.stack([XX.ravel(), YY.ravel()], axis=1))
    ZZ = Z.reshape(XX.shape + (3, ))
    ax.imshow(ZZ, alpha=0.2)

    plt.savefig("fig3.png")

if __name__ == "__main__":
    main()

fig3.png

　なにこれ？

　ああ、縮尺を合わせないといけないんですね。aspect, extentという引数でできそうです。

matplotlib.pyplot.imshow — Matplotlib 3.1.0 documentation

（これもちゃんと動かないので注意してください）

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.neighbors import KNeighborsClassifier

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    X = pca.fit_transform(iris.data)
    ax = plt.subplot()
    ax.set_xlim((-5, 5))
    ax.set_ylim((-2, 2))
    ax.scatter(X[:,0], X[:,1], c=iris.target, cmap="brg")
    
    clf = KNeighborsClassifier()
    clf.fit(X, iris.target)
    XX, YY = np.meshgrid(np.arange(-5, 5, 0.025),
                         np.arange(-2, 2, 0.025))
    Z = clf.predict_proba(np.stack([XX.ravel(), YY.ravel()], axis=1))
    ZZ = Z.reshape(XX.shape + (3, ))
    ax.imshow(ZZ, alpha=0.2,
              aspect="auto", extent=(-5, 5, -2, 2))

    plt.savefig("fig4.png")

if __name__ == "__main__":
    main()

　まず、先にax.set_xlimとax.set_ylimで図の範囲を指定し、そこにextentをあわせるようにしています。aspectはドキュメントを見た感じだとautoが無難そうに思います。

fig4.png

　どう見ても上下反転しているので、ZZを上下反転します。ついでに、マーカーの色を揃えることにします。

これが動くコードです

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.neighbors import KNeighborsClassifier

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    X = pca.fit_transform(iris.data)
    ax = plt.subplot()
    ax.set_xlim((-5, 5))
    ax.set_ylim((-2, 2))
    cm = ListedColormap(["b", "g", "r"])
    ax.scatter(X[:,0], X[:,1], c=iris.target, cmap=cm)
    
    clf = KNeighborsClassifier()
    clf.fit(X, iris.target)
    XX, YY = np.meshgrid(np.arange(-5, 5, 0.025),
                         np.arange(-2, 2, 0.025))
    Z = clf.predict_proba(np.stack([XX.ravel(), YY.ravel()], axis=1))
    ZZ = np.flip(Z.reshape(XX.shape + (3, )), axis=1)
    ax.imshow(ZZ, alpha=0.2,
              aspect="auto", extent=(-5, 5, -2, 2))

    plt.savefig("fig5.png")

if __name__ == "__main__":
    main()

fig5.png

　だいたい不満のない結果になりました。ここまで長かった。

他の分類器も試す

　せっかくなのでいろいろやってみます。SVM, ロジスティック回帰, ランダムフォレストを追加してやってみましょう。
　

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.neighbors import KNeighborsClassifier
from sklearn.svm import SVC
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier

def main():
    iris = load_iris()
    pca = PCA(n_components=2)
    X = pca.fit_transform(iris.data)

    fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(9, 9))
    knn = KNeighborsClassifier()
    svm = SVC(probability=True)
    lr = LogisticRegression()
    rfc = RandomForestClassifier(n_estimators=100)

    cm = ListedColormap(["b", "g", "r"])    
    XX, YY = np.meshgrid(np.arange(-5, 5, 0.025),
                         np.arange(-2, 2, 0.025))

    for ax, clf in zip(axes.ravel(), [knn, svm, lr, rfc]):
        ax.set_xlim((-5, 5))
        ax.set_ylim((-2, 2))
        ax.scatter(X[:,0], X[:,1], c=iris.target, cmap=cm)
    
        clf.fit(X, iris.target)
        Z = clf.predict_proba(np.stack([XX.ravel(), YY.ravel()], axis=1))
        ZZ = np.flip(Z.reshape(XX.shape + (3, )), axis=1)
        ax.imshow(ZZ, alpha=0.2,
                  aspect="auto", extent=(-5, 5, -2, 2))
        ax.set_title(clf.__class__.__name__)

    plt.tight_layout()
    plt.savefig("fig6.png")

if __name__ == "__main__":
    main()

fig6.png

　こんなもんか、という感じ。

まとめ

　やればできることはわかりました。もう少しかっこいい図にするには、さらなる工夫が要るのかもしれません。

【python】高次元の分離境界をなんとか2次元で見る

2019-05-24T04:01:31+09:00

はじめに

　分類器の特性を把握するために2次元データで分離境界を見るということが行われがちですが、高次元空間における分離器の特性を正確に表している訳ではありません。

　ということがずっと気になっていたので、なんとか高次元空間で分類させて2次元で見る方法を考えます。

方法

　PCAで2次元に落とせれば、線形変換で逆変換もできるので、それでやります。当然ながら情報は落ちますし、2次元でもなんとか見える程度のデータしか扱えませんが、妥協します。

　sklearnならinverse_transformという便利なメソッドがあるので、簡単です。

　というあたりまで考えた上で、こんなコードを書きました。

show_hyperplane.py

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

def show_hyperplane(dataset, clf, filename):
    pca = PCA(n_components=2)
    X = pca.fit_transform(dataset.data)
    plt.scatter(X[:,0], X[:,1], c=dataset.target)

    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.01),
                         np.arange(y_min, y_max, 0.01))

    clf.fit(dataset.data, dataset.target)
    Z = clf.predict(
        pca.inverse_transform(np.c_[xx.ravel(), yy.ravel()]))
    plt.pcolormesh(xx, yy, Z.reshape(xx.shape),
                   alpha=0.03, shading="gouraud")
    plt.savefig(filename)

　汎用的に作ったので、これでいろいろなものを見てみようという算段です。

実験

　まずirisとSVM。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_iris
from sklearn.svm import SVC

iris = load_iris()
svm = SVC(C=50, gamma="scale")    
show_hyperplane(iris, svm, "iris_svm.png")

iris_svm.png

　特に興味深い知見は得られませんでした。

　次、irisとランダムフォレスト。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_iris
from sklearn.ensemble import RandomForestClassifier

iris = load_iris()
rfc = RandomForestClassifier(n_estimators=500, n_jobs=-1)    
show_hyperplane(iris, rfc, "iris_rf.png")

iris_rf.png

　ランダムフォレストで斜めの分離超平面の図を出したサイトはここくらいしかないのでは？　だからどうしたって話ですが。

　簡単なのでAdaBoostも試します。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_iris
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import AdaBoostClassifier

iris = load_iris()
ada = AdaBoostClassifier(
    base_estimator=DecisionTreeClassifier(max_depth=4),
    n_estimators=200)    
show_hyperplane(iris, ada, "iris_ada.png")

iris_ada.png

　面白いんですが、性能はいまいち悪そう。

　ちなみに、base_estimatorのパラメータでコロコロ結果が変わります。パラメータ設定については、以下の2記事を参照してください。

【python】sklearnのAdaBoostをデフォルトパラメータで使ってはいけない - 静かなる名辞
 AdaBoostとRandomForestの比較 - 静かなる名辞

　ただの決定木もやっておきましょう。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_iris
from sklearn.tree import DecisionTreeClassifier

iris = load_iris()
dtc = DecisionTreeClassifier()
show_hyperplane(iris, dtc, "iris_tree.png")

iris_tree.png

　つまらない。

　さて、irisは飽きてきたのでdigitsで同じことをやります。こちらは何しろ元が64次元で、2次元に落とすとかなり重なり合うので、カオスな結果になってくれそうです。

　が、その前にshow_hyperplane.pyをいじります。元のままだといろいろうまくいかなかったからです。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

def show_hyperplane(dataset, clf, filename):
    pca = PCA(n_components=2)
    X = pca.fit_transform(dataset.data)
    plt.scatter(X[:,0], X[:,1], s=5, c=dataset.target)

    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.3),
                         np.arange(y_min, y_max, 0.3))

    clf.fit(dataset.data, dataset.target)
    Z = clf.predict(
        pca.inverse_transform(np.c_[xx.ravel(), yy.ravel()]))
    plt.pcolormesh(xx, yy, Z.reshape(xx.shape),
                   alpha=0.05, shading="gouraud")
    plt.savefig(filename)

　よし、やろう。

　まずSVM。今回からついでに学習データに対するスコアを見ます。コメントで記します。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_digits
from sklearn.svm import SVC

digits = load_digits()
svm = SVC(C=0.1, gamma="scale")    
score = svm.fit(
    digits.data, digits.target).score(
        digits.data, digits.target)
print(score) # => 0.9744017807456873
show_hyperplane(digits, svm, "digits_svm.png")

digits_svm.png

　あたりまえですが、64→2次元で情報落ちしているので、こんなふうにしか見えません。それでも、後々出てくるやつに比べればまともな方です。

　次。ランダムフォレスト。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_digits
from sklearn.ensemble import RandomForestClassifier

digits = load_digits()
rfc = RandomForestClassifier(n_estimators=500, n_jobs=-1)    
score = rfc.fit(
    digits.data, digits.target).score(
        digits.data, digits.target)
print(score) # => 1.0
show_hyperplane(digits, rfc, "digits_rfc.png")

digits_rfc.png

　これ面白いですね。ところどころ凹凸がありますが、それでもぱっと見SVMと同じくらい滑らかな分離超平面に見えます。高次元データほど強いというのもわかる気がします。

　アダブースト。

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_digits
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import AdaBoostClassifier

digits = load_digits()
ada = AdaBoostClassifier(
    base_estimator=DecisionTreeClassifier(max_depth=3),
    n_estimators=200)    
score = ada.fit(
    digits.data, digits.target).score(
        digits.data, digits.target)
print(score) # 0.9660545353366722
show_hyperplane(digits, ada, "digits_ada.png")

digits_ada.png

　大丈夫なんかこれ。決定木のアダブーストはランダムフォレストと比べて個人的にいまいち信頼していないのですが、こういうの見るとその思いが強まります。

　決定木もやりましょう。irisではつまらなかったけど、こちらではどうなるでしょうか。
　

from show_hyperplane import show_hyperplane
from sklearn.datasets import load_digits
from sklearn.tree import DecisionTreeClassifier

digits = load_digits()
dtc = DecisionTreeClassifier()
score = dtc.fit(
    digits.data, digits.target).score(
        digits.data, digits.target)
print(score) # 1.0
show_hyperplane(digits, dtc, "digits_tree.png")

digits_tree.png

　あははー、なにこれカオス。デフォルトのまま高次元で使うな、ということですね。

まとめ

　元が64次元くらいでもだいぶ情報落ちするので、本当の高次元データでも使えるかというと微妙なのですが、それでもなんとなく傾向はつかめますし、面白かったです。

　SVMとランダムフォレストはどっちも優秀ですね。

【python】PCAと非負値行列因子分解のバイプロットを見比べる

2019-05-14T20:48:35+09:00

はじめに

　非負値行列因子分解は負の値が出現しないような行列に対して行える分解で、主成分分析とか因子分析に似ています。

　参考：
非負値行列因子分解（NMF）をふわっと理解する - Qiita

　上の記事によると、いいところとしては、

非負なので現実のデータに向く
非負なので解釈が楽
さらにスパースになる

　というあたりらしい。

　なので、PCAと比べます。sklearnを使います。

比較実験

　irisでやります。なんとかの一つ覚えです。

import numpy as np
import matplotlib.pyplot as plt

from sklearn.decomposition import PCA, NMF
from sklearn.datasets import load_iris

def main():
    iris = load_iris()

    pca = PCA(n_components=2)
    nmf = NMF(n_components=2)

    fig, axes = plt.subplots(nrows=1, ncols=2)

    for i, mname, method in zip([0,1], ["PCA", "NMF"], [pca, nmf]):
        X_2d = method.fit_transform(iris.data)

        # title
        axes[i].set_title("{} {}".format("iris", mname))

        # scatter
        axes[i].scatter(X_2d[:,0], X_2d[:,1], c=iris.target)

        # arrows
        pc0 = method.components_[0]
        pc1 = method.components_[1]

        pc0 = pc0 * (np.abs(X_2d[:,0]).max() / np.abs(pc0).max()) * 0.8 
        pc1 = pc1 * (np.abs(X_2d[:,1]).max() / np.abs(pc1).max()) * 0.8

        for j in range(pc0.shape[0]):
            axes[i].arrow(
                0, 0, pc0[j], pc1[j], color='r')
            axes[i].text(
                pc0[j]*1.1, pc1[j]*1.1, iris.feature_names[j], color='r')

    plt.show()

if __name__ == "__main__":
    main()

結果

　こうして見るとそんなに違いませんが、原点の右上だけが図に含まれるのが見ての通り特徴です。相違点としては

スケールの違いなどがわかるような気がする
PCAでは重なっているpetal width(cm)とpetal length(cm)の違いが出ている

　などがあるでしょうか。また、NMFではpetal width(cm)のy方向成分はゼロのようです。

メリット

　上の図を見る限りでは、別にどっちでも大差はなさそうだし、PCAの方が慣れているので意味的な解釈もしやすい気がします。

　非負であることが要請されるような特殊ケース以外は別にPCAでもこまらないという気もするのですが、実際のところどうなんでしょうね。

まとめ

　とにかく使えます。

【python】sklearnのFeatureAgglomerationを使ってみる

2018-12-10T03:56:41+09:00

はじめに

　FeatureAgglomerationは階層的クラスタリングを用いた教師なし次元削減のモデルです。特徴量に対して階層的クラスタリングを行い（つまり通常のサンプルに対するクラスタリングと縦横の向きが入れ替わる）、似ている特徴量同士をマージします。マージの方法はデフォルトでは平均のようです。

　使用例をあまり見かけませんが、直感的な次元削減方法なので何かしらの役に立つかもしれないと思って使ってみました。

sklearn.cluster.FeatureAgglomeration — scikit-learn 0.20.1 documentation

使い方

　パラメータは以下の通り。

class sklearn.cluster.FeatureAgglomeration(
    n_clusters=2, affinity=’euclidean’, memory=None, connectivity=None, 
    compute_full_tree=’auto’, linkage=’ward’, pooling_func=<function mean>)

　色々いじれるように見えますが、主要パラメータは2つだけです。

n_clusters

　PCAでいうところのn_componentsです。変換先の次元数を表します。

pooling_func

　似ている特徴量をマージする方法。callableが渡せます。何もしなければ平均が使われるので、平均より気の利いた方法を思いつく人以外はそのままで大丈夫です。

　あとは階層的クラスタリングのオプションが色々あります。それはそれで大切なものだと思いますが、今回は無視することにします。

実験

　もう何番煎じかわかりませんが、irisの2次元写像で試します。

import matplotlib.pyplot as plt

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import FeatureAgglomeration

def main():
    iris = load_iris()

    pca = PCA(n_components=4)
    ss = StandardScaler()
    agg = FeatureAgglomeration(n_clusters=2)

    pca_X = pca.fit_transform(iris.data)
    agg_X = agg.fit_transform(
        ss.fit_transform(iris.data))

    print(pca.components_)
    print(agg.labels_)

    fig, axes = plt.subplots(nrows=1, ncols=2)
    axes[0].scatter(pca_X[:,0], pca_X[:,1], c=iris.target)
    axes[0].set_title("PCA")
    axes[1].scatter(agg_X[:,0], agg_X[:,1], c=iris.target)
    axes[1].set_title("FeatureAgglomeration\n{}".format(agg.labels_))
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

　動作原理、目的と用途を考えると、事前にスケーリングしておいた方が恐らく無難です。

　printされた出力。

[[ 0.36138659 -0.08452251  0.85667061  0.3582892 ]
 [ 0.65658877  0.73016143 -0.17337266 -0.07548102]]
[0 1 0 0]

　FeatureAgglomerationは圧倒的に結果の解釈性が良いことがわかります。写像先の0次元目は元の0,2,3次元目の平均で*1、写像先の1次元目は元の1次元目ですね。こういうのはシチュエーション次第ですが、ちょっと嬉しいかもしれません。

　出力される画像。

プロットの結果

　概ねPCAと同等に使えています。うまく言葉で表現はできませんが、FeatureAgglomerationの方はなんとなくギザギザ感？みたいなものがあります。平均するとそうなる、というのがなんとなくわかる気もするし、わからない気もする。

考察

　結果の解釈性が良いのと、まがりなりにすべての特徴量の情報が結果に反映されるので、PCAより使いやすいシチュエーションはあると思います。分類前の次元削減とかで使ったときの性能とかは今回検討していませんが、たぶんそんなに良いということはないはず。

　あとドキュメントをあさっていたら、こんなページがあったので、

Feature agglomeration — scikit-learn 0.20.1 documentation

　真似してPCAでも同じものを出してみたら（コードはほとんど書き換えていないので省略。agglo = の行で代入するモデルをコメントアウトで切り替えて、あとlabels_の出力を外しただけです）、やっぱりFeatureAgglomerationはヘボかった（低次元で元の情報を保持することに関しては性能が低かった）です。

　10次元に落として元の情報をどこまで復元できるかという実験。

PCA

FeatureAgglomeration

　まあ、これは仕方ないか。

まとめ

　とにかく結果の解釈性の良さを活かしたい、とか、なにか特別な理由があって使う分には良いと思います。

*1:厳密にはどれか2つが先に平均されて、更に残りと平均されるはず。つまり3つの比重が違う順番はチェックしていないのでわかりませんが、children_属性をちゃんと読み取ればわかると思います

【python】sklearnのSparsePCAを使ってみる

2018-11-17T22:30:03+09:00

はじめに

　SparsePCAというものがあることを知ったので、使ってみようと思います。

SparsePCAとは？

　その名の通り、スパースな主成分分析です。スパースな主成分ベクトルを推定します。

Sparse PCA - Wikipedia

　原理などは理解しないで、カジュアルに使えるかどうか試してみるだけという趣旨です。なので「どうやって動いているの？」という質問には答えられません。許してください。

sklearnの実装

　きっちり存在しています（存在しなかったらこんな記事は書きませんが）。

sklearn.decomposition.SparsePCA — scikit-learn 0.20.1 documentation

　主要なパラメータとしては、以下のものがあります。

n_components

　PCAのと同じです。

alpha

　スパースPCAのキモで、L1正則化の強さを調整できます。

ridge_alpha

　こちらはtransformの際に使われるリッジ回帰（L2正則化）の正則化パラメータです。なんでリッジを使うのかは、実のところよくわかりません。

max_iter

　このパラメータがあるということは、最適化とか勾配法的なもので推定するのだな、というくらいに思っておきます。

normalize_components

　主成分ベクトルのノルムを1にするかどうか。Trueにしておくと良いと思います。

　結果に大きな影響を及ぼすのは上くらいだと思います。他のパラメータについてはドキュメントを参照してください。

実験

　今回はwineデータセットでやってみました。素のPCAでやった場合、alphaを0.5と5にした場合の結果をバイプロットで示します。

import matplotlib.pyplot as plt
from sklearn.datasets import load_wine
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA, SparsePCA

def biplot(X_2d, components, target, ax):
    r1 = 5
    r2 = 1.01
    for i, coef in enumerate(components.T):
        ax.arrow(0, 0, coef[0]*r1, coef[1]*r1, color='r')    
        ax.text(coef[0]*r1*r2, coef[1]*r1*r2, i, color='b', fontsize=8)

    ax.scatter(X_2d[:,0], X_2d[:,1], c=target, cmap="rainbow")

def main():
    wine = load_wine()
    ss = StandardScaler()
    X = ss.fit_transform(wine.data)

    pca = PCA(n_components=2)
    spca = SparsePCA(n_components=2,
                     max_iter=3000,
                     n_jobs=-1,
                     normalize_components=True)
    
    fig, axes = plt.subplots(figsize=(12, 6), nrows=1, ncols=3)

    X_2d = pca.fit_transform(X)
    biplot(X_2d, pca.components_, wine.target, axes[0])
    axes[0].set_title("PCA")

    for i,alpha in zip([1, 2], [0.5, 5]):
        spca.set_params(alpha=alpha)
        X_2d = spca.fit_transform(X)
        biplot(X_2d, spca.components_, wine.target, axes[i])
        axes[i].set_title("SPCA alpha={:.2f}".format(alpha))
    plt.savefig("result.png")

    # 図と突き合わせて確認するために特徴量の名前を出力しておく
    for i, name in enumerate(wine.feature_names):
        print(i, name)

if __name__ == "__main__":
    main()

　max_iterをきもち高めにしましたが、結果は数秒程度で出ました。

result.png

0 alcohol
1 malic_acid
2 ash
3 alcalinity_of_ash
4 magnesium
5 total_phenols
6 flavanoids
7 nonflavanoid_phenols
8 proanthocyanins
9 color_intensity
10 hue
11 od280/od315_of_diluted_wines
12 proline

　とりあえず、PCAの結果とSparsePCAの結果で左右が反転しているのに注意。

　あとは見ての通りで、alpha=0.5で一部の係数が主成分にべたっと張り付くようになり、alpha=5では大半の係数が主成分に張り付いています。これがSparsePCAの効果で、結果の解釈が容易になるということらしいです（この次元数だとあまり威力はありませんが、高次元では活躍しそうです）。

　ワインにはあまり詳しくないので、今回は結果を細かく解釈することはしませんが……。

まとめ

　使えることがわかりました。

【python】複数の特徴をまとめるFeatureUnion

2018-05-15T14:41:24+09:00

　単一の入力データから、複数の処理方法で幾つもの異なる特徴量が得られる・・・というシチュエーションがある。

　この場合、「どれが最善か」という観点でどれか一つを選ぶこともできるけど、そうすると他の特徴量の情報は捨ててしまうことになる。総合的な性能では他に一歩譲るが、有用な情報が含まれている特徴量がある・・・というような場合は、ちょっと困る。

　こういう状況で役に立つのがFeatureUnion。特徴抽出や次元削減などのモデルを複数まとめることができる。

　結果はConcatenateされる。Concatenateというのがわかりづらい人もいると思うけど、たとえば手法1で10次元、手法2で20次元の特徴量ベクトルが得られたら、これをそのまま横に繋げて30次元のベクトルとして扱うということ。

sklearn.pipeline.FeatureUnion — scikit-learn 0.20.1 documentation

　ちなみに、こいつはsklearn.pipeline以下に存在する。Pipelineの兄弟みたいな扱い。引数の渡し方とかもほとんど同じである。

　簡単に試してみよう。digitsの分類を行うことにする。PCA+GaussianNB, LDA+GNB, FeatureUnion(PCA, LDA)+GNBの3パターンでスコアを見比べる。

import warnings
warnings.filterwarnings('ignore')

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.naive_bayes import GaussianNB
from sklearn.pipeline import Pipeline, FeatureUnion
from sklearn.model_selection import cross_validate, StratifiedKFold

def main():
    digits = load_digits()
    
    pca = PCA(n_components=30)
    lda = LDA()
    gnb = GaussianNB()
    
    pca_gnb = Pipeline([("pca", pca), ("gnb", gnb)])
    lda_gnb = Pipeline([("lda", lda), ("gnb", gnb)])
    pca_lda_gnb = Pipeline([("reduction", FeatureUnion([("pca", pca),
                                                        ("lda", lda)])),
                            ("gnb", gnb)])

    scoring = {"p": "precision_macro",
               "r": "recall_macro",
               "f":"f1_macro"}

    for name, model in zip(["pca_gnb", "lda_gnb", "pca_lda_gnb"], 
                           [pca_gnb, lda_gnb, pca_lda_gnb]):

        skf = StratifiedKFold(shuffle=True, random_state=0)
        scores = cross_validate(model, digits.data, digits.target,
                                cv=skf, scoring=scoring)
        
        p = scores["test_p"].mean()
        r = scores["test_r"].mean()
        f = scores["test_f"].mean()
        print(name)
        print("precision:{0:.3f} recall:{1:.3f} f1:{2:.3f}".format(p,r,f))

if __name__ == "__main__":
    main()

　結果は、

pca_gnb
precision:0.947 recall:0.944 f1:0.945
lda_gnb
precision:0.955 recall:0.953 f1:0.953
pca_lda_gnb
precision:0.959 recall:0.957 f1:0.957

　ちょっと微妙だけど、誤差ではないみたい。このように比較的手軽に性能を改善できることがわかる（効くかどうかはケースバイケースだけど）。

【python】SOMのライブラリSomocluはかなりおすすめ

2018-04-07T16:12:49+09:00

　SOM(Self-organizing maps：自己組織化写像)は割と古めの、データの可視化手法です（それ以外にも使えると思いますが）。

　今回はpythonのSOMライブラリSomocluを使ってみたら、けっこう良かったというネタです。

SOMの概要
ライブラリがない
それでも頑張ってググった
使ってみた
今どきSOMなんか使うの？（蛇足パート）
まとめ

SOMの概要

　昨今は深層学習が流行りですが、SOM、自己組織化写像は敢えて言えば単層学習とでも言うべきでしょうか。平面上だったり立体状（まあ理屈の上では何次元でも定義できる）に並べたニューロンにデータをマッピングします。それ以上の説明はwikipediaとか、ググれば色々出てくるページを読んでください。

wikipedia

自己組織化写像 - Wikipedia

九州工業大学大学院の先生が書いた読みやすかったページ

http://www.brain.kyutech.ac.jp/~furukawa/data/SOMtext.pdf

わかりやすい解説

子供でもわかる「自己組織化マップ」

ライブラリがない

　SOM、けっこう面白い性質があるみたいなのて使ってみたいのですが、ググってみるとpythonで使えそうなライブラリがとにかくあまり出てきません。

SOMPY

　申し訳ないけど、ちょっと使いづらかった。というかインストールしても挙動が変な感じだった。
GitHub - sevamoo/SOMPY: A Python Library for Self Organizing Map (SOM)

sompy

　日本人の方が実装されたようです。率直に言って「作ってみた」レベルで、実用にはどうかという感じ
自己組織化マップ(SOM)のPythonライブラリsompyを公開しました - 俺とプログラミング

PyMVPA

　多変量解析のためのそれなりに大きいライブラリで、SOMも実装されている。これが使えればよかったのだと思うが、python2系のサポートしかないので没・・・。
Self-organizing Maps — PyMVPA 2.6.1.dev1 documentation

　他にも色々あったのですが、割愛。古い手法なので、敢えて作ろうという人がいないのかな・・・。

　というか、SOMでググると「実装してみた」系の記事はたくさん出てくるのに、まともに使えるライブラリは出てこないというの、かなり異常というか残念というか・・・。

それでも頑張ってググった

　Somocluというのを見つけました。

Introduction — Somoclu 1.7.5 documentation

　ウリの部分を適当に訳したり訳さなかったりしつつ抜粋

OpenMPとCUDAがサポートされていてGPUでも計算できる
当然マルチプラットフォームでLinux, macOS, and Windowsでサポートされている
「Planar and toroid maps」平面とドーナツみたいな形のSOM両方が作れる
「Rectangular and hexagonal grids」四角と六角形がいける
「Gaussian or bubble neighborhood functions」近傍の計算を効率化する系のがある
「Visualization of maps, including those that were trained outside of Python.」
マップの初期化にはPCAが使える

　すごく良さそう。あと、pythonに依存しないツールでコマンドラインから直接コマンドで叩けます。pythonバインディングもあるよ、という位置づけ。真剣に開発されてる感じです。

使ってみた

　とりあえず使ってみました。SOMの可視化結果でよく見るU-matrixという奴を出します。以下のコードで動きました。

# coding: UTF-8
import numpy as np

from somoclu import Somoclu
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA

def main():
    # データを読み込む
    dataset = load_iris()
    X = dataset.data
    y = dataset.target
   
    # SOMに入れる前にPCAして計算コスト削減を測る（iris程度では無駄） 
    pca = PCA(n_components=0.95) 
    X = pca.fit_transform(X)

    # SOMの定義
    n_rows = 16
    n_cols = 24
    som = Somoclu(n_rows=n_rows, n_columns=n_cols,
                  initialization="pca", verbose=2)

    # 学習
    som.train(data=X, epochs=1000)

    # U-matrixをファイル出力
    som.view_umatrix(labels=y, bestmatches=True,
                     filename="umatrix.png")

if __name__ == "__main__":
    main()

　説明不要な感じ。コードも直感的だし、特に不満がないです。

　こんな画像が出てきます。

U-matrix

　この画像の見方は色の濃淡が重要で、色の明るい部分は相対的に縮尺が縮んでおり、逆に暗い部分は縮尺が相対的に大きい訳です。PCAで可視化した結果を参考に貼っておきます。

PCAによるirisの可視化結果

　紫がラベル0に、緑と黄色が1と2に対応している訳です。SOMを使うと、このようにデータの構造を捉えることができます。

　使いやすいし動作もまともだし、Somocluは素晴らしいライブラリです。SOMが必要になったら積極的に使っていきたいところ。

今どきSOMなんか使うの？（蛇足パート）

　t-SNEみたいなよくできた手法があるのに今更SOM？　と思う方もおられるかと思いますが、SOMはSOMでメリットがあると感じています。

　というのは、t-SNEはけっきょくパラメタに依存するし、ミクロな構造を捉えるのは得意でもマクロな構造はどこまで正しいのか？　という問題があるからです。

　例として、digitsを可視化してみます。

# coding: UTF-8
import numpy as np

from sklearn.datasets import load_digits
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
from somoclu import Somoclu
import matplotlib.pyplot as plt

def main():
    print("loading data")
    digits = load_digits()
    pca = PCA(n_components=0.95)
    pca_data = pca.fit_transform(digits.data)

    # tsneで可視化
    print("tsne")
    tsne = TSNE()
    X = tsne.fit_transform(pca_data)
    fig, ax = plt.subplots()
    plt.scatter(X[:,0], X[:,1], c=digits.target/10)
    
    i = 0
    for xy, l in zip(X, digits.target):
        if i%8 == 0: # 描画されるtextが多いと汚いので省く
            ax.annotate(l, xy=xy)
        i += 1
    plt.savefig("tsne_digits.png")

    # somで可視化
    print("som")
    # データを適当に省く
    sample_index = np.random.choice(X.shape[0], 400, replace=False)
    sample_X = pca_data[sample_index]
    sample_y = digits.target[sample_index]

    # som
    som = Somoclu(n_rows=30, n_columns=40,
                  initialization="pca")
    som.train(data=sample_X, epochs=1000)
    som.view_umatrix(labels=sample_y, bestmatches=True,
                     filename="som_digits.png")

if __name__ == "__main__":
    main()

t-SNEで可視化したdigits

SOMで可視化したdigits

　一見するとt-SNEは同じラベルごとにまとまっていて綺麗なんですが、形の似ている数字が近くに来るのはむしろSOMの方という気もします。0の周りに5,6,9が来るというのは（数字の形を考えると）妥当そうですね。主観的になってしまいますが、SOMも捨てたものではないという気がします。

まとめ

　SOMとSomocluは良いのでみんな使おう。

【python】sklearnで因子分析を試す

2018-03-31T00:22:11+09:00

　pythonで因子分析をやる人はあまりいないようだが、sklearnにはしっかりモデルが存在している。ついさっき気づいた。

sklearn.decomposition.FactorAnalysis — scikit-learn 0.20.1 documentation

　因子分析自体は前からどんなものなのか興味があり、かといってググるとRだったりSPSSだったりばっかり出てきて辟易していたのだが、sklearnにあると都合が良い。さっそく使ってみよう。

とりあえずirisをプロットする
とりあえずcomponentsを見る
使えることはわかった

とりあえずirisをプロットする

　私だけでも何十回もやってきた、世界中では何万回とやられてきたirisの二次元可視化をやってみる。

　次のようなコードを書いた。

# coding: UTF-8

from copy import deepcopy
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA, FactorAnalysis as FA
from sklearn.pipeline import Pipeline
import matplotlib.pyplot as plt

def decomp_and_plot(dataset, model, file_name):
    X = model.fit_transform(dataset.data)
    plt.figure()
    plt.scatter(X[:,0], X[:,1], c=dataset.target/len(dataset.target_names))
    plt.savefig(file_name)
    
def main():
    iris = load_iris()

    ss = StandardScaler()
    pca = PCA(n_components=2)
    pl = Pipeline([("scaler", ss), ("pca", deepcopy(pca))])
    fa = FA(n_components=2, max_iter=5000)

    decomp_and_plot(iris, pca, "pca_plt.png")
    decomp_and_plot(iris, pl, "spca_plt.png")
    decomp_and_plot(iris, fa, "fa_plt.png")

if __name__ == "__main__":
    main()

　PCA、変数をスケーリングしたPCA（相関行列を使うことと等価）、因子分析でそれぞれplotしてみる。

　結果はこれ。

PCAの結果

PCA(相関行列)の結果

　相関行列はぱっと見いまいち（この絵一枚でダメかどうかは判断できないが）。

因子分析の結果

　うーん、相関行列のとも違うし、なんとも言い難いというか、素人目にはぶっちゃけあんまり良くないように見えるのだが、確率モデルなのでノイズの存在を仮定して見るとこうなるということだろう。

とりあえずcomponentsを見る

　次のようなmain2を作り、実行した。

def main2():
    iris = load_iris()

    print(iris.feature_names)
    print("pca")
    pca = PCA(n_components=2)
    pca.fit(iris.data)
    print(pca.components_)

    print("fa")
    fa = FA(n_components=2, max_iter=5000)
    fa.fit(iris.data)
    print(fa.components_)

　結果

['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)']
pca
[[ 0.36158968 -0.08226889  0.85657211  0.35884393]
 [ 0.65653988  0.72971237 -0.1757674  -0.07470647]]
fa
[[ 0.72577591 -0.17754023  1.75733754  0.73196365]
 [-0.37036948 -0.24060118  0.02793388  0.04121372]]

　プロット結果から予想される通り、両者のcomponentsはよく似通っている。

　これがloadingなのかどうかはぶっちゃけよくわからないのだが（というか1を超えてくる時点でたぶん違うのだろうが）、とりあえずloadingだと思って解釈する。

　第一因子は花弁の長さと幅、がく片の長さに対応しているので花の大きさに対応しているっぽい。花の大きさとがく片の幅はなぜか若干反比例する。

　第二因子は花弁に関する係数が小さいので、がく片の大きさを表す因子と言って良さそうである。

　こんなところか。

使えることはわかった

　だから何？　って言われると、正直答えに窮しますが・・・とにかく使えます。主成分分析で良いじゃんと言われたら何も言い返せません。
　

【python】pythonで主成分分析のバイプロット

2018-03-28T23:13:05+09:00

　バイプロット（Biplot）という主成分分析（PCA）の結果の可視化方法があります。

　すごく大雑把に言うと、PCAによる写像の前の空間の各特徴（軸）が写像先の空間のどこに向いているかを可視化する方法です。

　具体的には、主成分ベクトル（因子負荷量などを使う場合もあります）と散布図を同じ図にplotします。これらを組み合わせることで、元の空間の性質が二次元（もしかしたら3次元）で手に取るようにわかります*1。

　バイプロットはR言語だと簡単に描けるらしいのですが、我らがpythonには（少なくとも一般的なライブラリには）そんな便利なものはありません。ちょっと困るのですが、幸い英語圏にはちらほらやりかたの情報があります。しかし、それはそれでページごとにやってることが違ったりして、（申し訳ないのですが）微妙に信用できなかったりします。

　で、けっきょく自分で書いてみることにしました。なお、参考にしたのはこの辺です。

方針

　まずsklearnの公式ドキュメントをできるだけ良く読み込みます。

sklearn.decomposition.PCA — scikit-learn 0.22.1 documentation

　PCA.components_が固有ベクトルであり、データをセンタリングしてこれと掛けるとPCAの出力が出てくることは前回の記事で確認しました。

　固有ベクトル行列が主成分*元のデータの特徴という形になっているとして、横に見ると負荷量（みたいなもの。本当は対応する固有値のsqrtを掛け算してやらないといけない）に、縦に見ると元の写像先で表現された特徴の軸になります。

　つまり、その軸をプロットするだけです。

　なお、この辺は微妙に議論があるようです。私もこれがどこまで正しい方法なのかは自信が持てません。

　参考：
色々と考えてみる: 文系のための「主成分分析の可視化」（２）

　だけど今回は、データをセンタリングしてPCAを学習させた上で、各軸に対応するone-hot vectorを渡してtransformしたら確かに上に書いた方法通りで上手く行きました（biplotの線の上に載った）。なので、「これで良いんだろう」と勝手に判断しました。どこまで妥当かはよくわからないんですけど。

実装

　こんな感じで書きました。

# coding: UTF-8

from sklearn.datasets import load_iris, load_wine
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

import matplotlib.pyplot as plt

def biplot(dataset, scale=False, arrow_mul=1, text_mul=1.1):
    if scale:
        ss = StandardScaler()
        X = ss.fit_transform(dataset.data)
    else:
        X = dataset.data

    if hasattr(dataset, "feature_names"):
        feature_names = list(dataset.feature_names)
    else:
        feature_names = ["F{0}".format(i)
                         for i in range(dataset.data.shape[1])]

    pca = PCA(n_components=2)
    X = pca.fit_transform(X)

    x_data = X[:,0]
    y_data = X[:,1]

    pc0 = pca.components_[0]
    pc1 = pca.components_[1]

    plt.figure()
    plt.scatter(x_data, y_data,
                c=dataset.target/len(set(dataset.target)),
                marker=".")

    for i in range(pc0.shape[0]):
        plt.arrow(0, 0, 
                  pc0[i]*arrow_mul, pc1[i]*arrow_mul,
                  color='r')
        plt.text(pc0[i]*arrow_mul*text_mul,
                 pc1[i]*arrow_mul*text_mul,
                 feature_names[i],
                 color='r')
    plt.show()

def main():
    iris = load_iris()
    wine = load_wine()

    biplot(iris, arrow_mul=2.5, scale=True)
    biplot(wine, arrow_mul=6, scale=True)

if __name__ == "__main__":
    main()

　今回はsklearnのデータセットを渡す形で関数にまとめました。ま、もしこのコードを流用したい人がいたら、必要なロジックだけ上手く切り出してください。

　結果は、こんな画像が出ます。

irisのバイプロット

wineのバイプロット

　上手く行ってる感じです。

　なお、上のコードでは変数をスケーリングしています（相関行列でPCAするのと等価）。スケーリングしなくてもできますが、やった方が矢印の長さが揃いやすいです（逆に変数のスケールを重視してPCAしたいときは、スケーリングしてはいけない。ケースバイケース）。

まとめ

　これくらい自作しなくても済めば良いのにと思いました。

*1:本当に手に取るようにわかるかはデータと見る人に依存しますが・・・

【python】numpyで主成分分析を実装してみた

2018-03-28T22:21:01+09:00

　numpyでPCA（principal component analysis：主成分分析）を実装してみました。自分の理解を深めるためです。

　sklearnに実装されているものと同じ結果を出すことを目標にしました。最終的には上手く行きました。

概要
実装
結果
まとめ

概要

　主成分分析のアルゴリズムの解説は他に譲ります。これは実装してみた記事です。

　実装のやり方は色々あるようですが、一番基本的な（だと思う）共分散行列の固有値と固有ベクトルを求める方法で行きます。

　やるべきこととしては、

データをセンタリングする（列ごとに平均を引く）
共分散行列を計算する
固有値と固有ベクトルを計算
データを固有ベクトルを使って写像する

　これらを実装すれば行けるはずです。というか、これで行くことにしました。

実装

　書いたソースコードを以下に示します。

# coding: UTF-8

import numpy as np

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA

import matplotlib.pyplot as plt

class MyPCA:
    def __init__(self, n_components=2):
        self.n_components = n_components

    def fit_transform(self, X):
        """横着してfit_transformしか実装してない
        """

        # 平均を0にする
        X = X - X.mean(axis=0)

        # 共分散行列を作る
        self.cov_ = np.cov(X, rowvar=False)
        
        # 固有値と固有ベクトルを求めて固有値の大きい順にソート
        l, v = np.linalg.eig(self.cov_)
        l_index = np.argsort(l)[::-1]
        self.l_ = l[l_index]
        self.v_ = v[:,l_index] # 列ベクトルなのに注意

        # components_（固有ベクトル行列を途中まで取り出す）を作る
        self.components_ = self.v_[:,:self.n_components].T

        # データとcomponents_をかける
        # 上と下で二回転置してるのアホ・・・
        T = (np.mat(X)*(np.mat(self.components_.T))).A

        # 出力
        return T

def main():
    iris = load_iris()

    pca = PCA(n_components=2)
    sklearn_X = pca.fit_transform(iris.data)

    my_pca = MyPCA()
    my_X = my_pca.fit_transform(iris.data)

    print(pca.explained_variance_)
    print(my_pca.l_)

    print(pca.components_)
    print(my_pca.components_)

    plt.figure()
    plt.scatter(sklearn_X[:,0], sklearn_X[:,1], c=iris.target/3)
    plt.savefig("sklearn_resut.png")

    plt.figure()
    plt.scatter(my_X[:,0], my_X[:,1]*-1, c=iris.target/3)
    plt.savefig("my_result.png")

if __name__ == "__main__":
    main()

　numpyを使ったので簡単に書けました。アルゴリズム部分はコメントで解説を入れたので、それを読めばどんな感じかは理解して頂けると思います。

結果

　mainのテキスト出力を見ると、次のようになっていました。

# 固有値
[4.22484077 0.24224357]
[4.22484077 0.24224357 0.07852391 0.02368303]

# components_
[[ 0.36158968 -0.08226889  0.85657211  0.35884393]
 [ 0.65653988  0.72971237 -0.1757674  -0.07470647]]
[[ 0.36158968 -0.08226889  0.85657211  0.35884393]
 [-0.65653988 -0.72971237  0.1757674   0.07470647]]

　固有値が余計に出ちゃってますが、これは別に構いません。また、componentsの2次元目が符号反転していますが、これも特に問題ないこと（のはず）なので無視します。

　自作の方は第二主成分を反転させてプロットしてみました。

sklearnのPCAでirisを可視化

自作PCAでirisを可視化

　同じ図を2つ載せるなって怒られそうですが・・・とにかく上手く行ったようです。

まとめ

　numpyで実装してみたら思ったより簡単だったので、これで当分は「わかった気」になれそうです。

　ただ、今回は特異値分解やらなかったので、それはまた宿題ということで・・・。

【python】カーネル主成分分析を試してみる

2018-03-28T00:19:12+09:00

　カーネル主成分分析（Kernel PCA）はカーネル法と主成分分析を組み合わせて用い、データを非線形次元圧縮する方法です（こんな説明で良いのか・・・）。

　カーネル法のことは勉強中・・・というか正直勉強しようとしてもよくわからないで跳ね返されるのをこれまで4回くらい繰り返してきたのですが、とりあえず使ってみました。

試してみた

　非線形データが手元にあると良いのですが、あいにくありません。輪っか状のデータなどを生成してやってみるのは簡単にできますが、面白くなさそうです。だいたいsklearnの公式サンプルにすらあります。
Kernel PCA — scikit-learn 0.21.2 documentation

　そこで、分類問題での適用を考えます。これならいつものようにPCA+CLFとKPCA+CLFで比較するだけなので、簡単そうです。更に、カーネルのgammaはグリッドサーチして最適値を探すだけ・・・。

　ただし、irisやdigitsで散々色々試してみましたが、ぶっちゃけ普通にやるとなかなかPCAを上回る性能が得られませんでした。最終的に、「digitsを3次元に次元削減し、LDAで分類する」という問題でどうにかそれなりに性能が上回ることがわかりましたが、実用的な意味はあまりありません。

　たぶん、sklearnのtoy datasetは低次元で線形分離できるタチの良いデータばっかりなのだと思います。それはそれで良いことですが、ちょっとタチの悪いデータも混ぜておいてもらえると嬉しいところです（かといって20newsgroupsのBoWだとタチが悪すぎるし・・・2000データ400次元くらいのちょうど良いデータはどこかにないものか）。

　コードを以下に示します。

# coding: UTF-8

import numpy as np

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA, KernelPCA as KPCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.pipeline import Pipeline
from sklearn.model_selection import GridSearchCV, StratifiedKFold as SKF
from sklearn.model_selection import cross_val_score

def main():
    dataset = load_digits()
    print(dataset.data.shape)

    pca = PCA(n_components=3)
    kpca = KPCA(kernel="rbf", n_components=3)
    lda = LDA()
    pl_pca = Pipeline([("pca", pca), ("lda", lda)])
    pl_kpca = Pipeline([("kpca", kpca), ("lda", lda)])

    parameters = {"kpca__gamma" : np.arange(0.00001, 0.003, 0.0001)}

    clf = GridSearchCV(pl_kpca, parameters, verbose=0, n_jobs=-1)

    print(cross_val_score(pl_pca, dataset.data, dataset.target, 
                          cv=SKF(shuffle=True, random_state=0),
                          scoring="f1_macro").mean())                      
    print(cross_val_score(clf, dataset.data, dataset.target,
                          cv=SKF(shuffle=True, random_state=0),
                          scoring="f1_macro").mean())    

if __name__ == "__main__":
    main()

　PCAでは0.68らい、KPCAでは0.71くらいのF1値が得られました。

　だから？　って言われると、返す言葉は思いつきませんが・・・。

まとめ

　やってみた記事ですが、何かの参考になればと思います。意外と上手く使うのは難しいと感じました。というか分類の次元削減としてはたぶんそんなに適当ではないです。

　どんな問題に応用されてるんだろうか。やっぱり可視化？

追記

　文字列の編集距離の可視化に使ってみました。

www.haya-programming.com

　文字列カーネルというのもあるらしいのですが、sklearnで対応していないし、未確認。編集距離を使う分には無難に使えます。

【python】sklearnのPCAで相関行列を使う

2018-03-27T02:41:44+09:00

　主成分分析には共分散行列を用いる方法、相関行列を使う方法がある。

　sklearnのPCAを見ると、これに対応するオプションは存在しない。

sklearn.decomposition.PCA — scikit-learn 0.20.1 documentation

　ずっと不思議に思っていたが、ググってたらこんなものを見つけた。

Enhance: PCA options for using Correlation or covariance matrix · Issue #2689 · scikit-learn/scikit-learn · GitHub

　要約：特徴量をスケーリングしてPCAすれば相関行列でやったのと同じことになるよ。PipelineでStandardScalerと組み合わせてね。おわり。

本当か確認する

　確認してみる。

>>> import numpy as np
>>> from sklearn.datasets import load_iris
>>> from sklearn.preprocessing import StandardScaler
>>> from sklearn.decomposition import PCA
>>> from sklearn.pipeline import Pipeline
>>> iris = load_iris()
>>> pca = PCA(n_components=2)
>>> pca.fit(iris.data)
PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
>>> pca.get_covariance()
array([[ 0.67919741, -0.03258618,  1.27066452,  0.5321852 ],
       [-0.03258618,  0.18113034, -0.31863564, -0.13363564],
       [ 1.27066452, -0.31863564,  3.11934547,  1.28541527],
       [ 0.5321852 , -0.13363564,  1.28541527,  0.58961806]])
>>> ss = StandardScaler()
>>> p = Pipeline([("scaler", ss), ("pca", pca)])
>>> p.fit(iris.data)
Pipeline(memory=None,
     steps=[('scaler', StandardScaler(copy=True, with_mean=True, with_std=True)), ('pca', PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False))])
>>> p.steps[1][1].get_covariance()
array([[ 0.9779242 , -0.10104477,  0.87069468,  0.86134879],
       [-0.10104477,  1.00395722, -0.41916911, -0.37286994],
       [ 0.87069468, -0.41916911,  1.04639367,  0.93676197],
       [ 0.86134879, -0.37286994,  0.93676197,  0.99857055]])
>>> np.corrcoef(iris.data, rowvar=False)
array([[ 1.        , -0.10936925,  0.87175416,  0.81795363],
       [-0.10936925,  1.        , -0.4205161 , -0.35654409],
       [ 0.87175416, -0.4205161 ,  1.        ,  0.9627571 ],
       [ 0.81795363, -0.35654409,  0.9627571 ,  1.        ]])

　違うじゃん。妥当そうなのはnumpyの結果だが（対角成分が1になってる）、とりあえずしょうがないのでスケーリングしたデータの共分散をnumpyで計算してみる。

>>> np.cov(ss.fit_transform(iris.data), rowvar=0, bias=1)
array([[ 1.00671141, -0.11010327,  0.87760486,  0.82344326],
       [-0.11010327,  1.00671141, -0.42333835, -0.358937  ],
       [ 0.87760486, -0.42333835,  1.00671141,  0.96921855],
       [ 0.82344326, -0.358937  ,  0.96921855,  1.00671141]])
>>> np.cov(ss.fit_transform(iris.data), rowvar=0, bias=1)
array([[ 1.        , -0.10936925,  0.87175416,  0.81795363],
       [-0.10936925,  1.        , -0.4205161 , -0.35654409],
       [ 0.87175416, -0.4205161 ,  1.        ,  0.9627571 ],
       [ 0.81795363, -0.35654409,  0.9627571 ,  1.        ]])

　標本分散はnp.corrcoefと等価だ。

　ここまでやったところでもう一回ドキュメントを読み、PCA.get_covariance()の結果が「Estimated covariance of data.」であり、厳密ではないことに気づいたので、問題は解決した。

　理論的にこうなる理由は、説明しようと思えばできるのだと思いますが、今回は大変なので触れません。

irisでやってみる

　irisの可視化にそれぞれを使ってみる。コードを以下に示す。

# coding: UTF-8

from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.pipeline import Pipeline

import matplotlib.pyplot as plt

def main():
    iris = load_iris()

    ss = StandardScaler()
    pca = PCA(n_components=2)
    p = Pipeline([("scaler", ss), ("pca", pca)])
    
    X = pca.fit_transform(iris.data)

    plt.figure()
    plt.scatter(X[:,0], X[:,1], c=iris.target/3)
    plt.savefig("iris_cov_pca.png")

    X = p.fit_transform(iris.data)

    plt.figure()
    plt.scatter(X[:,0], X[:,1], c=iris.target/3)
    plt.savefig("iris_corr_pca.png")

if __name__ == "__main__":
    main()

　結果は、

共分散行列で主成分分析したiris

相関行列で主成分分析したiris

　こうして見ると相関行列はあまりメリットがないように見えますが、実際には相関行列の方が良いタスクは色々あるようです。相関行列を使うことでbiplotが上手く行っているという例を出しているページを載せておきます。
PCA on correlation or covariance? - Cross Validated

まとめ

　とりあえずできることはわかったので良しとする。

　でも、「pipelineで出来るから要らねーよ」ってつもりらしいけど、ぶっちゃけオプション一つでできた方が親切だと思った（小並感）。

【python】sklearnのfetch_20newsgroupsで文書分類を試す(4)

2018-03-26T21:21:12+09:00

　前回は性能を追い求めると次元がでかくなりすぎて・・・というところで終わっていた。今回はもうちょっと頑張って次元を減らしてみる。

ストップワードの除去
PCA（主成分分析）とLDA（線形判別分析）
分類
ソースコード
結果とまとめ
次回
過去の回

ストップワードの除去

　とりあえずstop_wordsを指定していなかったので、指定してみる。

　stop_words="english"とすると、ストップワードを除去してくれる。

　結果だけ言うと、min_df=0.005のとき、

stop_words指定なし：3949次元
stop_words指定あり：3705次元

　だった。焼石に水。

PCA（主成分分析）とLDA（線形判別分析）

　PCAとLDAをかけ、次元削減をする。leakage怖いのでPipelineを使う（厳密なことを言い出すと、単語文書行列を作る段からPipelineに入れるべきなのだろうか？　きついのでパスさせて頂くが）。

　PCAは主にLDAの計算負荷削減と、変数の相関を除去することを意図してかける。1000次元まで落としてみたが、これでも累積寄与率は90%弱になる。まあ、正規化も何もしてないから、重要な情報を落としている可能性は否定できないのだが。

　LDAは次元削減に使う。有効性についてはこの前試してみたので、この記事を読んで欲しい。
【python】LDA（線形判別分析）で次元削減 - 静かなる名辞
　20newsgroupsは20クラスのデータなので、19次元に落とすことになる。相当早くなるだろうが、どこまで性能を維持できるかはデータの線形性にかかっている。

分類

　ランダムフォレストを使った。n_estimators=1000とし、他のパラメタはデフォルト。

ソースコード

　実験に使ったソースコードを以下に示す。

# coding: UTF-8

import numpy as np

from sklearn.datasets import fetch_20newsgroups
from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.feature_extraction.text import CountVectorizer as CV
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import precision_recall_fscore_support as prf
from sklearn.pipeline import Pipeline

def main():
    news20 = fetch_20newsgroups()    
    
    cv = CV(min_df=0.005, max_df=0.5, stop_words="english")
    matrix = cv.fit_transform(news20.data).toarray()

    pca = PCA(n_components=1000, svd_solver="randomized")
    lda = LDA()
    rfc = RFC(n_estimators=1000, n_jobs=-1)

    clf = Pipeline([("pca", pca), ("lda", lda), ("rfc", rfc)])

    trues = []
    preds = []
    for train_index, test_index in StratifiedKFold().split(matrix, news20.target):
        clf.fit(matrix[train_index], news20.target[train_index])
        trues.append(news20.target[test_index])
        preds.append(clf.predict(matrix[test_index]))
    scores = prf(np.hstack(trues), np.hstack(preds), average="macro")[:3]
    print("p:{0:.6f} r:{1:.6f} f1:{2:.6f}".format(scores[0],
                                                  scores[1],
                                                  scores[2]))

if __name__ == "__main__":
    main()

結果とまとめ

p:0.764012 r:0.760731 f1:0.761510

　前回の0.8を超えるスコアには届かなかったが、とりあえずそれなりに軽くはなった。もうちょっと真面目に追い込めばという話はあるが、追求しない。次回はもうちょっと違うことをやってみたい。

次回

　このシリーズずっと放置していましたが、気が向いたので書きました。
www.haya-programming.com

過去の回

【python】sklearnのfetch_20newsgroupsで文書分類を試す(1) - 静かなる名辞
 【python】sklearnのfetch_20newsgroupsで文書分類を試す(2) - 静かなる名辞
 【python】sklearnのfetch_20newsgroupsで文書分類を試す(3) - 静かなる名辞

【python】LDA（線形判別分析）で次元削減

2018-03-20T16:43:52+09:00

　一般によく使われる次元削減手法としてはPCA（主成分分析）がありますが、他にLDA（Linear Discriminant Analysis：線形判別分析）を使う方法もあります。

　これは本来は分類に使われる判別分析という古典的なアルゴリズムで、データが一番分離しやすくなる軸を求めていくものです。つまり教師ラベルを使います。教師ラベルを使うので、PCAのような教師なしの手法と比べて有利な可能性があります。

　線形判別分析の詳しい原理の説明などが欲しい方は、ググって出てくるwikipediaやqiitaなどを参考にしてください（投げやり）。この記事では、分類問題でこれを使ったとき、どのようなご利益があるのかを検証します。

実験

　sklearnのdigitsデータセットを使い、次元削減→分類というタスクを行って交差検証でスコアを出します。

　分類器は最初はSVMでやろうかと思ったけど、パラメタチューニングで幾らでも恣意的な結果になることに気づいたのでガウシアン・ナイーブベイズでやることにしました。

　実験に使ったコードは以下に示します。

# coding: UTF-8

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.naive_bayes import GaussianNB as GNB
from sklearn.pipeline import Pipeline
from sklearn.model_selection import StratifiedKFold as SKF
from sklearn.metrics import precision_recall_fscore_support  as prf

def main():
    digits = load_digits()

    gnb = GNB()

    df = pd.DataFrame([], columns=[
        "n_components",
        "pca-gnn precision", "pca-gnn recall", "pca-gnn f1",
        "lda-gnn precision", "lda-gnn recall", "lda-gnn f1"])
    for n_components in [5, 10, 15, 20, 25, 30, 40]:
        pca = PCA(n_components=n_components)
        lda = LDA(n_components=n_components)

        steps1 = list(zip(["pca", "gnb"], [pca, gnb]))
        steps2 = list(zip(["lda", "gnb"], [lda, gnb]))

        p1 = Pipeline(steps1)
        p2 = Pipeline(steps2)

        score_lst = []
        for decomp_name, clf in zip(["pca", "lda"], [p1, p2]):
            trues = []
            preds = []
            for train_index, test_index in SKF(
                    shuffle=True, random_state=0).split(
                    digits.data, digits.target):
                clf.fit(digits.data[train_index], 
                        digits.target[train_index])
                trues.append(digits.target[test_index])
                preds.append(clf.predict(digits.data[test_index]))
            scores = prf(np.hstack(trues), np.hstack(preds), average="macro")
            score_lst.extend(scores[:-1])
        df = df.append(pd.Series([n_components, *score_lst],
                                 index=df.columns),
                       ignore_index=True)
    print(df)
    plt.figure()
    df.plot(x="n_components", y=["pca-gnn f1", "lda-gnn f1"])
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

結果

　次のようになりました。

　テキスト出力

   n_components  pca-gnn precision  pca-gnn recall  pca-gnn f1  \
0           5.0           0.847918        0.841684    0.841109   
1          10.0           0.915834        0.911346    0.912563   
2          15.0           0.926992        0.923032    0.924061   
3          20.0           0.934522        0.930192    0.931194   
4          25.0           0.941886        0.938611    0.939205   
5          30.0           0.946139        0.944251    0.944669   
6          40.0           0.945330        0.943644    0.943960   

   lda-gnn precision  lda-gnn recall  lda-gnn f1  
0           0.917464        0.917144    0.917031  
1           0.953751        0.952588    0.952950  
2           0.953751        0.952588    0.952950  
3           0.953751        0.952588    0.952950  
4           0.953751        0.952588    0.952950  
5           0.953751        0.952588    0.952950  
6           0.953751        0.952588    0.952950

結果（n_components対F1値）

　LDAを使った方が低い次元で、より高い分類性能が得られているようです。

まとめ

　LDAは良い。

おまけ

　ソースコードをちゃんと読んだ方は、最初に書かれた以下の記述に気づいたかと思います。

import warnings
warnings.filterwarnings('ignore')

　これを付けないとLDAはけっこうな警告（主に以下の2つ）を吐いてくれます。

UserWarning: Variables are collinear
UserWarning: The priors do not sum to 1. Renormalizing

　上の警告はPCAで説明変数の多重共線性を除去してやると消えます（本末転倒っぽいけど）。下の警告は、正直調べてもよくわかりませんでした。

　とりあえず、警告が出てもちゃんと動いてるみたいなので別に良いか・・・。

追記

　LDAのn_componentsには上限があり、クラス数-1以上のn_componentsは指定しても無意味です。

　実際にやってみても、クラス数-1以上にはなりません。

>>> from sklearn.datasets import load_digits
>>> from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
>>> lda = LDA(n_components=15)
>>> lda.fit(digits.data, digits.target)
>>> lda.explained_variance_ratio_
array([0.28912041, 0.18262788, 0.16962345, 0.1167055 , 0.08301253,
       0.06565685, 0.04310127, 0.0293257 , 0.0208264 ])

　決定境界をクラス数-1個引くので（SVMで言うところのone-versus-the-rest）、n_componentsも必然的にそれだけ必要になります（逆にそれ以上は必要になりません）。

　上のグラフはそのつもりで眺めてください。また、LDAはけっきょくのところ線形変換なので、クラス数-1次元の線形空間にうまく張り直せないような入力に対しては無力なことも覚えておく必要があるでしょう（PCAも非線形構造はダメだが・・・カーネルでも持ってくる必要がある）。

【python】sklearnのPCAでsvd_solverによる速度差を比較

2018-03-19T17:23:15+09:00

　sklearnのPCA（主成分分析）がやたら遅くて腹が立ちました。計算コストを下げるために次元削減してるのに、次元削減で計算コスト食ったら意味がありません。

　とにかくこのPCAを高速化したかったので、svd_solverを変えてどうなるか試しました。なお、腹が立つくらい遅かった理由は最終的にちゃんとわかったので、この記事の最後に載せます。

svd_solverとは
実験
結果
まとめ
おまけ：腹が立った理由

svd_solverとは

　PCAは内部で特異値分解（SVD）を使っています。この特異値分解がコンピュータにやらせるにはそれなりに計算コストの高い処理で、とりあえずアルゴリズムが何種類かあるようです。

　sklearnのPCAで使える（指定できる）アルゴリズムは次の4つです。

auto

　デフォルト値。500*500以下の入力データならfullを、それ以上ならrandomizedを使うそうです*1

full

　standard LAPACK solverを使うそうです。とりあえずぜんぶ丸ごと特異値分解してから、n_componentsで指定した次元数だけ取ってくるそうな

arpack

　Truncate SVDという手法を使う。一次元ずつ寄与率の大きい主成分から計算していくらしい。n_componentsが小さければ速いことが期待されるんだと思う

randomized

　randomized SVDという手法で計算する。乱数使って速くした。乱数なので厳密解ではない

　なお、以上の情報はすべて公式ドキュメントから得ました。
sklearn.decomposition.PCA — scikit-learn 0.20.1 documentation

　とりあえずautoはどうでも良いので、残りの3つを比較することにします。

実験

　PCAをかけたくなるような高次元データといえばBag of Words、ということでこのブログですでに何回も取り上げたことのある、sklearnのfetch_20newsgroupsとCountVectorizerの組み合わせを使います。前者はテキストのデータセット、後者はBoWを生成するクラスです。

　次のような実験用コードを書きました。

# coding: UTF-8

import time
from itertools import product

from sklearn.datasets import fetch_20newsgroups
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.decomposition import PCA

def main():
    news20 = fetch_20newsgroups()

    for min_df in [0.02, 0.01, 0.008, 0.005]:
        cv = CountVectorizer(min_df=min_df, max_df=0.5,
                             stop_words="english")
        X = cv.fit_transform(news20.data).toarray()

        print("min_df:{0} X.shape:{1}".format(min_df, X.shape))
        for n_components, svd_solver in product(
                [100, 500],
                ["full", "arpack", "randomized"]):
            pca = PCA(n_components=n_components, svd_solver=svd_solver)
            t1 = time.time()
            pca.fit_transform(X)
            t2 = time.time()
            print("n_components:{0}  solver:{1:>10}  "\
                  "time:{2:>6.2f}  CP:{3:.4f}".format(
                      n_components, svd_solver, t2-t1, 
                      pca.explained_variance_ratio_.sum()))
        print("")

if __name__ == "__main__":
    main()

　BoWの次元数をmin_dfで変えていき、n_componentsを100と500、svd_solverを上記3つで変化させてPCAをかけたときの速度と累積寄与率（CP：Cumulative Proportion）をそれぞれ測ります。

結果

　次のようになりました。

min_df:0.02 X.shape:(11314, 866)
n_components:100  solver:      full  time:  3.60  CP:0.7455
n_components:100  solver:    arpack  time:  3.90  CP:0.7455
n_components:100  solver:randomized  time:  1.72  CP:0.7443
n_components:500  solver:      full  time:  3.89  CP:0.9528
n_components:500  solver:    arpack  time: 19.42  CP:0.9528
n_components:500  solver:randomized  time:  8.91  CP:0.9516

min_df:0.01 X.shape:(11314, 1916)
n_components:100  solver:      full  time: 22.38  CP:0.8029
n_components:100  solver:    arpack  time:  8.41  CP:0.8029
n_components:100  solver:randomized  time:  4.86  CP:0.8028
n_components:500  solver:      full  time: 22.06  CP:0.9304
n_components:500  solver:    arpack  time: 53.73  CP:0.9304
n_components:500  solver:randomized  time: 13.47  CP:0.9293

min_df:0.008 X.shape:(11314, 2391)
n_components:100  solver:      full  time: 34.24  CP:0.7899
n_components:100  solver:    arpack  time: 10.42  CP:0.7899
n_components:100  solver:randomized  time:  5.75  CP:0.7897
n_components:500  solver:      full  time: 34.88  CP:0.9193
n_components:500  solver:    arpack  time: 63.37  CP:0.9193
n_components:500  solver:randomized  time: 15.18  CP:0.9182

min_df:0.005 X.shape:(11314, 3705)
n_components:100  solver:      full  time:100.52  CP:0.7701
n_components:100  solver:    arpack  time: 16.46  CP:0.7701
n_components:100  solver:randomized  time:  8.70  CP:0.7699
n_components:500  solver:      full  time:100.73  CP:0.9000
n_components:500  solver:    arpack  time: 94.33  CP:0.9000
n_components:500  solver:randomized  time: 20.04  CP:0.8988

　要約すると、

fullは基本的に遅い。入力の次元数が増えるとびっくりするくらい遅くなる
arpackは100次元に落とすときは威力を発揮している。500次元に落とすケースではかえって遅くなる。ヘタするとfullより遅い
randomizedは速い。ただし厳密解ではないことがCPからわかる（full、arpackとは微妙に違う数字になっている）

　こういう状況です。わかりやすいですね。

　それぞれの使い分けは、

入力次元数の小さい入力ではfullで良い。というかヘタにそれ以外を指定するとかえって遅いケースもある
入力次元数が大きく、入力次元数>>出力次元数で厳密解がほしければならarpackの使用を検討する
厳密解じゃなくても良いのでとにかく速いのを！　ってときはrandomized

　ってことになるかと思う・・・。

まとめ

　けっこう変わる。頑張って使い分けよう。

おまけ：腹が立った理由

　sklearnのPCAではn_componentsに小数を指定できます。そうすると累積寄与率がその数字になるように勝手に次元数を決めてくれるので、こりゃ便利だわいと思って私はよく使っていました。

　しかし、実はarpack、randomizedではこの小数での指定は使えません。そのことはドキュメントにもちゃんと書いてあります。無理矢理に指定すると次のようなエラーを吐かれます。

ValueError: n_components=0.95 must be between 1 and n_features=866 with svd_solver='arpack'

　ということは何が起こるか？　勝手にfullにされます。遅い訳です。なんてこった。

　わかってしまえば下らない話で、要するに私が使いこなせていなかっただけなのですが、このことは「ちゃんとドキュメントをよく読んで使おうね」という教訓を私に残したのでした。

*1:300*800だったりしたらどうなるんだろう？　それとも共分散行列のサイズなのだろうか？