2012-02-12

BPRで遊んだよ！

機械学習

前回ちょっとだけ述べた、last.fmのデータがここに転がってたので、
こいつで、BPRアルゴリズムを試してみることにしました。
使用したのは、usersha1-artmbid-artname-plays.tsvファイルのデータのみです。
usersha1-profile.tsvは今回は使用してません。
usersha1-artmbid-artname-plays.tsvのフォーマットですが、
ユーザID、アーティストID、アーティスト名、再生数がタブ区切りで記述されています。
ただし、アーティストIDが付けれられていないアーティストもいるので、
その場合はアーティストIDの部分は空文字列です。

ただし、そのまま使うのは、色々と厳しいので、
以下のルールでデータを抽出することにしました。

・ユーザ $u$ が、アーティスト $i$ を100回以上再生している場合に、
$u$ は $i$ を好んでいるとみなす（ $S_{ui} = 1$ ）。
・100再生以上している場合、再生数に差があっても嗜好度に差をつけない。
・100再生未満のデータは学習に利用しない。
・アーティストIDがないアーティストは抽出しない。

従って、100回以上再生したアーティストがいないユーザや、
誰からも100回以上再生されていないアーティスト、
IDがないアーティストのデータは、今回の実験では除外され、
推薦されることはありません。

このルールによって抽出されたデータ以下の通り。
ユーザー数：7748
アーティスト数：28185

ただ、このデータをただアルゴリズムに突っ込んで、
精度出しましたって言うのは微妙な気がするので、
今回、次の3人の疑似ユーザを学習に紛れ込ませ、
どんなアーティストが推薦されるかをみることにしました。

・某アニメの元ネタくん（ユーザID：K-ON）
好きなアーティスト
p-model、電気グルーヴ、triceratops
soft ballet、筋肉少女帯、the pillows
計6組
（ニコニコ大百科の記事を参考にしました）

・暗いよ70年代後半〜80年代くん（ユーザID：gotham）
好きなアーティスト
joy division、virgin prunes、the birthday party、swans
sisters of mercy、theatre of hate、bauhaus、cristian death
killing joke、Siouxsie & The Banshees、Alien Sex Fiend
Mission、chameleons
計13組

・最近のジョギングのお供だよくん（ユーザID：joging）
好きなアーティスト
sonic youth、melt-banana、!!!、gang of four
echo & the bunnymen、the rapture、the stills、kula shaker
mars volta、doors、espers、nick drake、柴田淳
 光田康典、植松伸夫、L'Arc en ciel、凛として時雨
 king crimson、yes、screaming headless torsos
計20組

これで準備完了。
後は、このユーザ数7751（7748+3）、アーティスト数28185のデータを
BPRアルゴリズムで学習して、推薦結果を得るだけです。

てことで、ソースコードです。
今回もpythonで書きました。
以下は、行列 $S$ からBPRのMatrix Factorization（BPR MF）を行い、
WとHを取得する関数です。
引数ですが、Sは論文に示されている $\|U\| \times \|I\|$ の行列 $S$ です。
未観測の要素は0で補間します。
つまり、ユーザ $u$ がアーティスト $i$ を好んでいるなら $S_{ui}=1$ 、
そうじゃない（好みじゃない、または、知らない）なら $S_{ui}=0$ となります。
kは $W$ と $H$ の列数（隠れ変数の数ともいえます）、
alphaは学習係数 $\alpha$ 、Lamは正則化パラメータ $\lambda$ です。

def MatrixFactorizationBPR(S , k , alpha , Lam):
    #W , Hを初期化
    UserCount = S.shape[0]
    ItemCount = S.shape[1]
    W = numpy.random.randn(UserCount , k)
    H = numpy.random.randn(ItemCount , k)
    W_next = W.copy()
    H_next = H.copy()

    BPRopt_pre=0.0
    count = 0

    while(1):
        #ループ回数をカウント
        count = count + 1

        #u,i,jをランダムサンプリング
        (u,i,j) = Get_uijFromArray(S);

        #xuijを計算
        x_uij = numpy.dot(W[u , :] , H[i , :]) - numpy.dot(W[u , :] , H[j , :]);

        #次のWとHを計算
        G =  numpy.exp(-x_uij) / (1 + numpy.exp(-x_uij));
        W_next[u , :] = W[u , :] + alpha * (G * (H[i , :] - H[j , :]) - Lam * W[u , :]);
        H_next[i , :] = H[i , :] + alpha * (G * (W[u , :]) - Lam * H[i , :]);
        H_next[j , :] = H[j , :] + alpha * (G * (-W[u , :]) - Lam * H[j , :]);

        #WとHを更新
        W = W_next.copy();
        H = H_next.copy();

        #計算に時間がかかるので、50万回に1回のみ収束判定
        if count%500000 == 0:
            BPRopt = GetBprOpt(S , W , H , Lam);

            #BPRoptの値で収束判定
            if(abs((BPRopt - BPRopt_pre))<0.001):
                break;
            BPRopt_pre = BPRopt;

    return W , H

