Pythonで主成分分析をするコード


2012/02/29 修正しました

以前C++からOpenCVで主成分分析するコードを書きましたが、今回はOpenCVも使わず、Pythonで実装してみました。

import numpy
import numpy.linalg
class PCA(object):
    """主成分分析
        eigenvectors:行ベクトルとして固有ベクトルが格納されている
        eigenvalues:固有値の列.大きい順にソートされている
        mean:平均ベクトル(行ベクトル)
    """
    def __init__(self, dim):
        self.dim=dim

    def calc(self,src):
        """主成分分析の計算を行う
        src:各データが行ベクトルとして格納されている
        """
        shape=src.shape

        # 中心化して転置する.以降,各データは列ベクトル
        mu = (sum(numpy.asarray(src,dtype=numpy.float)) / float(len(src)))
        src_m = (src - mu).T

        if shape[0] < shape[1]: # データ数 < 次元数のとき
            # PCA後の最大次元数の決定
            if shape[0] < self.dim: self.dim=shape[0]

            print>>sys.stderr,"covariance matrix ...",
            n=numpy.dot(src_m.T, src_m)/float(shape[0])
            print>>sys.stderr,"done"

            print>>sys.stderr,"eigen value decomposition ...",
            l,v=numpy.linalg.eig(n)
            idx= l.argsort()
            l = l[idx][::-1]
            v = v[:,idx][:,::-1]
            print>>sys.stderr,"done"

            vm=numpy.dot(src_m, v)
            for i in range(len(l)):
                if l[i]<=0:
                    v=v[:,:i]
                    l=l[:i]
                    if self.dim < i: self.dim=i
                    break
                vm[:,i]=vm[:,i]/numpy.sqrt(shape[0]*l[i])

        else: # データ数 >= 次元数のとき

            if shape[1] < self.dim: self.dim=shape[1]
            cov=numpy.dot(src_m, src_m.T)/float(shape[0])
            l,vm = numpy.linalg.eig(cov)
            idx= l.argsort()
            l = l[idx][::-1]
            vm = vm[:,idx][:,::-1]

        self.eigenvectors=vm[:,:self.dim]
        self.eigenvalues=l[:self.dim]
        self.mean=mu

    def project(self,vec):
        return numpy.dot(self.eigenvectors.T,(vec-self.mean).T).T

    def backproject(self,vec):
        return (numpy.dot(self.eigenvectors,vec.T)).T+self.mean

画像処理でPCAをする場合、データ数に対して次元数がとても大きい時が多いのですが、その場合、普通に共分散行列を作るとメモリが溢れてしまいます。 例:ベクトルが100000次元、データ数1000の場合、共分散行列の計算を普通にやると100000\times 100000の行列になってしまいます。 そこで、共分散行列 XX^T の代わりに X^TXの行列を計算するやり方があります。 例:このやり方だと、1000\times 1000の行列で済みます。 上のコードでは、入力したデータの次元数、データ数に応じてこの2つの方法を切り替えています。 使い方の例は次のようになります。

 def main():
    import numpy     # データの読み込み(各データは行ベクトルとして入っている)
    data = numpy.loadtxt("data.csv",delimiter=",")
    pca=PCA(2)     # PCAにより固有ベクトルなどを求める
    pca.calc(data)     # pca.eigenvectors, pca.eigenvalues, pca.meanで固有ベクトル、固有値、平均が得られる
    # 固有ベクトルは列ベクトルとして得られる
    # 固有空間への投影(行ベクトルとして結果が得られる)
    print>>sys.stderr,"projecting ...",
    p1=pca.project(data)
    print>>sys.stderr,"done"

    # 固有空間からの逆投影
    print>>sys.stderr,"back projecting ...",
    v1=pca.backproject(p1)
    print>>sys.stderr,"done"

if __name__=="__main__":
    main()

これもシリーズ物で、次回以降、確率的主成分分析やカーネル主成分分析について紹介して行きます。

Ubuntu上のphpでcURLを使う


curlは指定したURLからデータを取ってくるコマンドですが,それをphp5から利用する場合,php5-curlが必要になります.

コマンドライン上で

sudo apt-get install php5-curl

