2012-06-02

NIPS2011斜め読み part2

機械学習

ICML2012もアクセプトされた論文のタイトルとアブストが公開され、
いよいよ盛り上がってまいりましたが、でも僕はNIPS2011！
懲りずに斜め読み第2弾いきます。
間違って理解している可能性も大いにあり得るので、
それを発見した際には、指摘していただければ幸いです。

Heavy-tailed Distances for Gradient Based Image Descriptors

コンピュータビジョン分野で、SIFTやHOGなどの勾配ベースの特徴量を利用して、
対応点マッチングやbag of keypointsの構築が行われますが、
その際、特徴量同士の比較には大抵ユークリッド距離が用いられます。
ユークリッド距離はノイズがガウス分布の場合に対応する距離ですが、
SIFTのような勾配ベースの特徴量は、ガウス分布に従わないよ！と主張し、
Gamma compound-Laplace (GCL)分布という新たな分布と、
それに基づく距離の計算方法を提案しています。
まず、GCL分布ですが、
$p(x \| \mu;\alpha , \beta)=\frac{1}{2}\alpha \beta^{\alpha}(\|x - \mu \| + \beta)^{- \alpha - 1}$
という式で定義され、高尖度かつheavy-tailedな特性を持つ分布になります
（1次元の確率分布であることに注意）。
ここで、 $\mu$ は特徴のプロトタイプ（サンプル $x$ は $\mu$ にノイズが乗って生成されるという考え方？）、
$\alpha,\beta$ は特性を制御する超パラメータです。
この特性は、SIFT特徴量の分布にフィットするようです。
じゃあ、この分布に基づく距離はどのように定義するか？ですが、
これには、尤度比検定で使われる尤度比を利用します。
具体的には、2つの独立したサンプル $x,y$ が与えられた時、
「2つのサンプルの距離が小さい」->「この2つは同じ分布から生成された」場合は、
$s(x,y)=\frac{p(x \| \hat{\mu}_{xy})p(y \| \hat{\mu}_{xy})}{p(x \| \hat{\mu}_{x})p(y \| \hat{\mu}_{y})}$
（ただし、 $\hat{\mu}_{xy}$ は2つのサンプル $\{x,y\}$ からの最尤推定解、
$\hat{\mu}_{x}$ と $\hat{\mu}_{y}$ はそれぞれ1つのサンプル $\{x\}$ 、 $\{y\}$ から推定される最尤推定解）
という尤度比が大きくなるはずです。
従って、この尤度比に基づき、1つの次元の距離を
$d(x,y)=\sqrt{-\log{(s(x,y))}}$
と定義しています。
これをGCL分布に当てはめると、
$\hat{\mu}_{xy}=x$ または $\hat{\mu}_{xy}=y$ （どちらを採用しても結果は同じ）、
$\hat{\mu}_{x}=x$ 、 $\hat{\mu}_{y}=y$ なので、結果的に距離は、
$d^2(x,y)=(\alpha + 1)(\log{( \| x - y \| + \beta)} - \log{\beta})$
となります。なお、 $\alpha,\beta$ はサンプル集合から最尤推定で求めます。
これは、1次元上の距離なので、多次元ベクトルの場合は、
各次元でこの距離の2乗を算出し、足し合わせてベクトル同士の距離とします。

確率分布を基に距離を定義するため、尤度比を利用するやり方は、
これ以外にも色々応用できそうです。

Efficient Methods for Overlapping Group Lasso