MatrixFactorizationBPR内にある、関数Get_uijFromArrayは、
(u , i , j)の三つ組をランダムサンプリングする関数、
また、GetBprOptは、 $BPR-Opt$ を計算する関数です。
これらの関数はそれぞれ、以下の通りです。

Get_uijFromArray

def Get_uijFromArray(S):
    #ユーザ数とアイテム数をカウント
    UserCount = S.shape[0];
    ItemCount = S.shape[1];

    #ユーザをランダムに選択
    u = numpy.random.randint(UserCount);

    #選択したユーザの嗜好情報を取得
    S_u = S[u,:]
    #集合I+とI-を生成
    IplusNo = numpy.where(S_u == 1)[0];
    IMinusNo = numpy.where(S_u == 0)[0];

    #I+集合からランダムに一つサンプリング
    RandMax = IplusNo.shape[0];
    I_idx = numpy.random.randint(RandMax);
    i = IplusNo[I_idx]

    #I-集合からランダムに一つサンプリング
    RandMax = IMinusNo.shape[0];
    J_idx = numpy.random.randint(RandMax);
    j = IMinusNo[J_idx]

    return u , i , j

GetBprOpt
なお、論文通りのBPRoptを出力すると、
かなりでかい値になって、収束判定しづらくなるので、
最後に足し合わせた要素数（ソース内のZu）で割るようにしています。

def GetBprOpt(S , W , H , Lam):
    Bpropt = 0.0;
    Zu = 0.0;

    for u  in range(W.shape[0]):
        #ユーザuの全アーティストに対する嗜好度算出
        SC =numpy.dot(W[u , :] , H.T);

        #好みのアイテムの嗜好度のみ抽出
        IpScore = SC[ S[u,:] == 1 ];

        #未観測のアイテムの嗜好度のみ抽出
        ImScore = SC[ S[u,:] == 0 ];

        #(i , j)全ての組み合わせを計算し足しあわせる
        AA = numpy.tile(IpScore , (ImScore.shape[0] , 1))
        BB = numpy.tile(ImScore , (IpScore.shape[0] , 1))
        X=(AA-BB.T);
        lnSigX = numpy.log( 1/( 1 + numpy.exp(-X) ) );
        Bpropt = Bpropt + ( sum( sum(lnSigX) ) );

        #足し合わせた要素数も計算（最後にこれで割る）
        Zu = Zu + float(ImScore.shape[0] * IpScore.shape[0])

    Bpropt = (Bpropt - Lam * ((numpy.linalg.norm(W)**2)
                         + (numpy.linalg.norm(H)**2)))/Zu;
    return Bpropt

以上が今回組んだBPR MFのソースコードとなります。
BPROptの計算は全組み合わせを1度に計算するんじゃなくて、
$W$ と $H$ の更新時にBPROptも更新すれば速かったんじゃないかと、
いまさら思ってるんですが・・・まぁ、いいや

さて、いよいよ実験です。
まずパラメータ設定はk=50、alpha=0.7、Lam=0.1としました。
そして、このプログラムに上述の実験用データ突っ込み一晩寝かせます。
すると・・・

350万回繰り返したところで収束してました。
その結果ですが、先ほどの疑似ユーザ君3名に対する
推薦結果上位10アーティストは以下の通りです。

K-ON	gotham	joging
soft ballet	christian death	2 many dj's
ali project	coil	coil
the yellow monkey	the sisters of mercy	johann sebastian bach
perfume	the gathering	spiritualized
メリー	saltatio mortis	the fall
gackt	melvins	melvins
capsule	bauhaus	colour haze
電気グルーヴ	front 242	melt-banana
depapepe	ulver	guided by voices
lula	summoning	roxy music

遊びなんで定量的な評価はしませんが、
個人的には結構うまくいってるんじゃないかなと思います。というか、かなり満足。

以降、考察と言うより自分の趣味の観点からの感想です。

まず、ユーザK-ONに対する推薦結果ですが、
あの6組からcapsuleが推薦されたのは、個人的には結構嬉しかったりします。
イエモンが上位にきてるのも何となくわかる気が・・・ｗ
あと、20位くらいまで見るとsyrup16gやbuck-tickなんかもいて、なかなか面白い。
まぁ、これ見てると、テクノとかロックとかジャンル的な要素よりも、
日本のアーティストであるということが、推薦に最も寄与しているような気がします。
隠れ変数50個の中に、「日本」という属性を示す要素があるかもしれませんね。

つぎに、gothamくんの結果をみてみると、予想通りゴス、インダストリアル、
ブラックメタル系のアーティストが上位に来ています。うん、暗い。怖い。
このメンツが上位にいるなら、Throbbing GristleやPsychic TVなんかも
上位にいてもよさそうですが、残念ながら100位圏外でした。
20位まで見てみると、jesuやelectric wizardも挙がっていてました。
好みのアーティストとして挙げたのは70年代後半〜80年代の
いわゆるニューウェーブ期のバンドでしたが、
推薦結果には、それらに影響を受けた90〜00年代の
アーティストもいて、なかなか面白い。
年代、国関係なしに、その手の系統のアーティストが推薦されるといった感じでしょうか。

