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

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