を実行することでphp5のcurlモジュールをインストールすることができ,利用できます.

何に使うかというと…wordpressにSimple Tweetというプラグインを入れたのですが,

Twitter OAuth は PHP5 以降でのみ利用できます。

というエラーメッセージが出て使えなかったからです.
このエラーメッセージはそのサーバのPHPでcurlが使えないときにも出るようです.

[shogun 1] shogunを使ってみよう


最近Shogun - A Large Scale Machine Learning Toolbox という機械学習のツールボックスを使い始めたので,shogunについて紹介します.

まず紹介ページの英語を適当に日本語に訳してみます.

SHOGUNは大規模なカーネル法と特にサポートベクタマシン(Support Vector Machine: SVM)[1]に重点を置いている機械学習ツールボックスである.いろいろなSVM実装(OCAS [21], Liblinear [20], LibSVM [2], SVMLight, [3] SVMLin [4] and GPDT [5])に対するインターフェースになる汎用的なSVMオブジェクトを提供する.SVMは様々なカーネルと組み合わせることができる.このツールボクスは線形カーネル,多項式カーネル,ガウシアンカーネル,シグモイドカーネルだけでなく,Locality Improved [6], Fischer [7], TOP [8], Spectrum [9], Weighted Degree Kernel (with shifts) [10] [11] [12]なども提供する.また効率的なLINADD [12]最適化法も実装している.線形SVMに対して,COFFIN framework [22][23]は疎行列,密行列他様々なデータタイプの混合に対してもオンデマンドで処理ができる.さらにSHOGUNは様々なpre-computedなカーネルで動くようになっている.重要なものとしては,サブカーネルの重み付き線形和による結合カーネルがある.この各サブカーネルは必ずしも同じドメインである必要はない.最適なサブカーネルの重みはMultiple Kernel Learning (MKL) [13][14][18][19]で学習できる.現在SVMは,1クラス,2クラス,多クラス分類,回帰問題を扱える.SHOGUNは線形判別分析(Linear Discriminant Analysis: LDA),線形計画マシン(Linear Programming Machine: LPM),パーセプトロン((Kernel) Perceptrons),隠れマルコフモデル(Hidden Markov Model)も実装されている.入力データは様々な型の疎・密行列を利用することができる.平均を引くなどの一連の前処理を各入力データに適用することもできる.
SHOGUNはC++で実装されており,Matlab, R, Octave, Python, Java, Ruby向けのインタフェースが用意されており,機械学習のオープンソースソフトウェアとして公開されている.

長々と書きましたが,重要なのは以下のポイントになります.

  • 様々なSVM実装が利用できる
  • 様々なカーネルを使える
  • SVM以外にもLDAとかHMMとかいろいろ使える
  • C++/Matlab/R/Octave/Python/Ruby/Javaなどで使える
  • オープンソースソフトウェアである

関西CV・PRML勉強会で紹介した時のスライドはこちらにあります.

僕はこのツールボックスを,Pythonから利用しはじめました.スライド中ではapt-getやMacPortsでインストールできると言っていますが,サンプルコードなどを入手したい場合は,gitを使ってソースコードをダウンロードし,コンパイルしてインストールするのがお勧めです.

gitがインストールされているならば

git clone git://github.com/shogun-toolbox/shogun.git
git clone git://github.com/shogun-toolbox/shogun-data.git

のコマンドで,shogunディレクトリとshogun-dataディレクトリができます.shogunディレクトリにはshogunのソースコードとサンプルコードが,shogun-dataディレクトリにはサンプルのデータが入っています.

ソースコードをダウンロードしたら,

cd shogun/src
./configure
make

し,管理者権限で

make install

することでインストールすることができます.
./configureの前に,swigがインストールされていないとpythonなどからの利用ができないので,先にインストールしておく必要があります.これはapt-getやMacPortsなどでインストールすることができます.

サンプルコードを動かすためには,サンプルデータのディレクトリの設定をする必要があります.shogunディレクトリ中で

rmdir data
ln -s ../shogun-data data

としてリンクを張っておく,もしくはshogun-dataの中身をshogun/dataへ移動すると完了です.

次のエントリではPython版のサンプルコードを動かすためのノウハウについて書きます.