最後にユーザjoging。
これに関しては、本当にこの結果を利用して、
新たなジョギングソングを発掘しようと目論んでいました。
で、結果をみて最初に目についたのが、推薦結果の第3位。
まさかのバッハ。
確かにG線上のアリア聞きながら走るのも、
なかなか乙なものがあるのかもしれません。
他にも、monoという日本のポストロックバンドを知ることができたりと、
かなり使えそうです。

2012-01-09

情報推薦の論文

機械学習

明けましておめでとうございます。

ブログ開設当初は、月1ペースの更新を目標にしてたんですが、
案の定、2ヶ月以上放置。
気がつけば年が明けてました・・・
こんな感じで、今年もテキトーにやっていこうと思うので、
よろしくお願いします。

というわけで、今回は情報推薦です。
amazonの「この商品を買った人はこんな商品も買っています」や、
聴いた曲の履歴からお勧めのアーティストを推薦してくれるLastFM等々
世の中、情報推薦サービスが増えてきています。
研究トピックとしても結構盛り上がってる分野で、
ACM主催のrecsysのように情報推薦の国際会議もあります。
学生のころ所属していた研究室でも、情報推薦の研究をしている人がいて、
評価実験の被験者になったりもしました。
こういう評価実験の被験者をやってても、何が推薦されるか結構楽しみだったりするので、
自分もちょっと挑戦してみたい分野でした。

といった興味も、最近は忘れかけていたのですが、
先日こんなwikiを発見して、情報推薦熱がちょっとだけ再燃。
ということで、適当に最近の論文をチラ見してみたり、
リファレンス論文を追ってみたりして、以下の論文を発見しました。

S. Rendle, C. Freudenthaler, Z. Gantner, and L. Schmidt-Thieme
BPR: Bayesian Personalized Ranking from Implicit Feedback (UAI2009)

いつものように、細かい説明や手法の正当性なんかの説明は論文に委ねるとして、
ここでは、手法の大まかな考え方と流れだけ述べていこうと思います。
タイトルにある「Implicit Feedback」とは、
クリック情報や商品を購入したという情報などのような、
数値的には0以上の値であるユーザの振る舞いの情報と言えます。
以降は、ユーザがある商品を購入した/しないという情報を基に推薦する
システムを例に挙げて、話を進めていきます。

まず、従来のMatrix Factorizationに基づく推薦手法を述べます。
ユーザ $u$ の商品 $i$ に対する嗜好度を $\hat{x}_{uj}$ 、
また、 $\hat{X}$ を $\hat{x}_{ui}$ 要素に持つ $\| U \| \times \| I \|$ の行列とします。
なお、ユーザの集合を $U$ 、商品の集合を $I$ としています。
更に、 $W$ を $\| U \| \times k$ の行列、 $H$ を $\| I \| \times k$ と定義し、
$\hat{X} = WH^T$
とします。 $\hat{x}_{ui} =<{\bf w}_u , {\bf h}_i>$ です。
既に分かっているユーザの嗜好度情報 $\hat{x}_{uj}$ を基に、
$\hat{X}$ を近似する $W$ と $H$ を求めて、推薦を行おうというのが、
従来のMatrix Factorizationによる推薦手法のようです。
正則化なしの最小二乗近似の場合、特異値分解により $W$ と $H$ を得られますが、
オーバフィッティングにより、うまくいかないようです。
そもそも、この論文のタスクImplicit Feedbackの場合、
商品を購入した/していないの情報だけしかないので、
ある商品に対するユーザの嗜好度 $\hat{x}_{ui}$ が分かっているわけではありません。
（嗜好度はユーザの商品に対する評価（amazonで言う星の数）などのExplicit Feedbackといえます。）

というわけで、この論文の提案手法です。
まず、この論文では「既に購入した商品は、購入していない商品に比べて嗜好度合いは高い」
という仮定を置き、モデルパラメータをMAP推定する手法となっています。
この手法で肝となるのは、アイテム（商品）のペアがあり、ユーザはそのうちのどちらが好みであるか
という情報を学習することです。
「商品のペアのうち、どちらが好みであるか」というのは、仮定より
ユーザが既に購入している商品の方が好みで、購入していない方は好みでないということになります。
両方の商品を購入済みの場合は、嗜好度に差をつけません。
従って、このペアに関しては学習に用いません。
逆に、商品のペアのうち両方とも購入していないとき、
このユーザはどちらの商品が好みであるか、
これを予測して最終的には推薦まで持っていくわけですね。

ユーザ $u$ が購入した商品の集合を $I_u^+$ とします。
ここで、ユーザ $u$ は、商品 $i$ の方が商品 $j$ より好みだ、
と言う三つ組を集めた集合を $Ds$ とします。
$Ds = \left\{ \left(u , i , j\right) \| i \in I_u^+ \wedge j \in I \backslash I_u^+ \right\}$ です。
この $Ds$ からモデルパラメータを推定します。

