べき乗法を使った主成分分析の Ruby 実装

The implementation of ruby for principal component analysis(PCA) that uses power method.

◆Ruby には NArray なんてもんがある

　仕事で Ruby を使うようになってから、ちょっとした試作は Ruby を使うようになりました。まぁなんといっても総じてコーディング量が少なくなるから、後で見返す時に見通しが良いからです。でも、Ruby は実行速度が遅いのでパフォーマンスを求められるようなところにはあんまし向きません。

　と思っていたら、NArray なる拡張ライブラリがあって、派生クラスの NMatrix が行列演算をすることを知りました。拡張ライブラリは C 言語で書かれていますので、計算をこれに任せてしまえば、そこそこのスピードで動作するはずです。また、めんどくさい行列計算用のコードをあまり書かなくて良いところも魅力的ですね。

◆ちょっと気になるべき乗法

　で、こちらのページを参考にさせて頂き主成分分析を実装してみました。主成分分析は、主成分を導くときに固有値と固有ベクトルを求めなければなりません。よく知られているヤコビ法などを使う場合が多いと思いますが、参考ページにはべき乗法（パワー法）による固有値と固有ベクトルの求め方の説明があります。この方法は Google が page ランクを計算するときに採用している手法でもあり少し気になっていました。理屈はともかくとして、計算は難しくないのでこれを実装することにしました。

◆書いたソースコード

#!/usr/bin/env ruby
require 'narray'

# correlation matrix
def corr(m)
  n=m.shape[0]
  ret = NMatrix.float(n,n)

  v=[]
  (0...n).each{|i|
    v << m[i,true].sbt!(m[i,true].mean)
  }
  
  (0...n).each{|i|
    (0...n).each{|j|
      ab = v[i][].mul!(v[j]).sum
      aa = v[i][].mul!(v[i]).sum
      bb = v[j][].mul!(v[j]).sum
      ret[i,j]= ab / Math::sqrt(aa * bb)
    }
  }
  ret
end

# power method
def eig(a)
  # create a initial value
  x = NMatrix.float(a.shape[0],1).fill(1).transpose
  n = 1
  max = 1.0
  pmax = 10.0
  while((max - pmax).abs > 1e-5)
    a_x = a * x
    pmax = max
    max = a_x.max
    x = a_x / max
    n += 1
  end

  # [Eigenvalue,Eigenvector]
  [a_x.max,x / Math::sqrt(x[].mul!(x).sum)]
end

# normalization
def nrm(a)
  ret = a.dup
  a.shape[0].times{|i|
    mean = a[i,true].mean
    stddev = a[i,true].stddev
    ret[i,0] = (a[i,true].sbt!(mean) / stddev)[0,true]
  }
  ret
end

m = NMatrix[ [8,9,4],
             [2,5,7],
             [8,5,6],
             [3,5,4],
             [7,4,9],
             [4,3,4],
             [3,6,8],
             [6,8,2],
             [5,4,5],
             [6,7,6.0] ]
a = corr(m)

lambda1,x1 = eig(a)

p lambda1
p x1

# next a
a = a - lambda1 * x1 * x1.transpose
lambda2,x2 = eig(a)

p lambda2
p x2

# next a
a = a - lambda2 * x2 * x2.transpose
lambda3,x3 = eig(a)

p lambda3
p x3

x = NMatrix.float(3,3)
x[0,0] = x1[0,true]
x[1,0] = x2[0,true]
x[2,0] = x3[0,true]
p x

# principal component score
pcs = nrm(m * x)
p pcs

puts "done"

◆計算手順

　参考ページに従って変量 m の行列を定義します。例題は10紙の新聞のニュース、ビジネス、スポーツ記事について10点評価で調査する想定です。従って、10行3列の行列を作ります。

　行列 m を関数 corr に渡し相関行列 a を作ります。corr では 3記事分の相関係数を計算しますので 3行3列の対象行列を返します。

　相関行列 a の固有ベクトルを求めることで、主成分得点(変量を主成分の軸に射影した時の場所)が計算できるようになります。固有ベクトルは Google も使っているべき乗法で計算します。相関行列 a を引数に関数 eig を呼び出します。

　関数 eig は、初期値 x = (1,1,1)t を作り、ax = a・x を計算します。 x = ax / ax.max を次回の x として再び ax = a・x を計算します。このとき、前回の ax.max と今回の ax.max を比較し、その差が 1e-5 程度になったら計算を終了します。ベクトル x を正規化(長さを 1 にする計算)して戻り値の固有ベクトルとします。このとき、ax.max が固有値になります。

　2番目からは、a = a - λ・x・xt を関数 eig の引数として計算を繰り返します。

　各固有ベクトルが求まったら、変量の行列 m と固有ベクトルをまとめて作った行列 x を掛けて主成分得点を計算しプログラムは終了します。

　最後に、このソースコードは Ruby 1.9.1 で動作しています。
Latest update at 2009/4/13