scikit-learnを使ってPCAによる次元削減を7行で書いてみた

高次元のデータを解析する場合、データの分布などが見てもわからないので、PCAを使って次元削減をしてから可視化する場合があります。
PCAについてはこれまでにこのブログでも何度も触れているのでそちらを参照してください。

今回はフランスのINRIAがFundingしているscikit-learnというライブラリを使って7行で(空行含まず)書いてみました。言語はPythonです。

import numpy
import sklearn.decomposition

dim=2
data=numpy.loadtxt("data.csv", delimiter=",")
pca=sklearn.decomposition.PCA(dim)

result=pca.fit_transform(data);
numpy.savetxt("result.csv", result, delimiter=",")

たったこれだけで、data.csvというファイル(行ごとに多次元データがカンマ区切りで入っている)を読み込み、PCAで次元削減(2次元へ圧縮)したあとresult.csvというファイルに書き出しています。

白色化する場合は pca.whiten=True という行を、fit_transform の前に入れればOKだと思います。

PCAというより、numpy.loadtxt / numpy.savetxt 素晴らしい!という内容になりました。

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()

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

OpenCVで主成分分析をしてみた

主成分分析とは

主成分分析:PCA (Principal Component Analysis)は,多次元データの解析法の一種で,多次元空間中のデータ分布のうち,最も分散の大きくなる方向から順に基底を取っていく手法です.これをすることにより,データにおける主な変化の傾向を知ることができます.具体的にはいろんなところ(朱鷺の杜Wiki , タコでもわかる主成分分析 , etc.)で紹介されているので,ここでは割愛します.

主成分分析をやるには

大量の高次元ベクトルに対して上記のように分散を計算していくのは手計算では無理です.普通はmatlabやRを使いますが,画像処理屋はOpenCVを使って解析をします.(僕は主に,自前で実装したIncremental PCAを使っているので,以下のコードを実際には使っていません.あと,ちゃんとしたデータに対して動作確認もしてないので,以下のコードを使うのは自己責任でお願いします.)

このコードは,C++とOpenCVを使って,200次元のデータ1000個に対して主成分分析を行う例です.この例では,主成分分析後,上位3次元の基底に着目し,200次元のデータを3次元の空間へ投影した時の座標を出力します.また,各基底がどのくらい強くその分布を表しているかを意味する,累積寄与率も計算して出力しています.入力データは,スペースで区切られた200次元,1000個のベクトルです.

#include <cv.h>
#include <highgui.h>
#include <iostream>
#include <fstream>
#include <iomanip>
int main(){
	const int DIM=200; // 200次元
	const int SAMPLES=1000; // 1000サンプル
	const int RDIM=3; // 圧縮後3次元
	cv::Mat src(SAMPLES,DIM,CV_32FC1);
	cv::Mat result; // DIMとSAMPLESのうち小さい方次元xサンプル数次元
	double val;

	std::ifstream ifs("data.txt");

	for(int j=0;j<SAMPLES;j++){
		for(int i=0;i<DIM;i++){
			ifs >> val;
			((float*)src.data)[j*src.cols+i]=val;
		}
	}

	cv::PCA pca(src,cv::Mat(),CV_PCA_DATA_AS_ROW,RDIM);

	result=pca.project(src);

	for(int j=0;j<SAMPLES;j++){
		for(int i=0;i<RDIM;i++){ // 上位RDIM次元だけを表示
			std::cout << std::dec << ((float*)result.data)[j*result.cols+i] << ",";
		}
		std::cout << std::endl;
	}

	cv::Mat evalues=pca.eigenvalues;
	float sum=0.0f;
	for(int i=0;i<pca.eigenvalues.rows;i++){
		sum+=((float*)evalues.data)[i];
	}

	float contribution=0.0f;
	for(int i=0;i<RDIM;i++){// 上位RDIM次元の寄与率を計算
		contribution+=((float*)evalues.data)[i] / sum;
		std::cout << i+1 << "次元の累積寄与率:" << contribution << std::endl;
	}

	return 0;
}

このように,C++/OpenCVで簡単に主成分分析をすることができます.