この論文では、求めたいパラメータ $\Theta$ の事後分布を $p \left( \Theta \| >_u \right)$ としています。
$>_u$ は、ユーザ $u$ の好みの構造を表しているようです。
事後分布は尤度と事前分布の積に比例するので、
$p \left( \Theta \| >_u \right) \propto p \left( >_u \| \Theta \right) p \left( \Theta \right)$
となります。

ユーザの好みは、それぞれ独立としています。従って、
$\displaystyle p \left( >_u \| \Theta \right) = \prod_{\left(u , i , j\right) \in Ds} p \left( i >_u j \| \Theta \right)$
となります（途中計算は省いているので、詳細は論文の方をご覧ください）。
ここで、シグモイド関数 $\sigma \left(x \right) = 1/\left(1 +e^{-x} \right)$ を用いて、
$p \left( i >_u j \| \Theta \right) = \sigma\left(\hat{x}_{uij} \left( \Theta \right) \right)$
となります。 $\hat{x}_{uij} \left( \Theta \right)$ は、どのようなモデルを適用するかによって変わる関数です。
具体的なモデルについては、Matrix Factorizationに基づくものと、
k-nearest neighborの2つが論文内で示されています。
具体的なモデルはちょっと置いておいて、今度は $p \left( \Theta \right)$ を考えます。
こちらは、おなじみ（？）平均0、共分散行列 $\lambda {\bf I}$ の正規分布としています。
ということで、目的関数は、
$\displaystyle BPR-Opt = \max_{\Theta} \ln p \left( >_u \| \Theta \right) p \left( \Theta \right)$
です。これと上記の式から、以下の最適化問題を得ます。

$\displaystyle BPR-Opt = \max_{\Theta} \sum_{\left(u , i , j\right) \in Ds} \ln \sigma\left(\hat{x}_{uij} \left( \Theta \right) \right) - \lambda || \Theta ||^2$

この最適化問題を解くことで、パラメータを推定し、推薦を行うというのが
この論文で提案しているBayesian Personalized Ranking（BPR）の枠組みとなります。

$BPR-Opt$ の微分を計算すると以下のようになります。

$\displaystyle \frac{\partial}{\partial \Theta} BPR-Opt \propto \sum_{\left(u , i , j\right) \in Ds} \frac{-e^{-\hat{x}_{uij}}}{1+e^{-\hat{x}_{uij}}} \cdot \frac{\partial}{\partial \Theta} \hat{x}_{uij} - \lambda \Theta$
（2012/2/3追記
　自分で計算してみたら、 $e$ の前のマイナスが出ませんでした。
　詳しくは、この記事の一番下参照）

ここでは、 $\hat{x}_{uij}(\Theta)$ を略して $\hat{x}_{uij}$ と書いています。
通常の最急勾配法は、
$\Theta \leftarrow \Theta + \alpha \frac{\partial}{\partial \Theta} BPR-Opt$
という更新式になるのですが、
一般的に $\left(u , i , j\right) \in Ds$ は膨大な数であり、
その数だけの合計を更新のたびに計算しなければいけないのは、
計算量や収束性の面でかなり難があります。

そこで、最急勾配法ではなく、以下のような確率的勾配法を用いて最適解を推定します。

1． $\left(u , i , j\right) \in Ds$ から1組をランダムにサンプリングする。
2．この選択した1組を用いて、以下式で $\Theta$ を更新する。
$\Theta \leftarrow \Theta + \alpha \left( \frac{-e^{-\hat{x}_{uij}}}{1+e^{-\hat{x}_{uij}}} \cdot \frac{\partial}{\partial \Theta} \hat{x}_{uij} - \lambda \Theta \right)$
（2012/2/3 追記
　ごめんなさい！！ここも間違えて書いちゃってます。
　詳しくは下記参照。）

収束するまで、これを繰り返すのが、BPRの学習アルゴリズムです。

ここまで来ると、じゃあ $\Theta$ と $\frac{\partial}{\partial \Theta} \hat{x}_{uij}$ は何ぞやとなってきます。
その具体的な形の１つがMatrix Factorizationです。
まず、Matrix Factorizationで求めたいのは、
$\hat{X}$ を $WH^T$ の形で近似する $W$ と $H$ です。
従って、 $\Theta = \left(W , H \right)$ です。
次に、 $\hat{x}_{uij} \left( \Theta \right)$ は、
$\hat{x}_{uij} \left( \Theta \right) = \hat{x}_{ui} \left( \Theta \right) - \hat{x}_{uj} \left( \Theta \right) = <{\bf w}_u , {\bf h}_i> - <{\bf w}_u , {\bf h}_j>$
とみなせます。
これにより、 $\frac{\partial}{\partial \Theta} \hat{x}_{uij}$ は次のように場合分けができます。

サンプリングにより、 $\left( u , i , j \right)$ を選択したとして、

${\bf w}_u$ の更新
$\frac{\partial}{\partial \Theta} \hat{x}_{uij} = {\bf h}_i - {\bf h}_j$

${\bf h}_i$ の更新
$\frac{\partial}{\partial \Theta} \hat{x}_{uij} = {\bf w}_u$

${\bf h}_j$ の更新
$\frac{\partial}{\partial \Theta} \hat{x}_{uij} = - {\bf w}_u$