Group Lassoやproximal operatorなど僕が知らなかった&アツそうなトピック満載の論文。
正直、あまり理解できていません。難しい…
まず、Group Lassoですが、それぞれの特徴がいくつかのグループにわけられていた時、
そのグループ毎に正則化を行うことで、グループ単位でパラメータが0になる解が
得られるという手法です。式で表すと、
$\displaystyle \min_{\bf{x} }\hspace{1ex} l(\bf{x})+ \phi_{\lambda_2}^{\lambda_1}(\bf{x})$ 　ただし、 $\phi_{\lambda_2}^{\lambda_1}(\bf{x})=\lambda_1 \parallel \bf{x} \parallel_1 + \lambda_2 \sum_{i=1}^g w_i \parallel \bf{x}_{G_i} \parallel_2$
( $l$ は連続で凸な損失関数、 $\bf{x}_{G_i}$ はグループ $i$ に属する特徴のベクトル、
$w_i$ はグループ毎の重みで、実験ではそのグループに所属する特徴数の平方根を採用)
を解く問題になります。
このGroup Lassoの最適化問題ですが、特徴がそれぞれ1つのグループに所属している場合に比べ、
所属するグループがオーバーラップしている場合、解くのがより難しい問題になるそうです。
これを解くために、この論文では、accelerated gradient descentという最適化の手法を
適用することを提案しています。
手法としては、
$f_{L, {\bf x} } ({\bf y}) = l( {\bf x} )+{<} l'( {\bf x} ) , {\bf y}-{\bf x}\ {>} + \phi_{\lambda_2}^{\lambda_1}({\bf y})+ \frac{L}{2} \parallel {\bf y} - {\bf x} \parallel_2^2$
とし、 ${\bf x}_{i+1} = \arg \min_{{\bf y}}f_{L_i,{\bf s}_i}({\bf y})$ という計算を繰り返していくというものです。
ここで、 $\bf{s}_i={\bf x}_i + \beta_i ({\bf x}_i - {\bf x}_{i-1})$ であり、 $L_i$ と $\beta_i$ は毎回適切なものを探索する必要があります。
さて、
$\displaystyle \pi_{\lambda_2}^{\lambda_1}({\bf v}) = \arg \min_{{\bf x}} \hspace{1ex} \frac{1}{2} \parallel {\bf x} - {\bf v} \parallel_2^2 + \phi_{\lambda_2}^{\lambda_1}({\bf x})$
とすると、毎ステップの最適化式は、次のように変換できます。
${\bf x}_{i+1} = \pi_{\lambda_2/L_i}^{\lambda_1/L_i} \hspace{1ex} ({\bf s}_i - \frac{1}{L_i}l'({\bf s}_i))$
この、 $\pi_{\lambda_2}^{\lambda_1}({\bf v})$ は、proximal operatorと呼ばれ、いろいろと便利な性質を有しています。
例えば、 $\lambda_2=0$ （つまり、通常のlasso）なら、
proximal operatorの最適化問題の解は、解析的に求まります。
この論文でも、group lassoの問題におけるproximal operatorの性質を解析し、
そこから、問題のサイズを小さくするための前処理方法の提案や
その双対問題を解いて更に効率化できると主張しています。
Proximal Operatorを用いた最適化手法は、
損失関数が微分可能かつ凸な関数で、正則化項が凸だが微分不可能という場合に
有効なようで、調べてみると、「解くべき最適化問題のサイズが大きすぎて困っちゃうわ」->
「そんなときは、これ！Proximal Operator」->
「よーし、この問題におけるProximal Operatorを解析しちゃうぞ」
という流れの論文をちらほら見かけました。
また、どちらかというと、Group Lassoという言葉に惹かれてこの論文を読んだので、
そちらに関しても、関連文献をサーベイしたいですね。

Transfer Learning by Borrowing Examples for Multiclass Object Detection

その、Group Lassoを転位学習に応用した論文。
物体検出を目的とした論文で、その検出器の学習を対象となるクラスのサンプル（画像）だけでなく、
他の似たクラスのデータも利用しましょうという手法です。
例えば、ソファの検出器を作る際、色々なクラスの画像データセットを与えると、
ソファクラスの画像だけでなく、肘掛け椅子などの似た画像のクラスも
学習に利用して検出器を構築します。
どのクラスのデータをどれだけの重みで利用するかということも学習します。
今、n枚の画像サンプルがあり、それぞれ $\cal C=\{1 \cdots C \}$ 中の
いずれかのクラスが割り当てられているとします。
加えて、バックグラウンドの画像サンプルもb枚用意し、それには-1と言うクラスを割り当てます。
従って、n+b枚の画像サンプル $({\bf x}_i , y_i)$ ( $i=1,\cdots,n+b$ 、 ${\bf x}_i$ は特徴ベクトル、 $y_i$ はクラスラベル)から
検出対象の検出器の識別関数を学習することになります。
ということで、クラス $c$ における検出器の目的関数は以下の通り。
$\displaystyle \sum_{c \in \cal C} \min_{{\bf \beta}^c} \hspace{1ex} \min_{{\bf w}^{*,c}} \sum_{i=1}^{n+b}(1-w_i^{*,c}) \hspace{2ex} l(<{\bf \beta}^c , {\bf x}_i> , \rm sign(y_i)) + \lambda R({\bf \beta}^c) + \Omega_{\lambda1 , \lambda2}({\bf w}^{*,c})$

( $\min$ の左の $\sum$ は必要なんでしょうか？)
ここで、 ${\bf \beta}^c$ は検出器の識別関数の重みベクトル、
${\bf w}^{*,c}$ はどのサンプルをどの程度の重みで利用するか
決めるベクトル（各要素が $0 \leq w_i^{*,c} \leq 1$ ）です。
従って、 $w_i^{*,c}$ が小さければ小さいほど、そのサンプルの重みは高くなります。
また、検出対象であるクラスcに属するサンプルとバックグラウンドである
クラス-1に属するサンプルの重みは $w_i^{*,c}=0$ という制約を設けます。
（なお、論文内では、 ${\bf w}^{c}=1-\bf{w}^{*,c}$ を重みとして話を進めているのですが、
式の議論を明瞭にするため（?）、目的関数内では ${\bf w}^{*,c}$ を使用しているようです。）
また、 $R(\cdot)$ は ${\bf \beta}^c$ は損失関数に対する正則化項、
$\Omega_{\lambda1 , \lambda2}(\cdot)$ は先ほどみたGroup Lassoの正則化項です。
ここで言うgroupはクラスです。従って、Group Lassoの効果により、
クラスごとにサンプルを学習に利用するか否かの決定がなされます。
あとは、 ${\bf \beta}^c , {\bf w}^{*,c}$ を得るために目的関数の最適化を行う必要がありますが、
それには、一方を固定して一方を最適化する操作を繰り返すという手法をとっています。
この論文のさわりとしてはこんな感じだと思います。
他にも、ソファのクラスの検出器学習に車クラスを利用しなかった場合は、
車クラスの検出器学習にソファクラスを利用しないという制限を加える方法や
対象クラスに適合するための画像変換方法も示されていますが、ここでは言及しません。

しかし、なぜ ${\bf w}^{c}$ のまま式を構築せず、わざわざ ${\bf w}^{*,c}$ を使ったんでしょう？
Group Lassoの正則化項を $\Omega_{\lambda1 , \lambda2}({\bf w}^{c})$ にすれば、
対象に似たクラスのサンプルには高い重みが与えられ、それ以外の重みは0になると思うので、
こちらの方がシンプルだし、理にかなっているような気がするんですが・・・

と言ったところで今回はお開き。
最近、論文読んでばっかりで、全然実装してないですね。

2012-05-20

NIPS2011斜め読み

機械学習

NIPSは毎年12月頃に開催される機械学習関連の国際会議です。
Proceedingsはweb上に公開されているので、
今回は昨年末に開催されたNIPS2011の中から、
興味を惹かれた論文を何本か適当にチョイスして読んでみました。
斜め読みした程度なので、かなり理解が甘いです。
でも、自分用メモ書としてブログに書いちゃう。

Sparse Filtering

タイトルのシンプルさに惹かれ、思わず保存した論文。
M個のサンプル ${\bf x}^{(i)} \hspace{1ex} (i=1 \cdots M)$ が与えられて、
教師なし学習でそれぞれのサンプルをより良い
特徴ベクトル ${\bf f}^{(i)}$ にしてあげましょうというものです。
（この枠組みがUnsupervised Feature Learningと言われるそうですが、
次元削減とかとは違うのでしょうか？）
で、各サンプルの特徴ベクトルの要素 $f_{j}^{(i)}$ は、
$f_{j}^{(i)}={\bf w}_j^T {\bf x}^{(i)}$ という線形和で得るんですが、
なるべく、得られる特徴ベクトルはスパースにしたい（0要素を多くしたい）。
もう少し言うと、特徴抽出後のM個のサンプル ${\bf f}^{(i)}$ は
Population Sparsity、Lifetime Sparsity、High Dispersalという
3つの性質を有した集合にしたい。
そのためには、
${\rm minimize} \sum_{i=1}^M \parallel \hat{{\bf f}}^{(i)}\parallel_1$
ただし、
$\hat{{\bf f}}^{(i)} = \frac{\tilde{{\bf f}}^{(i)}}{\parallel \tilde{{\bf f}}^{(i)}\parallel_2}$
$\tilde{{\bf f}}^{(i)} = \frac{{{\bf f}}^{(i)}}{\parallel{{\bf f}}^{(i)}\parallel_2}$
を解けばいいのよーというものです。
超パラメータがすくないぜ（抽出後の次元数だけ？）、
実装楽チンだぜ、収束速いぜという3拍子がそろってるらしいので、
試してみたいですね。

Learning a Distance Metric from a Network

ネットワークの解析にリンク情報だけでなく、
各ノードの特徴ベクトルも考慮しましょうという論文。
例えば、twitterなんかのSNSを利用していると、おすすめのユーザが提示されますが、
その推薦のために、リンク関係（twitterだとフォロー関係）だけでなく、
各ユーザの特徴ベクトル（twitterでは、ツイートから構築される特徴ベクトル）も
考慮すれば、より洗練された推薦ができるようです。
じゃあ、具体的にどうするのって話ですが、これには計量学習の枠組みを利用します。
これは特徴ベクトル集合 ${\bf X}$ とネットワークの接続行列 ${\bf A}$ から
特徴ベクトル同士の距離の計算方法 $D_M ({\bf x}_i,{\bf x}_j )$ を学習するというものです。
この距離は、半正値行列 ${\bf M}$ を用いて、
$D_M ({\bf x}_i,{\bf x}_j )=({\bf x}_i-{\bf x}_j )^T {\bf M}({\bf x}_i-{\bf x}_j )$
として表されるので、 ${\bf X}$ と ${\bf A}$ から最適な ${\bf M}$ を求めることになります。
さて、どのような ${\bf M}$ を求めればよいのかですが、
${\bf X}$ の各特徴ベクトル同士の（ $D_M$ で計算される）距離関係のみから
得られる接続行列（例えば、あるサンプルのk近傍のサンプルは接続しているとみなす）が
実際のネットワークの接続行列 ${\bf A}$ と一致するように、 ${\bf M}$ を求めるというものです。
その（特徴ベクトルの距離関係から計算される）接続行列を得る方法として、
k-nearest neighborを使う場合と、maximum weight subgraphを使う場合の
2つの場合の最適化アルゴリズムが、論文内に示されていますが、
その辺りはまだよくわかりません。

Metric Learning with Multiple Kernels

計量学習をマルチカーネル化する手法を示した論文。
多くの計量学習はこの論文に示されている方法で、
マルチカーネルに拡張できそうです。
まず、マルチカーネル学習ですが、
$k_{\mu}({\bf x}_i , {\bf x}_j)=\sum_{l=1}^{m} \mu_lk_l({\bf x}_i , {\bf x}_j) , \mu_l \geq 0 , \sum_{l=1}^{m} \mu_l = 1$
のように、 $m$ 個あるカーネル関数を重みつきで足し合わせることで、
新たなカーネル関数を作るための重み $\mu_l$ を学習する枠組みです。
このカーネル関数 $k_{\mu}({\bf x}_i , {\bf x}_j)$ で、内積が定義される空間を $\cal H_{\mu}$ とし、
元の特徴ベクトル ${\bf x}$ を $\cal H_{\mu}$ に写像したベクトルを ${\bf \Phi}_{\mu}( {\bf x} )$ とするとその空間上での計量距離は、
$d_{{\bf M},\mu}^2 ({\bf \Phi}_{\mu}( {\bf x}_i ) , {\bf \Phi}_{\mu}( {\bf x}_j )) = ({\bf \Phi}_{\mu}( {\bf x}_i ) - {\bf \Phi}_{\mu}( {\bf x}_j ) )^T {\bf M} ({\bf \Phi}_{\mu}( {\bf x}_i ) - {\bf \Phi}_{\mu}( {\bf x}_j ) )$
となるので、学習サンプルから最適な半正値行列 ${\bf M}$ とカーネルの重み $\mu_l$ を求めるのが、
マルチカーネルの計量学習となります。
多くのカーネル空間上の計量学習アルゴリズムは、
$\displaystyle \min_{{\bf M}_l}F_h({\bf M}_l) \hspace{1ex} s.t. \hspace{1ex} C_h(d_{{\bf M}_l}^2 ({\bf \Phi}_l( {\bf x}_i ) , {\bf \Phi}_l( {\bf x}_j ))) , {\bf M}_l \succeq 0$
という最適化問題を解いていて、その最適解は $\eta_h {\bf I}+{\bf \Phi}_l( {\bf X} )^T {\bf A}_l {\bf \Phi}_l( {\bf X} )$ になるようです
（ ${\bf \Phi}_l( {\bf X} )^T$ は、各行が写像後の学習サンプルの行列）。
しかも、多くの計量学習アルゴリズムは $\eta_h = 0$ になります。
そして、マルチカーネル化したいアルゴリズムの最適解が ${\bf \Phi}_l( {\bf X} )^T {\bf A}_l {\bf \Phi}_l( {\bf X} )$ という形なら、
マルチカーネル化した時の、最適解も ${\bf \Phi}_{\mu}( {\bf X} )^T {\bf A} {\bf \Phi}_{\mu}( {\bf X} )$ となるそうです。
あとは、これを行列 ${\bf A}$ を固定して $\mu$ を求めて、次に $\mu$ を固定して ${\bf A}$ を求める
という操作を収束するまで繰り返します。
このまま、 $\mu_l$ と ${\bf A}$ を求めるのが、 $ML_h-MKL_{\mu}$ と論文内で呼んでいるもので、
これは、採用する計量学習アルゴリズムが凸なら、 ${\bf A}$ を固定して $\mu_l$ の最適解を解く問題も
凸になりますよという手法です。
また、 ${\bf P}={\bf \mu}{\bf \mu}^T$ として、最適化式をちょっと変形した $ML_h-MKL_{P}$ も
提案しているのですが、これは、計量学習アルゴリズムが ${\bf P}$ に対して線形で、
カーネル数が少ない時には、 ${\bf A}$ を固定して ${\bf P}$ を求める問題が、
等価な凸最適化問題に変換できるというものです。
更に、 $\mu$ を半正値行列 ${\bf A}''$ 内にひっくるめてしまうことで、
その ${\bf A}''$ のみを推定する問題になるので、
既存の計量学習アルゴリズムをそのまま利用できる
$NR-ML_h-MKL$ というものも提案しています。
マルチカーネル学習+計量学習がこれまで提案されていなかったって言うのがちょっと意外でした。

Trace Lasso: a trace norm regularization for correlated designs

回帰問題などで正則化項として $l1$ ノルムを利用すると、
解がスパースが得られやすくなります。
これはlassoと呼ばれ、多くの研究で特徴選択や回帰・分類問題に利用されています。
ところが、特徴間に強い相関がある場合、lassoはあまりよろしくない結果をもたらします。
ということで、特徴間の相関に適応する新たな正則化項を提案したのがこの論文です。
まず、学習サンプル ${\bf x} \in \cal R^p$ とし、n個のサンプルをまとめた行列を
${\bf X}=({\bf x}_1,\cdots,{\bf x}_n)^T \in \cal {R}^{n \times p}$ とします。
一般に線形の教師あり学習は、
$\displaystyle \hat{{\bf w}}=\arg \min_{{\bf w}} \frac{1}{n}\sum_{i=1}^{n}\hspace{1ex}l(y_i,{\bf w}^T{\bf x}_i)+\lambda f({\bf w})$ 　（ $y_i$ は学習のラベル、 $l$ は損失関数）
を解くことになるんですが、ここで $f({\bf w})$ として
$l2$ ノルムの2乗 $\parallel {\bf w}\parallel_2^2$ を採用すればリッジ回帰、
$l1$ ノルム $\parallel {\bf w}\parallel_1$ を採用すればlassoになります。
このとき、 $\bf{X}$ の各列が正規化（たぶん、 $l2$ ノルムを1にする正規化）されていれば、
$\parallel {\bf w}\parallel_2=\parallel {\bf X} \rm{diag}({\bf w}) \parallel_F$ , $\parallel {\bf w}\parallel_1=\parallel {\bf X} \rm{diag}({\bf w}) \parallel_{2,1}$ と変形できます。
（ $\parallel \cdot \parallel_F$ はフロベニウスノルム、 $\parallel \cdot \parallel_{2,1}$ は各列の $l2$ ノルムの和）
リッジ回帰もlassoもノルムが異なるだけで、その中身は共通しています。
そこで、このノルムをトレースノルム $\parallel {\bf X} \rm{diag}({\bf w}) \parallel_*$ を採用しよう
というのが、この論文で提案しているトレースlassoです。
（トレースノルムは $\parallel{\bf A} \parallel_* =\rm{trace}( \sqrt{ {\bf A}^T {\bf A}})$ というもので、 ${\bf A}$ の特異値の和になります。）
このトレースlassoですが、 ${\bf X}$ の各列同士が直交してる（各特徴量同士無相関）なら、
通常のlassoと等価になります。
一方、 ${\bf X}$ の各列が全て同じベクトルの場合、これはリッジ回帰と等価になります。
従って、トレースlassoは学習サンプルにおける特徴間の相関が強ければリッジ回帰、
弱ければlassoに近いふるまいをするといえます。
また、もっと一般的に各列が単位ベクトルである適当な行列 ${\bf P}$ を用いて、
$\parallel {\bf P} \rm{diag}({\bf w}) \parallel_*$ という正則化項も提案しています。
この ${\bf P}$ 次第で、lassoもリッジ回帰も今話題の（?）グループlassoも記述できるようです。

今回はここまで。
全体を通して感じた印象としては、

やっぱり皆スパースがお好き
数年前に比べて、Multi ViewとかMulti Taskとかの論文が少なくなった？
どの論文読んでも、最適化や近似のテクニックが眩しすぎる
グループlassoうまそう
タイトルがシンプルな論文が多い

などなど

あと、自分のチョイスが偏りすぎてちょっとよろしくないのかなって気がします。
しかも、チョイスしたのはいいけど、読んでてちんぷんかんぷんの論文や、
パッと見ただけでそっと閉じた論文も多々あり、
自分の知識や能力、モチベーション的な何かの欠如を改めて実感しました。
研究を行う機関に所属していない自分にとっては、
NIPSやICMLなど、ただでProceedingsをweb公開してくれる国際会議は本当にありがたいです。

2012-04-17

エンベデッドシステムスペシャリスト試験

お久しぶりです。
4月15日、エンベデッドシステムスペシャリスト試験を受験してきました。
ブログの更新が滞ってたのは、これの試験勉強に励んでいたため・・・
というわけではないです。ただ単に面倒くさかっただけです。

受験の理由ですが、仕事で工場設備等のソフトウェア設計に従事してるにも拘らず、
組込みシステムに関する基礎知識が欠如していると感じたので、
とりあえず、資格試験勉強で知識を賄おうと思ったのがきっかけです。
ちなみに設備といっても、マイコンのソフトウェアをやってるわけではなく、
Windowsが入ってるPCが設備に組み込まれてて、
そいつがボード経由で設備の制御をしたり、
通信で電子機器の制御をしたりというシステムが主なので、
Windowsアプリの開発が多く、
ときどきリアルタイムOS上で動くソフトの開発もやるって感じです。

あと、学生の頃の友人達の多くが学生時代に高度を取ってて、
「えー、高度持ってないのが許されるのは小学生までだよねー」
なんてことを言われたか言われてないかは置いといて、
まぁ、周りに流されたってのも受験の1つの理由です
（僕は学生の頃、IPA関連の資格は何一つ持ってませんでしたが・・・）。

で、手ごたえですが

午前1
去年、応用情報を取得したので免除でした。

午前2
解答がその日のうちに提示されたので、
答え合わせをしてみると25問中20問正解。
6割正解していれば通過なので、まぁ大丈夫でしょう。
でも、単語から伝わるニュアンスや消去法、完全なる勘で
答えた問題がかなりあったので、
解答が提示されるまで、ここで落ちたと思ってました。

午後1
午後からは記述式。
大問が3題あって、問1は必須、問2と問3はどちらか選択する方式です。
例年、問2はソフトウェア系の問題、問3はハードウェア系の問題になってるようなので、
僕は、問3をみることなく問2を選択しました。
問1はレーザ加工機、問2はスマートアダプタシステムがテーマで、
おなじみの通信速度に関する問題やタスク間、装置間の
メッセージ送受信に関する問題がありました。
とりあえず、全問埋めることができ、それなりにできた感じはするんですが、
過去問をやったときも、自信があっても採点してみると
ボロボロってことがよくあったので、
本番はそうなってはいないことを祈るばかりです。

午後2
大問が2題で、そのうち1問選択する形式です。
問1はハードウェア、問2はソフトウェアなので問2を選択。
今年は介護用電動ベッドシステムの問題でした。
こちらも、午後1と同様、手ごたえとしてはまあまあ、
でも、過去問の経験上受かってる気はあまりしないってとこです。

全体的に今回の試験は、過去に比べると計算問題が少ないような気がします。
2桁以上の数字の掛け算だけでも萎えてしまう僕にとっては、
その点ラッキーでした。

リアルタイムOSのタスク設計やメッセージのやり取り、
装置間の通信制御辺りは、自分の仕事にも直結するので、
試験勉強してても、「なるほど、確かにそういう設計すべきだよな」とか
「あー、こうやっちゃたらこういう問題が起きちゃうわけか」
みたいに参考になる部分が多くあったので合否にかかわらず、
受験しておいてよかったんじゃないかなと思います。

あと、関係ないですがインフルエンザにかかってしまいました。
熱が37度程度の微熱でもインフルエンザってことがあるんですね。
食欲も旺盛で、こうしてブログ書くくらいの元気はあるんですが・・・
会社には医者の許可が下りるまで来るなと言われているので、
今週は全部お休みをいただくことになったのですが、
正直、仕事が結構ピンチなので、不安で不安で仕方ありません。

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$

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

甲斐性なしのブログ

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

NIPS2011斜め読み part2

Heavy-tailed Distances for Gradient Based Image Descriptors

Efficient Methods for Overlapping Group Lasso

Transfer Learning by Borrowing Examples for Multiclass Object Detection

NIPS2011斜め読み

Sparse Filtering

Learning a Distance Metric from a Network

Metric Learning with Multiple Kernels

Trace Lasso: a trace norm regularization for correlated designs

エンベデッドシステムスペシャリスト試験

BPRで遊んだよ！

情報推薦の論文

Eigenvalue Sensitive Feature Selection

Principal Anglesの日本語訳って何？