それ以外
$\frac{\partial}{\partial \Theta} \hat{x}_{uij} = {\bf 0}$

つまり、 $\left(W , H \right)$ のうち、１回のステップで更新されるのは、
サンプリングされた $\left( u , i , j \right)$ に対応する ${\bf w}_u$ , ${\bf h}_i$ , ${\bf h}_j$ の計３行のみです。

これがBPRによるMatrix Factorizationアルゴリズムですが、
正直、本当に収束するんかいな？って感じです。
確かめるためには、実装して試してみるのが１番なんで、
早いうちに実装したいなと思いますが、
とりあえず、今日のところはここまで。
この実装結果は、別の記事に書きたいと思います。
なるべく今年中に・・・

2012/2/3追記
収束しない！？
と思ったら、いくつか間違いを発見。

まず、 $BPR-Opt$ 微分ですが、
自分で計算してみたら、
$\displaystyle \frac{\partial}{\partial \Theta} BPR-Opt \propto \sum_{\left(u , i , j\right) \in Ds} \frac{e^{-\hat{x}_{uij}}}{1+e^{-\hat{x}_{uij}}} \cdot \frac{\partial}{\partial \Theta} \hat{x}_{uij} - \lambda \Theta$
になりました。
　
ということで、更新式は
$\Theta \leftarrow \Theta + \alpha \left( \frac{e^{-\hat{x}_{uij}}}{1+e^{-\hat{x}_{uij}}} \cdot \frac{\partial}{\partial \Theta} \hat{x}_{uij} - \lambda \Theta \right)$
になる。
と思いきや、論文をもう一度読み直してみると、
$\Theta \leftarrow \Theta + \alpha \left( \frac{e^{-\hat{x}_{uij}}}{1+e^{-\hat{x}_{uij}}} \cdot \frac{\partial}{\partial \Theta} \hat{x}_{uij} + \lambda \Theta \right)$
と、正則化項のところがプラスになってます。

自分の計算が間違えてるのかなと思いましたが、
同じ研究グループの別の論文
Learning Attribute-to-Feature Mappings for Cold-Start Recommendations (ICDM2010)
内でもBPRのアルゴリズムが記述されており、
そちらでは、更新式が
$\Theta \leftarrow \Theta + \alpha \left( \frac{e^{-\hat{x}_{uij}}}{1+e^{-\hat{x}_{uij}}} \cdot \frac{\partial}{\partial \Theta} \hat{x}_{uij} - \lambda \Theta \right)$
になってるので、これが正解のようです。

2011-10-17

Eigenvalue Sensitive Feature Selection

機械学習

クラスタリングや分類などを行う際、なるべく冗長となる特徴量は取り除いてから
行った方が精度や実行時間の観点で有効だとされています。
この「冗長となる特徴量を取り除く」方法は特徴選択と呼ばれ、
これまでに多数の手法が提案されています。
ということで今回は、新たな特徴選択手法を提案している論文

Yi Jiang , Jiangtao Ren
Eigenvalue Sensitive Feature Selection (ICML2011)

を読んでみました（以下、この論文の提案手法をESFSとします）。

具体的な方針としては、ある特徴を変化させた（取り除いた）とき、
Laplacian Eigenmaps（Belkinらによって提案されたグラフラプラシアンを利用した次元削減手法。以下、LE）
の結果にどれだけ影響を及ぼすか計ろうというものです。
例えば $i$ 番目の特徴の変化が、LEの結果に多大な影響を及ぼすようなら、
その特徴は重要な特徴であり、
その逆に、あまり変化を与えないならば、その特徴は冗長な特徴 -> 取り除いてOK
という考え方です。

細かいことは置いておいて、早速、ESFSのアルゴリズムを。

と、その前に・・・

1つのサンプルの特徴ベクトルを ${\bf x}=(x_1 , x_2 , \cdots , x_K)^T$ と表します。
従って、その次元数は $K$ です。
サンプルは今、 $n$ 個あるとします。

すいません。
今度こそ、本当にアルゴリズムです。

1．各特徴に重み ${\bf w}=(w_1 , w_2 , \cdots , w_K)^T$ をかけたときの、
類似度行列を $S({\bf w})$ を以下のようにします。

$\displaystyle S({\bf w})_{ij} = \exp\left(-\frac{\sum_{t=1}^{K}(w_t({x}^i_t - {x}^j_t))^2}{2\sigma^2}\right)$

実際、アルゴリズム上では ${\bf w}=(1,1,\cdots,1)^T$ のときの $S({\bf w})$ が必要になります。
つまり、ここではLE等でおなじみの「普通の」類似度行列 $S$ を計算しておきます。

2．一般化固有値問題 $L({\bf w}){\bf q}({\bf w})_r=\lambda_r({\bf w})D({\bf w}){\bf q}({\bf w})_r$ を解きます。
$D({\bf w})$ は、 $D({\bf w})_{ii} = \sum_{j}S({\bf w})_{ij}$ となる要素を持つ対角行列、
$L({\bf w})$ は、 $L({\bf w})=D({\bf w})-S({\bf w})$ です。
これまた全変数が ${\bf w}$ の関数となっていますが（ $S$ が ${\bf w}$ の関数であるため）、
実装上は、 ${\bf w}$ は考える必要ありません。
従って、ここまではLEとほぼ同じ操作を行うことになります。

3．各特徴のeigenvalue sensitive criterion（EVSC）を求めます。
このEVSCが、各特徴のLEの結果に与える影響、
つまり、特徴の重要度となります。
式の導出過程等の解説は論文の方に譲って、
ここでは、最終的な計算式のみ記述します。

$\displaystyle {\rm EVSC}(t)=\sum_{r=1}^n\|\frac{\partial \lambda({\bf w})_r}{\partial w_t}|{\bf w}=\bf{1}\|$

$\displaystyle \frac{\partial \lambda({\bf w})_r}{\partial w_t}=\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}g(i,j,t)(q({\bf w})_{ri}-q({\bf w})_{rj})^2$
$\displaystyle \hspace{5em} - \sum_{i=1}^{n-1}\sum_{j=i+1}^{n}g(i,j,t)(\lambda({\bf w})_r(q({\bf w})_{ri}^2+q({\bf w})_{rj}^2))$

$g(i,j,t)=-w_t\frac{(x_t^i - x_t^j)^2}{\sigma^2}S({\bf w})_{ij}$

ここで、 $q({\bf w})_{ri}$ は固有ベクトル ${\bf q}({\bf w})_{r}$ の $i$ 番目の要素を表します。
ちょっとだけ解説すると、 $t$ 番目の特徴量の重み $w_t$ を $1 \rightarrow 0$ へ変化させ、
それ以外の重みは $w_i=1 (i \neq t)$ と固定した場合の $\lambda({\bf w})_r$ の変化率は、
$\|\frac{\partial \lambda({\bf w})_r}{\partial w_t}|{\bf w}=\bf{1}\|$ で近似されます。
従って、全固有値についてその変化率を計算し合計すれば、
$t$ 番目の特徴の重要度が分かるという寸法です。
正直、この辺りの理解が微妙で、ちょっと納得できていない部分があるんですが・・・

何はともあれ、あとはEVSC値が高い特徴のみ残し、低い特徴は取り除けばOK。

実装して適当な疑似データ作って遊んでみました。
疑似データは、3次元のガウス分布 $N_1({\bf \mu}_1 = (2 , -1 , 0)^T , {\bf I})$ に従うサンプル100個、
ガウス分布 $N_2({\bf \mu}_2 = (-2 , 1 , 0)^T , {\bf I})$ に従うサンプル100個の計200個のデータです。

こんな感じ

やってみると、算出されるEVSC値には結構差が出ます。
1つ目の特徴が最も大きく、2つ目の特徴に比べ3倍くらいの値、
3つ目の特徴のスコアが最も小さく、2つ目の特徴に比べ1/3くらいの値になります。
そのため、上位2つの特徴のみ使用するとなった場合、
3つ目の特徴量が取り除かれることになります。
クラスタリング等を考慮すると、それっぽい結果なので、
うまくいっているんじゃないかなと思います。
満足満足。

とはいえ、特徴選択手法は疑似データでやっても、ちょっと面白みに欠けるかなと。
個人的には、特徴選択は実データで行って、
「この特徴とこの特徴の組み合わせが残るのかぁ、意外だな〜」などと、
発見する使い方が面白いと思います。

あと、python（+numpy、scipyモジュール）で実装したんですが、
馬鹿正直に上記の式をfor文使いまくって実装すると、やっぱりかなり遅いです。
高速化のためには、for文で計算している部分を行列計算になおすなど、
繰り返し構文をなるべく使わない工夫が必要になります。
実際にどのように実装したかは、気が向いたら書こうと思います。

2011-09-08

Principal Anglesの日本語訳って何？

と、タイトルからいきなり疑問形なわけですが・・・
先日の記事で書いたGrassman Discriminant Analysisの論文内に
このPrincipal Anglesについて記述されていたので疑問に思ったわけです。
線形代数をきちんと学んだ人にとっては常識なのかもしれませんが、
お行儀があまりよろしくない学生だった僕には初耳の単語です。

さて、このPrincipal Anglesの定義は以下の通りです
（記号・表記の定義は前回の記事をご覧ください）。
2つの部分空間の間に下式を満たす $0 \leq \theta_1 \leq \theta_2 \leq \cdots \leq \theta_m \leq \pi/2$ が定義され、
これをPrincipal Anglesという。

$\displaystyle \cos{\theta_k}=\max_{{\bf u}_k} \: \max_{{\bf v}_k} \: {\bf u}_k^T {\bf v}_k$

${\rm subject} \: {\rm to}$
${\bf u}_k\in span({\bf Y}_1)$ , ${\bf v}_k\in span({\bf Y}_2)$ ,
${\bf u}_k^T{\bf u}_k=1$ , ${\bf v}_k^T{\bf v}_k=1$
${\bf u}_k^T{\bf u}_i=0$ , ${\bf v}_k^T{\bf v}_i=0 \; (k \neq i)$

つまり、 $\theta_1$ に関しては、なす角が最小（ $\cos{\theta}$ が最大）になるように、
両部分空間からそれぞれベクトル（ただしノルムが1）を持ってきて、
そのなす角を $\theta_1$ としましょう。
同様に $\theta_2$ も両部分空間からベクトルを選び、そのなす角を $\theta_2$ とするが、
ベクトルは $\theta_1$ の時に選んだそれと直交するように選びましょう。
ということです。 $\theta_3$ 以降も同様です。

では、このPrincipal Anglesを得る方法ですが、
実は、 ${\bf Y}_1^T{\bf Y}_2$ の特異値が $\cos\theta_1,\cos\theta_2,\cdots,\cos\theta_m$ と書かれています。
・・・なぜ？？？

ということで、本当にそうなのか確かめてみたいと思います。
とはいえ、全部のPrincipal Anglesが特異値だと示すのはちょっと面倒なので、
ここでは、Principal Anglesの１つは ${\bf Y}_1^T{\bf Y}_2$ の特異値ですよと示す方針で行きたいと思います。

まず、 $span({\bf Y}_1)$ 及び $span({\bf Y}_2)$ の任意の要素は、
任意のベクトル ${\bf a},{\bf b} \in {\bf {\cal R}}^m$ を用いて、
それぞれ以下のように表すことができます。

${\bf Y}_1 {\bf a} \in span({\bf Y}_1)$
${\bf Y}_2 {\bf b} \in span({\bf Y}_2)$

これらをPrincipal Anglesの定義式に代入すると、

$\displaystyle \cos{\theta}=\max_{{\bf a}} \: \max_{{\bf b}} \: {\bf a}^T {\bf Y}_1^T {\bf Y}_2 {\bf b}$ 　　（１）

${\rm subject} \:{\rm to}$
${\bf a}^T {\bf a}=1$ , ${\bf b}^T {\bf b}=1$

となります（今、１つのPrincipal Angleしか考えないので、直交の制約は省略）。
この最適化問題をラグランジュの未定乗数法を解くために、
ラグランジュ乗数 $\lambda_1, \lambda_2$ を導入し、次の式を考えます。

$L= {\bf a}^T {\bf Y}_1^T {\bf Y}_2 {\bf b} - \lambda_1({\bf a}^T {\bf a}-1) - \lambda_2({\bf b}^T {\bf b}-1)$

これを ${\bf a}$ 及び ${\bf b}$ で微分し、それぞれ0と置きます。
$\frac{\partial L}{\partial {\bf a}} = {\bf Y}_1^T {\bf Y}_2 {\bf b} - 2 \lambda_1 {\bf a}=0$
$\frac{\partial L}{\partial {\bf b}} = {\bf Y}_2^T {\bf Y}_1 {\bf a} - 2 \lambda_2 {\bf b}=0$

２つの式から、次の２式が導出されます。
${\bf Y}_1^T {\bf Y}_2 {\bf Y}_2^T {\bf Y}_1 {\bf a} = \lambda' {\bf a}$
${\bf Y}_2^T {\bf Y}_1 {\bf Y}_1^T {\bf Y}_2 {\bf b} = \lambda' {\bf b}$
ただし、 $\lambda' = 4\lambda_1\lambda_2$ 。

おお！これを見ると、 ${\bf a}$ と ${\bf b}$ の停留点はそれぞれ、
${\bf Y}_1^T {\bf Y}_2$ の（同じ特異値に対応する）左特異ベクトルと右特異ベクトルじゃないか！
ゴールが見えてきました。
ここで、 ${\bf a}$ と ${\bf b}$ の最適解（停留点の１つ）を $\hat{{\bf a}}$ と $\hat{{\bf b}}$ 、
${\bf Y}_1^T {\bf Y}_2$ の特異値分解を ${\bf Y}_1^T {\bf Y}_2={\bf U}{\bf \Sigma}{\bf V}^T$ とします。
これを（１）の式に代入すると、

$\cos{\theta}=\hat{{\bf a}}^T {\bf U}{\bf \Sigma}{\bf V}^T \hat{{\bf b}}=({\bf U}^T\hat{{\bf a}})^T {\bf \Sigma}{\bf V}^T \hat{{\bf b}}$
となります。
$\hat{{\bf a}}$ と $\hat{{\bf b}}$ はそれぞれ左右の特異ベクトルであることを思い出すと、
${\bf U}^T \hat{{\bf a}}={\bf V}^T \hat{\bf{b}}$ は１つの要素のみ1、他の要素は全て0のベクトルだとわかります。
従って、 $\cos{\theta}$ は ${\bf Y}_1^T{\bf Y}_2$ の特異値であるといえます。
ちょっと怪しい気がしますが・・・これで大丈夫なはず。

このPrincipal Anglesは、Grassmann Discriminant Analysis含め様々な
サンプルが部分空間である機械学習手法で要となっています。
例えば、Face recognition using temporal image sequence[Yamaguchi98]の論文では、
$\sin{\theta_1}$ を距離として用い、最近傍法を行うことで、
部分空間集合の学習・識別を実現しているらしい（元論文は読んでいません）。

最後にPrincipal Anglesとグラスマン多様体の関係についてちょっとだけ述べておくと、
グラスマン多様体上の２点間の測地線距離 $d_G({\bf Y}_i , {\bf Y}_j)$ は

$\displaystyle d_G^2({\bf Y}_i , {\bf Y}_j) = \sum_{k=1}^{m} \theta_k^2$

になるそうです。
多様体については、もっとちゃんと勉強せねば・・・と修士１年のころから思いつつ、
いつの間にか、社会人２年目なわけですよ

2011-09-02

Grassmann Discriminant Analysisの論文を読む

機械学習

ときどき、学生のころ勉強していたこと（機械学習など）が懐かしくなって、
それに関連しそうな論文を探したり読んだりするんですが、
今日はその論文の1つを紹介しようかなと思います。

J. Hamm and D. D. Lee.
Grassmann Discriminant Analysis: a Unifying View on Subspace-Based Learning (2008)

従来の教師あり学習だと特徴ベクトルとして表現されているデータを
学習・判別するというのが一般的だといえます。
しかし、この論文では、複数の特徴ベクトルをひとまとめにした集合を、
1つのデータとして学習・判別する方法について言及しています。
例えば、同じ被写体を異なる角度から写した複数の画像を入力とし、
その被写体が何であるか判別するといったタスクへの応用が考えられます。

「複数の特徴ベクトルをひとまとめにした集合」をこの論文では、部分空間として取り扱います。
つまり、 $N$ 個の学習サンプルがそれぞれ $m$ 次元部分空間として表現されているときの、
教師あり学習手法を提案しています。
取り扱う問題をもう少し詳しく記述すると、
まず、特徴ベクトルの次元を $D$ 、各サンプルの部分空間の次元は $m$ とします（当然、 $D>m$ です）。
${\bf Y}\in{\bf {\cal R}}^{D\times m}$ を部分空間の正規直交基底を並べた直交行列とし、
${\bf Y}$ で張られる部分空間を $span({\bf Y})$ と表現します。
従って、 $N$ 個の学習サンプルの集合は、 $\{span({\bf Y}_1),span({\bf Y}_2),\cdots,span({\bf Y}_N)\}$ になります。
なお、このような部分空間を元とする集合をグラスマン多様体 ${\bf {\cal g}}(m,D)$ と呼ぶそうです。
なので、グラスマン多様体上の $N$ 個の学習サンプルをどう学習するかという問題になります。

グラスマン多様体とかなにやら物騒な単語が出てきました。
と言っても、その辺りの理論は置いていおいて実装したいというなら、
やることは部分空間同士のカーネルを定めて、
カーネル拡張したフィッシャー判別分析法を適用するというものなのでさほど難しくないと思います。
てことで、早速実装方法に移ります。

まず、特徴ベクトルの集合 $\{{\bf x}_1,{\bf x}_2,\cdots,{\bf x}_p\}$ から部分空間の直交基底を計算して $\bf{Y}$ を得ます。
その方法としては、 $\bf{X}=[{\bf x}_1,{\bf x}_2,\cdots,{\bf x}_p]$ という行列を特異値分解して、
その左特異ベクトルのうち $m$ 個を ${\bf Y}$ とすればよいはずです。
この操作を $N$ 組の学習サンプル全てに施し、 $\{{\bf Y}_1,{\bf Y}_2,\cdots,{\bf Y}_N\}$ を得ます。

次に、グラム行列を計算します。
この論文で提案しているカーネル関数は以下の2種類です。

Projection kernel
$k_p({\bf Y}_i,{\bf Y}_j) = ||{\bf Y}_i^T{\bf Y}_j||_{F}^2$

Binet-Cauchy kernel
$k_{bc}({\bf Y}_i,{\bf Y}_j) = (\det{\bf Y}_i^T {\bf Y}_j)^2$

このどちらかのカーネル関数を用いてグラム行列を計算した後は、
カーネルフィッシャー判別分析を適用し次元削減すればOK。
カーネルフィッシャー判別分析については、この論文内にも記述されていますし、
Fisher Discriminant Analysis with Kernels[Mika99]が元論文のはずなので、
そちらも参考になると思います。
次元削減後は普通にユークリッド空間上のデータとして扱えます。
煮るなり焼くなりしちゃってください。

とりあえず、目的と手法の実装方法を書いていみましたが、
上記2つのカーネル関数がどこから出てきたのかとか、グラスマン多様体がどうたらとか
肝心なことをスルーしちゃってます。
その辺りのことは、また別の記事で書くかもしれません（面倒になって放置する可能性の方が高いです）。
なお、ここまで書いてきましたが、実はまだ自分で実装して試していません。
これまた、いずれ実装してその結果をブログに書ければなぁと思います。

甲斐性なしのブログ

うどんくらいしか食べる気がしない

BPRで遊んだよ！

情報推薦の論文

Eigenvalue Sensitive Feature Selection

Principal Anglesの日本語訳って何？

Grassmann Discriminant Analysisの論文を読む