甲斐性なしのブログ

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

l1ノルム正則化でスパースな離散フーリエ変換

機械学習 スパース学習

はじめに

離散フーリエ変換は離散化された時間軸の信号を周波数軸に写像する基底変換の一種で、周波数解析などによく使われます。この離散フーリエ変換の係数は、基底ベクトルの重み付き線形和で表したときの誤差の2乗和を最小にする重みと見なせます。つまり、離散フーリエ係数は一種の線形回帰問題の解として得られます。ということは、l1ノルム正則化項つけたLasso回帰として解けばスパースな離散フーリエ変換ができるのでは？ということでやってみました。

離散フーリエ変換の線形回帰的定式化

元の信号のデータ点数をT個、基底ベクトルの個数をF個とし、上述の通り回帰問題として離散フーリエ係数を求める問題を定式化すると下記のようになります。

\displaystyle \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 \tag{1}

ここで、{\bf y}は元の時間軸の信号で、{\bf y} = \left[ y(0), y(1), \cdots, y(T - 1) \right]^Tです。
{\bf W} = \left[{\bf w}_0, {\bf w}_1, \cdots, {\bf w}_{F - 1} \right]はF個の基底ベクトルを劣ベクトルとして並べたT \times Fの行列です。
一般的に離散フーリエ変換では{\bf W}は{\bf w}_k = \left[\exp(i\frac{2\pi k \times (0)}{T}) ,\exp(i\frac{2\pi k \times (1)}{T}), \cdots, \exp(i\frac{2\pi k \times (T - 1)}{T} ) \right]^Tという複素数の要素を持つ行列で、各列ベクトル同士は直交します({{\bf w}_i^{*} {\bf w}_j} = 0 \quad i \neq j)。
そして、\bf{x}が求めたい離散フーリエ係数であり、{\bf x} = \left[ x(0), x(1), \cdots, x(F - 1) \right]^Tです。
通常、F = Tです。そのため、(1)の最適解は{\hat{\bf x}} = {\bf W}^{-1} {\bf y} = {\bf W}^{*} {\bf y}で表され、最小2乗誤差は0になります。

l1ノルムによるスパース化

ということで、回帰問題としてみなせるならl1ノルムの正則化項をつける（つまりは、Lasso回帰を行う）とスパースな解が得られるだろうという考えです。(1)の目的関数に正則化項を付け加えます。

\displaystyle \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 + \lambda \| {\bf x} \|_1 \tag{2}

こうなると最早この問題は解析的には解けず、正則化項は\|{\bf x} \|_1 \leq \gammaという制約条件を付けるのと等価なので、最適解による最小2乗誤差も一般的には0にはなりません。ということで、cvxpyというソルバを用いて上記の最適化問題を解きます。
なお、この定式化は入口こそ離散フーリエ変換を回帰問題と見なしてLasso回帰を行うという考えですが、結局のところこれは、この資料等に記載されいるノイズあり圧縮センシングの一種です。　

複素数の最適化問題から実数のみの最適化問題へ

離散フーリエ変換は上述のように、基底ベクトルが複素数で得られる解も一般的に複素数となります。ところが、このままだと実装上面倒くさい（というより、ソルバが動いてくれるかわからない）ので実数のみの式に変換します。といっても、\exp({i \theta}) = \cos\theta + i\sin\thetaより実部と虚部に分けてそれぞれ個別の変数として考えるだけです。
具体的には、{\bf W} を{\bf W} = \left[ {\bf w}_0^c, {\bf w}_1^c, \cdots, {\bf w}_{T - 1}^c, {\bf w}_0^s, {\bf w}_1^s \cdots, {\bf w}_{T - 1}^s \right]というT \times 2Tの行列として再定義します。ここで、
{\bf w}_k^c = \left[\cos(\frac{2\pi k \times (0)}{T}) ,\cos(\frac{2\pi k \times (1)}{T}), \cdots, \cos(\frac{2\pi k \times (T - 1)}{T} )\right]^T
{\bf w}_k^s = \left[\sin(\frac{2\pi k \times (0)}{T}) ,\sin(\frac{2\pi k \times (1)}{T}), \cdots, \sin(\frac{2\pi k \times (T - 1)}{T} )\right]^T
です。あわせて説明変数{\bf x}を要素数2Tのベクトルとして式（２）を解けばよいです。この変換で基底ベクトルの数も説明変数の要素数も倍になったように見えますが、元々も複素数でそれぞれが実部と虚部を持っていたため本質的には変わりません。
また、離散フーリエ変換の場合、x(k) + {\bar x(k + \frac{T}{2}})という関係があるため、実質0\sim\frac{T}{2}-1番目の要素までを求めれば、残りの\frac{T}{2}\sim T-1番目は自動で決定されます。これは、l1ノルム正則化項をつけても変わりません（実際に試しました）。そのため、説明変数の要素数を2T\rightarrow Tへとシュリンクすることも可能ですが、今回は理解を深めるためにもそのまま解くようにしました。

cvxpyを用いた実装

上述の通り、今回は最適化問題を解くためにcvxpyというツールを利用しました。このcvxpyの使い方はこちらの記事を参考にしました。cvxpyの使い方は他の解説記事に譲るとして、早速式（２）を解く関数のソースコードを示します。

from cvxpy import *
import numpy as np
from numpy.random import *
import matplotlib.pyplot as plt

#y:原信号、T:信号のデータ点数、lam:ハイパーパラメータλ
def sparse_dft(y, T, lam):
    #行列Wを生成
    V = np.meshgrid(np.arange(T),np.zeros((T,1)))[0]
    u = np.arange(T)
    V = (V.T * u) * (2 * np.pi /T)
    W = np.c_[np.cos(V), np.sin(V)]

    #変数xを定義
    x = Variable(T * 2)

    #ハイパーパラメータを定義
    lamda = Parameter(sign="positive")
    lamda.value =lam

    #目的関数を定義
    objective = Minimize(sum_squares((y - W * x))/T + lamda*norm(x, 1))
    p = Problem(objective)

    #最適化計算
    result = p.solve()

    return result, x

特に特別なことはしておらず、最初に{\bf W}を生成して、後はcvxpyの解説に従って変数、ハイパーパラメータ、目的関数を定義し、p.solve()で問題を解きます。

実験

以下のようなコードで疑似データを発生させて実験してみました。

    T = 2**10
    y = (3 * np.cos(28* 2*np.pi * np.arange(T)/T) +2 * np.sin(28* 2*np.pi * np.arange(T)/T) +
            np.cos(400* 2*np.pi * np.arange(T)/T) + np.sin(125 * 2*np.pi * np.arange(T)/T) + 5 * randn(T) )

いくつかの周波数が異なるsin波とcos波を重み付きで重ね合わせて、更に正規分布に従うノイズを乗せた1024点のデータです。疑似データの波形は以下の通り。
f:id:YamagenSakam:20180129000732p:plain

通常の離散フーリエ変換

まずは上記の疑似データにFFTにより離散フーリエ変換を行いました。

    Y = np.fft.fft(y)/T #スケーリングを合わすため要素数で割る
    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-1.5, 1.5))
    plt.plot(Y.real , "r")
    plt.hold(True)
    plt.plot(Y.imag , "b")

結果は下記の図です。
f:id:YamagenSakam:20180129000744p:plain
重ねたsin波、cos波に対応する周波数成分が大きな値になるのは当たり前ですが、ノイズが重畳しているせいで全周波数成分まんべんなく値が入ります。

Lassoによる離散フーリエ変換

次は今回記した手法での離散フーリエ変換（\lambda = 0.5で実験）。上述のsparse_dft関数で計算します。

    (result, x) = sparse_dft(y, T, 0.5)

    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-1.5, 1.5))
    plt.plot(x.value[0:T,0], "r")
    plt.hold(True)
    plt.plot(x.value[T:T*2,0], "b")

f:id:YamagenSakam:20180129000758p:plain
通常の離散フーリエ変換とは異なりスパースな解が得られ、余分な周波数成分は0になっていることがわかります。

逆変換

得られたフーリエ係数から逆変換を行い元の信号を再構築してみました。こうすることで、余計なノイズ成分が除去された信号が得られるはずです。
なお、逆変換のため、得られた{\bf x}から\cosにかかる係数を実部、\sinにかかる係数を虚部として複素数に変換してから、ifftメソッドで逆変換を行っています。

    X=(x.value[0:T] + (x.value[T:T*2] * 1j))*T
    y_hat = np.fft.ifft(X.T)

    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-20, 20))
    plt.plot(y_hat.real.T, "r")

以下が原波形から今回の方法でノイズを除去したデータです。
f:id:YamagenSakam:20180129000809p:plain
また、以下が上述の疑似データ生成のコードからノイズ項を取り除いたデータの波形です（つまり、真に得たい波形）。
f:id:YamagenSakam:20180129002851p:plain
確かにノイズが除去された波形が得られています。実際のノイズ項を除いた信号と比べると若干スケーリングが小さくなっていますが、悪くない結果だと思います。

まとめ

離散フーリエ変換を回帰問題と見なし、目的関数にl1ノルムの正則化項をつけて誤差の最小化を行うことで、スパースな離散フーリエ係数が得られました。またそれを逆変換することで原信号からノイズを除去した信号を得ることもできました。また、最適化問題を解くためにcvxpyというツールを利用しました。
今回は何も考えずcvxpyに問題を突っ込んで解いただけでしたが、今後は、Lasso回帰の最適化等で使われる近接勾配法やcvxpyの詳しい使い方、凸解析など最適化の理論、実装両面で理解を深めていきたいと思います。

はてなダイアリーからはてなブログへ移行

このたび、本ブログをはてなダイアリーからはてなブログへ移行しました。旧はてなダイアリーの記事にアクセスすると、対応する本ブログの記事へのリダイレクトされるようなので、ブックマーク登録等の変更は不要です。
さて、この移行作業に当たって、元のダイアリーで書いていた数式がうまく表示されず修正に手間取った箇所がいくつかあったので、備忘録の意味でもここに記しておきます。（大体は元の書き方がよろしくなかっただけですが）

小中大括弧でエラー

  • 元々の書き方
[tex:\(x + y \)]
[tex:\{x + y \}]
[tex:\left{x + y \right}]
[tex:\[x + y \]]
[tex:\left[x + y \right]]

見え方
\(x + y \) ：おかしい。
\{x + y \}：うまく表示される。
\left{x + y \right}：おかしい。
[x + y ]：うまく表示される。
\left[x + y \right]：おかしい。

小括弧は\エスケープだけではエラーになり、中大括弧は\left, \rightをつけるとエラーになります。そのため、小括弧には\left、\rightをつける、中大括弧は\left\{ \right}\}のように\left、\rightの後、更に\でエスケープする方法で対策しました。

  • 修正後の書き方：
[tex:\left(x + y \right)]
[tex:\{x + y \}]
[tex:\left\{x + y \right\}]
[tex:\[x + y \]]
[tex:\left\[x + y \right\]]

見え方
\left(x + y \right)：\left, \right追加でOKになった。
\{x + y \}：変更なし。
\left\{x + y \right\}：\エスケープ追加でOKになった。
[x + y ]：変更なし。
\left[x + y \right]：\エスケープ追加でOKになった。

sum、min、max等の上置き,下置きが上付き、下付きになる

  • 元々の書き方
[tex:\sum_{i = 1}^{N}x_i]
[tex:\min_{x} f(x)]

見え方
\sum_{i = 1}^{N}x_i
\min_{x} f(x)

本当は\Sigmaや\minの真上、真下に添え字が来るようにしたいのに、右上、右下に来ています。これは\displaystyleにすることで解決しました。

  • 修正後の書き方
[tex:\displaystyle \sum_{i = 1}^{N}x_i]
[tex:\displaystyle \min_{x} f(x)]

見え方
\displaystyle \sum_{i = 1}^{N}x_i
\displaystyle \min_{x} f(x)

\bfにより意図しないところまでも太字化される

  • 元々の書き方
[tex:(\bf{x}_N)]

見え方
(\bf{x}_N)：xの先にある添え字のNや)まで太字化される。

わかりにくいですが、xだけ太字にしたいのに、それ以降の文字まで太字になっています。これは元々の書き方がそもそも間違っており、中括弧の位置を\bf{x}→{\bf x}と変更すると修正できます。

  • 修正後の書き方
[tex:({\bf x}_N)]

見え方
({\bf x}_N)：xのみが太字。


ということで、今回は以上。せっかくはてなブログに移行しましたが、次回更新は遠い未来になるかも知れません。

NIPS2017論文メモ

機械学習

お久しぶりです。2014年9月以来、実に3年ぶりの更新です。しばらく機械学習関連の勉強は滞っていたのですが、思うところがあって再開。とりあえず、昨年12月に開催されたNIPS2017のプロシーディングス中から面白そうな記事をピックアップして斜め読みしました。備忘録＆リハビリのためにも、記事を更新してみよう思った次第です。

Deep Subspace Clustering Networks

部分空間クラスタリング+オートエンコーダで非線形なクラスタリングを行う手法を提案。

部分空間クラスタリングは、同じクラスなら同じ（低次元）部分空間埋め込まれている（つまり、同じクラスに所属している別のサンプルの線形結合で表せる）という仮定の下クラスタリングを行う枠組みで、SSCやLRRが有名。サンプル集合を\bf{X}として、いずれの手法も\bf{X}=\bf{X}\bf{C}という制約条件のもと、線形結合の行列\bf{C}のスパース化や低ランク化を行う（前者はl1ノルムの最小化、後者はトレースノルムの最小化）。\bf{C}を求めたら|\bf{C}|+|\bf{C}^T|をaffinity matrixとしてスペクトラルクラスタリングを行うのが一般的（他にも色々やり方があるらしい）。

これらの方法では非線形なクラスタリングは行えない。そこで、オートエンコーダを利用しようというのがこの論文の提案。オートエンコーダはニューラルネットを入力と出力を同じデータとして学習させることで、中間層の出力\bf{z}が非線形な特徴抽出となるという手法。入力から特徴を抽出する中間層まで層をエンコーダ層、特徴を抽出する中間層から出力までの層をデコーダ層と呼ぶ。この論文では、このエンコーダ層とデコーダ層の間に\bf{Z}\bf{C}とするSelf-Expressive層を新たに設けることを提案(\bf{Z}はエンコーダ層で特徴抽出されたサンプルの集合)。これは、線形の全結合となる。この\bf{C}と他の層の重みを求めるため、下記の損失関数を最小化。
L(\Theta, {\bf C} ) = \frac{1}{2} \|{\bf X} - {\bf \hat{X}_{\Theta}} \|^2_{F}+\lambda_1 \|{\bf C} \|_p + \frac{\lambda_2}{2} \|{\bf Z_{\Theta_e}} - {\bf Z_{\Theta_e}} {\bf C} \|^2_{F} \quad {\rm s.t.} \quad {\rm diag}( {\bf C} )={\bf 0}
\Thetaはニューラルネットの重み係数集合、{\bf Z_{\Theta_e}} はエンコーダ層の出力集合（バッチ）を表す。ノルム項はl1とl2の両方を実験。

上述のように、Self-Expressive層は線形の全結合層とみなせるので、バックプロパゲーションにより一括で学習できる。この方法で\bf{C}を求めたら後は従来通りにスペクトラルクラスタリング。結局のところ、オートエンコーダによって特徴抽出したデータに対して部分空間クラスタリングを行うが、その部分空間クラスタリング用の\bf{C}もニューラルネット上に組み込んで、一緒に学習してしまおうという手法。

顔画像クラスタリングのベンチーマークデータセット（Extended Yale BとORL）や物体画像クラスタリングのベンチマークデータセット（COIL20とCOIL100）で実験。結果、他の手法（カーネル部分空間クラスタリングなど非線形クラスタリング手法含む）よりも高精度なクラスタリングが可能。ノルム項ついては、何故かl2ノルムを使った方が良い結果に。

Mixture-Rank Matrix Approximation for Collaborative Filtering

協調フィルタリングでよく使われるLow-rank matrix approximationにおいて、低ランクと高ランクの近似行列を混合することにより精度を上げた手法。
実際、レイティング数が少ないアイテムやユーザに対する近似行列は低ランクにした方が推定結果は良く、逆にレイティング数に対する近似行列は高ランクにした方が推定結果は良い。

従来の手法はユーザ・アイテムの実際のレイティングマトリックスを\bf{R}として、\bf{R}\approx \bf{\hat{R}} = \bf{U}\bf{V}^T（\bf{U}, \bf{V}ともにk列の行列）となる\bf{U}, \bf{V}を求め、i番目のユーザのj番目のアイテムに対する嗜好スコアを\bf{U}_i\bf{V}^T_j（それぞれ、\bf{U}のi行目、\bf{V}のj行目のベクトル）として推定する。PMFという手法では、この問題を確率的にモデル化してMAP推定により\bf{U}, \bf{V}を求める手法を提供している。ここでkは固定の値とするため、ユーザおよびアイテムのレイティング数にかかわらず、ランクkの近似行列を推定する。

本論文では、このPMFをベースに複数のランクkを重み付きで混合させたモデルを提案。具体的には下記の確率モデルを構築している。
\displaystyle\Pr({\bf R}|{\bf U},{\bf V},{\bf \alpha},{\bf \beta},\sigma^2)=\prod_{i=1}^m \prod_{j=1}^n \left[\sum_{k=1}^K \alpha_i^k \beta_j^k N(R_{i,j}|{\bf U}_i^k{{\bf V}_j^k}^T,\sigma^2) \right]^{ \mathbb{1}_{i,j}}
ここで、{\textbf{U}_i^k}はランクkの行列\textbf{U}^kのi行目（ユーザiの特徴ベクトル）、{\textbf{V}_j^k}はランクkの行列\textbf{V}^kのj行目（アイテムjの特徴ベクトル）を表す。また、\alpha_i^kはユーザiにおけるランクkモデルの重み、{\beta_j^k}はアイテムjにおけるランクkモデルの重み。これにより、ランクkのモデルがユーザi、アイテムjにマッチするなら\alpha_i^k,\beta_i^kは大きくなるはずであり、マッチしないなら\alpha_i^k,\beta_i^kは小さくなるはず（期待的には0になる）。これらの\bf{U},\bf{V},\bf{\alpha},{\bf \beta}をMAP推定により求めようというのが提案手法。

MAP推定のため事前分布を与える。\textbf{U}_i^k、\textbf{V}_j^kの事前分布はそれぞれ、平均\textbf{0}、共分散\sigma_U^2 \textbf{I}、\sigma_V^2 \textbf{I}のガウス分布としている（行列\textbf{U},\textbf{V}の事前分布はこれらの総乗）。また、{\bf \alpha}^k, {\bf \beta}^kはそれぞれスパースなベクトルにしたいので、\alpha_i^k,\beta_i^kの事前分布はラプラス分布とする（それぞれ、平均\mu_{\alpha},\mu_{\beta}、分散b_{\alpha},b_{\beta}というパラメータを持つ）。これらを積算することによって導き出される事後分布のlogにマイナスをつけた損失関数を最小化する。実際には、そのままの損失関数を最小化するのは難しいので、イェンセンの不等式で導き出される下限を最小化する。

\textbf{U},\textbf{V},{\bf \alpha},{\bf \beta}は\sigma,\sigma_U^2,\sigma_V^2,\mu_{\alpha},\mu_{\beta},b_{\alpha},b_{\beta}を固定して確率的勾配法で解く。一方、　\sigma,\sigma_U^2,\sigma_V^2,\mu_{\alpha},\mu_{\beta},b_{\alpha},b_{\beta}も学習したいので、今度は\textbf{U},\textbf{V},{\bf \alpha},{\bf \beta}を固定して最尤推定する（こちらは、解析的に求まる）。これらの操作を収束するまで繰り返す。

MovieLensのデータセットやNetflix Prizeのデータセットで実験。結果、kを固定するPMFに比べて高い推定精度を得られるし、他の手法と比べても有意な差が得られる。また、混合するモデルの個数（k個数）はたくさんかつ細かくした方が良いようだ（計算時間はかかるけど）。


以下2つの論文については、まだほとんど理解できていない。とりあえず、備忘録として残す。

Kernel Feature Selection via Conditional Covariance Minimization

カーネル法を使って特徴とラベル間の非線形関係をとらえて特徴抽出する手法。条件付き共分散作要素のトレースを基準にすることを提案。（相互）共分散作要素は2つのヒルベルト空間間の共分散を表しており、そこから条件付き共分散作要素が導き出される。
論文では条件付き共分散作要素のトレースがラベルとの独立性を評価する指標となると主張し、その理由を説明している。また、一般に特徴選択の最適化をまともに行おうとするとNP困難な問題になってしまうので、ここでは離散最適化から連続最適化へ問題を緩和している。また、それでも計算量的にしんどいので、いくつか計算量を抑えるための目的関数、制約条件の変換も行っている。

共分散作要素をもう少し理解してから再チャレンジが必要。今のところ共分散行列のヒルベルト空間への拡張と理解。日本語だとこの辺りが良さげ。

Unsupervised Image-to-Image Translation Networks

2つのドメインの画像セットがあり、一方のドメインの画像をもう一方のドメインの画像に変換する手法を提案した論文。例えば、低解像画像から高解像画像へ変換するタスクや、朝の風景画像を夜の風景画像に変換するタスクなどがあげられる。Unspervisedなので、学習用に2つのドメイン間で対応する画像のセットが与えられるわけではなく、ただ単にドメインAの画像セットとドメインBの画像セットが与えられるのみ。
これを実現するために2つのドメインのデータは同じlatent spaceを共有するという仮定を置く。この仮定をもとに、variational autoencoder（VAE）とgenerative adversarial network（GAN）によりドメイン変換を行う。

まず、VAEでは2つのドメインのインプットが同じlatent spaceに変換されるように、エンコーダ層を学習。この同じlatent spaceへ変換というのを実現するために、エンコーダ層の最後の数層とデコーダ層の最初の数層はドメイン間で重みを共有するという制約をつける（なぜこれで、うまくいくのかがイマイチわからない）。
次にGANは、adversarial discriminatorが画像生成（ドメイン変換）のネットワークによって生成された画像か実際の画像かを高精度に判定できるように学習させる。画像生成ネットワーク側はそれに負けじと高精度なドメイン変換を行う。
これら2つを両立させる目的関数（VAEとGANと制約項の和）を最適化することにより、ドメイン変換ネットワークを学習する。

VAEやGANは名前しか聞いたことなかったが、面白い枠組みだし応用範囲も広そう。まだまだブームは続く？

所感

久々にNIPSの論文を読んでみましたが、やはりニューラルネット関連が多いという印象。それでもニューラルネットばかりだけではなく、未だにカーネル法に関する論文も結構あるので（「kernel」で論文タイトルの検索かけると20件ヒット）まだまだこの分野も廃れていない、今後盛り返す可能性も十分ありと感じました。
あと、VAEやGANなど生成モデル+ニューラルネットの枠組みは、今後いろいろな分野に波及しそうな気がします（GAN+トピックモデルとかVAE+HMMとか、もう誰かやってる？）。

【Python】生存報告がてら昔書いたFFTのコードを晒す

お久しぶりです。ちゃんと生きてます。色々と立て込んでたり何だりでモチベーションが低下していました。このまま、フェードアウトしていくのもよろしくないと思いつつもネタがない・・・なんて思っていたところ、PCのデータを整理してたら、以前作成したFFT（高速フーリエ変換）のプログラムを発見。てことで、このソースコードを今回のネタにしようかなと思います。
なお、フーリエ変換やDFT（離散フーリエ変換）についての解説は割愛します。また、FFTについても多少は説明しますが、DFTからFFTアルゴリズムの導出過程など詳しい説明は省きます。また、データのサンプル数は2のべき乗個に限定します。

FFTのソースコード

ということで早速ソースコードです。まずは、モジュールのインポートとFFTの関数です。

import numpy as np
import matplotlib.pyplot as plt
import time

def fft(x):
    n = x.shape[0]
    X = np.zeros(n ,  dtype = 'complex')
    if n > 2:
        y = x[0:n:2] #(1)
        z = x[1:n:2]

        Y = fft(y) #(2)
        Z = fft(z)

        Wn = np.exp(-1j * (2 * np.pi * np.arange(0 , n/2))/float(n)) #(3)
        X[0:n/2] = Y + Wn * Z
        X[n/2:n] = Y - Wn * Z
        return X
    else:
        X[0] = x[0] + x[1] #(4)
        X[1] = x[0] - x[1]
        return X

モジュールはいつも通り、numpyをインポート。後は、後々の動作確認で使うモジュールです。
関数fftが今回の肝であるFFTを行う関数ですが、見ての通り、再帰を使って書いてます。通常、スタックオーバーフローとかを避けるために、ビットリバース処理なんかをして繰り返し構文で書くのですが、今回はそこまでやってません。こっちの方が簡単だし直感的なので再帰で書く方法を採用しています。
ソースコード内のコメントにあるように、ポイントは4つ。
まず、(1)で元の信号xから偶数番目だけ取り出した信号yと奇数番目だけ取り出した信号zを作ります。
(2)では、fft関数を再帰呼び出ししてyとzのFFTの結果Y,Zを得ています。このY,Zを利用すると元の信号xのFFTの結果Xは次の式で、表されます。
\left\{ \begin{array}{l} X(i) = Y(i) + \exp \left( -j \frac{2 \pi i}{n} \right) Z(i) \: (0 \leq i < \frac{n}{2}) \\ X(\frac{n}{2} + i) = Y(i) - \exp \left( -j \frac{2 \pi i}{n} \right) Z(i) \: (0 \leq i < \frac{n}{2}) \end{array} \right.

(3)はこの処理を行っています。n点のDFTをこの式のようにn/2点のDFTの和に書きなおすことができるってのが、FFTの勘所です。
そして、(1)、(2)の再帰呼び出し処理を繰り返していくと、最終的にn = 2に行きつきます（元のデータサンプル数が2のべき乗の場合）。ここが再帰の終着点で、2点のDFTは(4)のように凄く簡単な計算で求められます。
これだけです。numpy使えば簡単にFFTのプログラムが書けます。ただし、numpyには元々FFTを行うメソッドが用意されているので、普通はわざわざ書く必要はありません。処理時間もこっちの方が断然速いです。

動作確認

てことで、一応動作確認。

def main():
    #テスト用の信号生成
    n =2**10
    x = (3 * np.cos(250* 2*np.pi * np.arange(0,n)/n) +
         2 * np.sin(250* 2*np.pi * np.arange(0,n)/n) +
         np.cos(30* 2*np.pi * np.arange(0,n)/n) )

    tic = time.clock()
    X = fft(x) #FFTの計算
    toc = time.clock()
    print toc - tic #処理時間表示

    #以下、結果グラフ表示
    plt.figure()
    plt.plot(X.real , "r")
    plt.xlim((0, n))
    plt.xticks(np.arange(0, 1024 , 100))
    plt.hold(True)
    plt.plot(X.imag , "b")
    plt.show()

テスト用の信号は、データのサンプル数が1024個で、周波数250のcos波、sin波、周波数30のcos波をそれぞれ適当なゲインをかけて足し合わせたものを使います。で、結果はこちら。

赤が実数部、青が虚数部です。うん、それっぽい。numpyに用意されているfftメソッドも同様の結果を吐き出します。また処理時間は、1024点なら30[ms]程度で処理できます。しかし、既存のメソッドの方は0.2[ms]程度。速い・・・

まとめ

今回はFFTのプログラム（再帰版）について話しました。上述のように、FFTはループ処理で書くことの方が多いようで、僕が知っている参考書では、大抵そちらの方法のサンプルプログラムが記述されています。しかし、それらは全部C言語で書かれているものなので、他の言語のサンプルコードが載っているものだと、どうなのか気になるところです。
あと、今回はサンプル点数が2のべき乗個である場合という最も基本的なFFTのみ記述しましたが、サンプル点数が素数の場合のFFTアルゴリズムが提案されていたりと応用的なFFTアルゴリズムも数々あります。これらを調べてみるのも、面白いかもしれません。

関係ないですが、転職により愛知から大阪に引っ越ししてきました。新天地からの更新でした。

マルチ○○学習まとめ

機械学習

機械学習の分野では、マルチ○○学習という名付けられた枠組み・手法が色々提案されています。僕は、接頭辞が共通だと、すぐごっちゃになって何が何だか分からなくなってしまうので、ちょっと整理したいと思っていました。ということで、今回は「マルチカーネル」、「マルチビュー」、「マルチタスク」、「マルチラベル」、「マルチインスタンス」をメモ書き程度にまとめました。

マルチカーネル学習

複数の特徴ベクトル表現、または、類似の尺度が考えられるタスクに対し、複数のカーネル関数を組み合わせたカーネル学習を考えましょうという手法の枠組みです。例えば、画像データの場合、特徴ベクトルとしてbag-of-keypointsを構築したり、色情報を利用したりと、様々な類似尺度を考えることができます。そこで、それぞれの類似尺度に対応するカーネル関数を構築して組み合わせれば、ハッピーになれるだろうという考え方です。
多くの場合、次のようにカーネル関数の線形和を考えます。
k({\bf x} ,{\bf y}) = \sum_{i=1}^{N} \beta_i k_i({\bf x} ,{\bf y})（Nはカーネル関数の数）
このカーネル関数を学習の目的関数に組み込んで、最適な\betaを求めるのが、メジャーな手法といえるでしょう。必要なカーネル関数の重み（\beta）だけ大きくなるので、どれが学習に寄与する類似尺度かもわかるというメリットもあります。
有名なのは、大規模データをマルチカーネル+SVM等で学習する手法を提案したLarge Scale Multiple Kernel Learningという論文でしょうか。また、Dimensionality Reduction for Data in Multiple Feature Representationsでは、マルチカーネルを用いた次元削減手法を提案しています。他にも「multiple kernel learning」とかで検索すれば様々な論文が引っかかり、Multiple Kernel Learning Algorithmsというサーベイ論文なども見つかるので、勉強はしやすいと思います。

マルチビュー学習

1つのサンプルが様々な特徴量（ビュー）の組み合わせで構成されるデータセットを学習する問題の枠組み。例えば、動画データは画像シーケンスと音声信号という2つのビュー組み合わせとみなせますし、画像が含まれるウェブサイトなんかはテキストデータ+画像+リンク情報というように考えることができます。このようなデータセットに対し、効果的な学習手法を考えるのがマルチビュー学習の枠組みです。
手法としては、Combining labeled and unlabeled data with co-trainingで提案された共訓練（co-training）と呼ばれる手法が代表的でしょうか。これは、元々半教師あり学習の一手法を提案した論文で、マルチビューとみなせるデータセットならば（ただし、ビュー間はラベル情報のもと条件付き独立である必要あり）、各ビュー別々に分類器を学習し、ラベルなしデータをそれに判別させた結果を用いて更に分類器を学習するという反復手法がうまく行くことが示されています。
一方で、各々のビューが潜在的な部分空間を共有しているという仮定の下、その部分空間を見つけて写像してあげようというアプローチもあります。その中でも、正準相関分析（CCA）を利用した手法がよく知られています。他にも、フィッシャー判別分析を利用した手法、Neighborhood Preserving Projectionsを利用した手法などがあります。
また、上述のマルチカーネル学習の手法を適用した論文も多数あります。というか、マルチビュー学習の問題に、類似尺度が複数考えられる場合に効果的なマルチカーネル学習の手法を適用するのは自然なアプローチと言えます。もう少し言うと、マルチカーネルを適用すること自体、一種のマルチビュー学習とみなせるんじゃないでしょうか*1。そのため、マルチカーネル学習を適用して解いた問題を他のマルチビュー学習手法で解けば、もっと良くなるってこともあるかも・・・
なお、この項目を書くにあたって、A Survey on Multi-view Learningというサーベイ論文を参考にしました。

マルチタスク学習

複数の関連するタスク同士で情報を共有することにより、予測精度を向上させようという考え方の枠組み。特徴ベクトルやラベルの定義域はタスク間で共通の場合が多いです。考え方は転移学習と類似しているので、その辺りの論文も参考になります。両者の違いは、マルチタスク学習の方は、タスク間で協調し合うことで全てのタスクの精度を向上させようという目的に対し、転移学習はある目標のタスクがあって、その目標タスクの精度向上が目的である点で違います。
マルチタスク学習の例は、書き手が異なる学習データを用いた手書き文字認識です。学習データに本人が書いたデータが少ない場合、他人の書いた文字データを利用することで精度を向上させようという試みは、マルチタスク学習とみなせます（タスク1はAさんの手書き文字認識器を学習、タスク2はBさんの手書き文字認識器を学習・・・）。それ以外にも、テキストデータを分類する問題で分野の違うテキストデータを用いたり、シーンの異なる画像データで物体検出を行ったりといった場面でマルチタスク学習が適用できます。
手法としては、異なるタスクのサンプルを学習する際、そのサンプルがどれだけ適合するかという重みを付けるアプローチや、正則化項などを用いてタスク間でパラメータが類似するように学習するアプローチなどがあります。例えば、前者には、タスク間の関係が共変量シフトという状況下であると仮定し確率密度比を重みづけした手法が、Covariate Shift Adaptation by Importance Weighted Cross Validation等の論文で提案されています*2。後者も行列ノルムの正則化を利用し、各タスクのパラメータがスパース、かつ、0になる要素がタスク間で共通（jointly sparse）になるように学習する手法Robust visual tracking via multi-task sparse learning等々様々な手法があります。また、上述のマルチビュー学習とマルチタスク学習の問題設定を組み合わせたA Graph-based Framework for Multi-Task Multi-View Learningなんてものもあります。
サーベイ論文としては、転移学習のサーベイですがA Survey on Transfer Learningが有名でしょう。検索すれば日本語の資料なんかもたくさん出てきます。

マルチラベル学習

その名の通り、1つのサンプルに複数のラベルが割り当てれる分類問題の枠組み。例えば、小説のカテゴリ分けを考えた場合、「SF＋ミステリー」、「ホラー＋恋愛＋ファンタジー」等のように、1つのジャンルに定められないことが多々あります。このような問題をマルチラベル学習では考えており、学習用サンプルに複数ラベル（カテゴリ）が割り当てられ、テスト用データの分類結果も複数ラベルが出力されるような手法が様々提案されています。
手法としては、大まかに分けると既存の学習手法をマルチラベルに拡張するアプローチと、マルチラベル問題をシングルラベルの問題に変換するアプローチの２つがあります。前者はSVMの拡張やAdaBoostの拡張、k近傍法の拡張など様々なものが提案されています。
後者の方もたくさんのアルゴリズムが提案されていています。例を挙げると、そのラベルが付いているか否かの2値分類問題に変換する手法（ジャンルがSFか否か、ミステリーか否かといった分類器をそれぞれ構築）、存在するラベルの組み合わせを1つのラベルとしてみなす手法（と言う1つのラベル、という1つのラベル、<ホラー＋恋愛＋ファンタジー>という1つのラベル）、複数ラベルが付いているサンプルは特徴ベクトル同一だがラベルが異なる別のサンプルとして学習する手法、などなど
マルチラベル学習の手法はほとんど知らなかったので、朱鷺の杜Wikiのマルチラベルの項目やここで紹介されていたチュートリアル論文Mining Multi-label Dataが非常に参考になりました。

マルチインスタンス学習

他の枠組みより多少問題設定がややこしいです。まず、用語として「bag」と呼ばれるものがあります。これは、複数のサンプル*3がひと塊りになったものです。マルチインスタンス学習では、学習データとして複数のbagと、それぞれのbagに正例か負例かというラベル情報が与えられます。ここで言う正例とは、bag内に1つでも正例に属するインスタンスがあること（負例インスタンスがいくつあってもよい）、対して負例はbag内全てのインスタンスが負例に属することを言います。bag内個々のインスタンスのラベルは分かりません。このような、学習データを与えられた上で、正例か負例かわからないbagを分類するというのがマルチインスタンス学習問題の枠組みです。
このような問題設定は、画像の分類問題に利用できます。例として、与えられた1枚の画像がビーチの画像か否かを判定するタスクを考えます。教師あり学習なので、学習用データとしてビーチの画像であると分かっている画像と、ビーチの画像ではないと分かっている画像が複数枚事前に与えられているとします。このような場合に、画像1枚の画像から領域ごとにわけて特徴ベクトルを抽出すると、1枚の画像=bag、抽出した各特徴ベクトル=インスタンスとみなせます。学習用データには画像（bag）に対して正例か負例かというラベル情報が割り当てられているので、このような問題設定は、まさにマルチインスタンス学習と言えます。それ以外にも、ある分子が薬として適切か否かを判定するタスク、文書データの分類などにも使われるようです。
アプローチとしては、与えられたbagセットからインスタンス単位での分類器を構築する、新規のbagに対しては、各インスタンスを分類器で判別してその結果を統合してbagがどちらに属するか判断するというアプローチ、bagレベルで学習するアプローチ（bag同士の距離を定めたりカーネル法を利用したり）、bagから1つの特徴ベクトルを構築（bag-of-keypointsなど）して学習するアプローチなどがあります。
まだちゃんと読んでいませんが、Multiple instance classification: Review, taxonomy and comparative studyというレビュー論文を参考にこの項目を書きました。

*1:例えば、カーネルの超パラメータを変えるだけでも、ある意味でビューを変えているというのが僕の考えです

*2:ここにあげた論文自体は、マルチタスク学習や転移学習を想定しているわけではありませんが、それらの一種とみなせると思います

*3:マルチインスタンス学習の枠組みでは「インスタンス」と呼びます。

動的計画法でテキストセグメンテーション

アルゴリズム 動的計画法

いやいや、お久しぶりです。
実に半年以上も更新が滞ってました。
ちゃんと生きてます。

以前、動的計画法によるシーケンシャルデータのセグメンテーションという記事を書きましたが、
今回はそれを応用して、テキストセグメンテーションを行おうと思います。

テキストセグメンテーションって？

テキストセグメンテーションとは、文章を意味的な段落に分割することを言います。
例えば、ある1つの文章が、野球の話題→サッカーの話題→尿道結石の話題→また野球の話題
といった具合に話題が推移したとします。
この文章を、これを下図のように話題ごとのセグメントに分割しようってのが、
テキストセグメンテーションです。

-----------------------------------------------------
今年もペナントレースが始まりました。
去年、セリーグは巨人、パリーグは楽天が優勝。
・・・
何はともあれ頑張ってもらいたいものです。
-----------------------------------------------------
そういえば、今年はサッカーワールドカップの年。
・・・
どの国が優勝するか注目です。
-----------------------------------------------------
ところで、私、尿道結石になってしまいました。
・・・
とても辛いので、早く治したいです。
-----------------------------------------------------
話を元に戻しますが、去年優勝した楽天は、
エースだった田中投手が抜け、
・・・
今から楽しみです。
-----------------------------------------------------

テキストセグメンテーションは様々な手法が提案され、
LDAのような潜在トピックモデルを利用した手法や接続詞に着目した手法などがあります。
今回は、セグメント内の単語のまとまり具合をスコアとし、
その合計スコアを動的計画法により最大化するという手法でセグメンテーションを行います。
探してみると以下の論文がやりたいことと一致しているので、こちらを参考にしました。
A. Kehagias, P. Fragkou, V. Petridis. 2003. Linear text segmentation using a dynamic programming algorithm. In Proceedings of the EACL, 171–178.

動的計画法でセグメンテーションを行う関数自体は以前の記事で作成したものを踏襲するので、
テキストデータクラスを定義し、そのメソッドとして、セグメント内のスコアを計算するh(t1,t2)と
長さを返すGetLength()さえ定義すれば実現できます。

前提

前提として、セグメントを考える際の最小単位は、センテンスとします。
つまり、シーケンスデータで言うところの1つのデータポイントは、
今回の場合は、1つのセンテンスとなります。
従って、文章に含まれるセンテンス数をTとすると、GetLengthメソッドは、
このTを返すことになります。
また、セグメントのスコアを計算する際、単語単位で計算しますが、
今回はストップワードを除いた名詞のみを抽出して利用します。
文章に含まれる名詞数（ストップワード除く）をLとします。

セグメントのスコア

まず、センテンス間の類似度を考えます。
この類似度は、センテンス間で共通の単語を含んでいれば、その個数分類似すると考えます。
そのため、次のようなL \times Tの行列{\bf F}を定義し（元論文と行と列が逆転してるので注意）、

F_{lt} = \begin{eqnarray} \left\{ \begin{array}{ll} 1 & \mbox{l-th word is in t-th sentence} \\ 0 & \mbox{else.} \end{array} \right. \end{eqnarray}

次の式で計算される行列{\bf D}はセンテンス間の類似度を表す行列となります。

{\bf D} = {\bf F}^T {\bf F} （ただし、D_{ii} = 0）

つまり、センテンスiとセンテンスjの類似度はD_{ij}となります。
なお、元論文はセンテンス間で共通の単語を複数含んでても、
その個数に限らずD_{ij} = 1となっています。

この行列{\bf D}を用いて、セグメント内のスコアを計算します。
今、セグメントの始点をt1、終点をt2とすると、
そのセグメントに含まれるセンテンスの番号tはt1 \leq t < t2になります
（t2番目のセンテンスはそのセグメントに含まれない）。
そしてそのスコア計算は以下のようになります。

J = \frac{\sum_{i=t1}^{t2-1} \sum_{j=t1}^{t2-1} D_{ij}}{(t2 - t1)}

要は、セグメントに含まれるセンテンスそれぞれの類似度を全部足し合わせ、
セグメントの長さで割ったものをスコアとしています。
元論文はセグメント長の平均や分散を予め用意した学習データから学習したり、
分母をr乗したりとか、もっと工夫していますが面倒なのでそこまでやりません。

ソースコード

単語の抽出

何はともあれ、単語を抽出しなければ始まりません。
MeCabを利用して形態素解析を行い、名詞を抽出します。
MeCabなど形態素解析は初挑戦で戸惑いましたが、
TaggerクラスのParseToNodeメソッドの戻り値であるnodeインスタンスの
featureを利用すれば品詞とその細分類が簡単に得られるとわかったので、これを利用します。
ストップワードとしては、品詞の細分類が代名詞、非自立、数、接尾のものとしています。
ということで、この操作を行い、抽出した単語のリストを返す関数を以下に示します。

import MeCab

def GetListOfWord(text):
    tagger = MeCab.Tagger('-Ochasen')
    encoded_text = text.encode('utf-8')
    node = tagger.parseToNode(encoded_text)
    WordList = []
    while node:
        if(node.feature.split(",")[0] == "名詞" and
            node.feature.split(",")[1] != "代名詞" and
            node.feature.split(",")[1] != "非自立" and
            node.feature.split(",")[1] != "数" and
            node.feature.split(",")[1] != "接尾"):
            if WordList.count(node.surface.decode('utf-8')) == 0:
                WordList.append(node.surface.decode('utf-8'))
        node = node.next
    return WordList

調べたところによると、半角記号が名詞のサ変接続になってしまうようですが、
とりあえず今回は気にしないことにします。

センテンスの抽出

特別なことはしません。
読点「。」、感嘆符「！」、疑問符「？」があれば、
そこがセンテンスの区切りとして抽出します。
センテンスのリストを返す関数は以下の通りです。

def GetSentenceList(text):
    sentence = u""
    SentenceList = []
    for t in text:
        sentence = sentence + t
        if t == u"。" or t == u"？" or t == u"！":
            SentenceList.append(sentence)
            sentence=u""
    return SentenceList
行列Fの計算

行列\bf{F}は疎な行列なので、
scipyのスパース行列クラスの１つlil_matrixとして扱います。
まず、上記の2関数を用いて、単語のリストとセンテンスのリストを取得します。
その後、forループを用い、各々センテンスで単語を抽出、
そのセンテンスに含まれる単語のインデックスを取得、
行列\bf{F}におけるそのインデックスの行に1を代入という流れになります。

import scipy.sparse as sp

def GetFMatrix(text):
    WordList = GetListOfWord(text)
    SentenceList = GetSentenceList(text)
    F = sp.lil_matrix((len(WordList) , len(SentenceList) ))
    i = 0
    for sentence in SentenceList:
        idx = [WordList.index(w) for w in GetListOfWord(sentence) if w in WordList]
        F[idx , i] = 1
        i += 1
    return (F , WordList , SentenceList)
テキストデータクラス

先述のとおり、テキストデータクラスのメソッドとして、
GetLength()とh(t1 , t2)さえ定義すれば、
前回作成した動的計画法プログラムがそのまま使えます。
まず、__init__メソッド内で、上述のGetFMatrix(text)関数を呼び、
行列\bf{F}と単語リスト、センテンスリストを得ます。
そして、行列\bf{D}を計算します。
後は、GetLength()はセンテンスの個数を返し、h(t1 , t2)は上述のスコア計算値を返せばOK。

import numpy as np

class TextSequence:
    def __init__(self , text):
        (F , self.WordList , self.SentenceList) = GetFMatrix(text)
        D = (F.T * F).toarray()
        D -= np.diag(np.diag(D))
        self.F = F
        self.D = D

    def GetLength(self):
        return len(self.SentenceList)

    def h(self , t1 , t2):
        return np.sum(np.sum(self.D[t1:t2][: , t1:t2])/(t2-t1))

実験

青空文庫から芥川龍之介の羅生門、鼻、トロッコ、黒衣聖母の4作品を引っ張ってきて、
これらを連結した文章を一つ作成しました（akutagawa.txtとしています）。
全部で543個のセンテンスで構成されています。
この文章に対し、セグメント数を4としてセグメンテーションを行い、
キチンと作品と作品の境にセグメントの区切りが来れば成功。
ということで、次のように実行しました。

>>> import DPSegmentation
>>> text =codecs.open('akutagawa.txt', 'r', 'utf-8').read()
>>> Seq=TextSequence(text)
>>> (T , score) = DPSegmentation.SequenceSegmentation(Seq , 4)

結果としては、綺麗に各作品セグメントで分割することに成功。
分割したテキストデータをここに貼るわけにもいかないので（貼ってもあまり意味ないし）、
各セグメントに対する単語のスコア（その単語がどれだけ、そのセグメントに関わるか）を
算出し、その上位5件を示すことにします。
このスコア算出のために以下に示すGetWordScore(t1 , t2)メソッドを
TextSequenceクラスに追加しました。

    def GetWordScore(self, t1 , t2):
        W=((self.F[: , t1:t2] * self.F[: , t1:t2].T).toarray()).sum(0)
        WordScoreList=[]
        for idx in W.argsort()[-1::-1]:
            WordScoreList.append((self.WordList[idx] , W[idx]))
        return WordScoreList

やってることは、ある単語がセンテンス間で共通する場合1、しない場合0となる操作を、
それがセグメント内の全センテンス間の組み合わせで求め、
その合計値を各々の単語のスコアとするというもの。

で、そのスコアを算出した結果は以下の通りです。

・セグメント1（羅生門セグメント）
(u'下人', 283.0),(u'老婆', 170.0),(u'死骸', 88.0),(u'雨', 79.0),(u'男', 70.0)

・セグメント2（鼻セグメント）
(u'供', 484.0),(u'鼻', 445.0),(u'弟子', 165.0),(u'僧', 165.0),(u'顔', 112.0)

・セグメント3（トロッコセグメント）
(u'良平', 240.0),(u'トロッコ', 200.0),(u'線路', 88.0),(u'土工', 77.0),(u'今', 46.0)

・セグメント4（黒衣聖母セグメント）
(u'利', 233.0),(u'麻', 233.0),(u'祖母', 218.0),(u'耶観', 191.0),(u'栄', 146.0)

それっぽい単語がキーワードとして抽出できているかなと思います。
やはり、タイトルの名詞や主人公の名前がセグメントに寄与してるようです。
また、黒衣聖母の「麻利耶観音」はどうやら形態素解析で、
「麻」と「利」と「耶観」と「音」にわけられてしまってますね・・・

所感

今回用いたセグメンテーション手法の特徴を挙げてみると、

・事前学習用データを一切使わない
・利用した単語は名詞のみ
・セグメント間の非類似度などは考慮していない
・スコア計算も至ってシンプル

と比較的お手軽な方法ですが、割と良い感じにセグメンテーション出来てるのかなと思います。
まぁ、別々の文章を無理やりくっつけて、さぁわけろと言ってるので、
うまい具合にできてもらわないと、困ると言えば困るのですが・・・

この手法の欠点を挙げると、kT^2オーダーの計算量なので、
長い文章だとかなり時間がかかります。
また、トピック分析をしてるわけでもないので、
各セグメントがどの話題について言及してるかは、見ることができません。
上記のように、セグメント内で共通する単語とかを見れば、
それなり話題が見えますが、やっぱり微妙かなと思います。

と言う訳で、pLSIやLDAなどの潜在トピックモデルを用いたセグメンテーション手法も
いずれは試してみたいと思います。

NIPS2012自棄読み　その2

機械学習

なんやかんやで前回から1ヶ月以上たってしまいました。
こりゃいかんということで、NIPS2012 part2いきます。
前回は機械学習の基礎的な手法を提案した論文を紹介しましたが、今回は応用よりの論文を2本取り上げてみました。

Learning Image Descriptors with the Boosting-Trick

画像のローカルな領域の特徴量（特徴ベクトル）記述手法としてSIFTやSURFなど勾配ベースの手法があります*1。これらは、回転・スケール変化などに不変であるため、物体検出等によく利用されますが、不変と言ってもやはり限界はつきもので、これらの特徴量を使っても検出に失敗することも多々あります。
そこで、ラベル付きのデータを用いて、この特徴量自体を教師あり学習的に求めて精度を上げようという枠組みがあります。ここで紹介する論文は、その枠組みの中で、領域の勾配情報を基にした弱学習器 + AdaBoostによって高次元に写像し、そこから特徴量に落とし込むという手法を提案しています。この弱学習器については、論文を参照していただくこととして、ここでは特徴量抽出までの大まかな流れを紹介します。
\bf{x}及び\bf{y}を画像のパッチとします。また、f(\bf{x} , \bf{y})を両画像パッチの類似度を示す関数とし、教師データのラベルl_iは与えられた両パッチ{\bf x}_i ,{\bf y}_iが類似しているか (l_i = 1)、類似していないか(l_i = -1)を表すとします。そして、以下の関数
\cal{L} = \exp \left(-l_i f({\bf x}_i , {\bf y}_i)\right)
を最小にするfは何ぞや、ということを考えます。
この論文の手法（Low-Dimensional Boosted Gradient Maps : L-BGM）は、M個の勾配ベースの弱学習器h_1({\bf x}) , h_2({\bf x}) , \cdots , h_M({\bf x})を用意し、それを用いてfを次のような関数として定義しています。

f({\bf x} , {\bf y}) = f_{LBGM}({\bf x} , {\bf y})=\sum_{i,j} \alpha_{ij}h_i({\bf x}) h_j({\bf y}) = {\bf h}^T({\bf x}) {\bf A} {\bf h}({\bf y})
これを上記の目的関数に当てはめて、弱学習器{\bf h}とその重み行列{\bf A}を学習すればよいのですが、両方を同時に最適化するのは難しい。そこで、最初にAdaBoostの手法を使って各{h}_iを学習し、その後に目的関数を最小にするような{\bf A}を確率的勾配法で求めるという2stepの手法を提案しています。また、AdaBoostにより各学習器の重み係数が出てきますが、その多くは0になるようです。そのため、M個の学習器の一部であるP個だけを使うことになります。ということで、実際にはM \times Mの行列{\bf A}ではなく、使用する学習器に対応する要素で構築したP \times Pの行列{\bf A}_Pを求め、使用することになります（当然、P < M）。
これらを求めると次は、画像パッチ{\bf x}に対する各弱学習器の分類結果及び{\bf A}_Pから特徴量を構築します。{\bf A}_Pをランクdの対称行列とすると、この行列は次のように分解できます。
{\bf A}_P = {\bf B}{\bf W}{\bf B}^T = \sum_{k=1}^d w_k {\bf b}_k {\bf b}_k^T

{\bf B}はP \times dの行列、{\bf W}は要素が1 or -1の対角行列です。すると、類似度は
f_{LBGM}({\bf x} , {\bf y}) = \left({\bf B}^T {\bf h} ( {\bf x}) \right)^T {\bf W} \left({\bf B}^T {\bf h} ( {\bf y}) \right)

というように、特徴ベクトル同士の符号付き内積の形で表せます。
従って、{\bf B}^T {\bf h}( {\bf x})を画像パッチ{\bf x}の特徴量として用いましょうというのがこの論文で提案している手法です。{\bf A}_pに正定値行列の制約をつければ、{\bf W}消せるしもっと良くなるんじゃないの？と考えましたが、論文によるとそのような制約を入れても結果はさほど変わらないようです。
いったん高次元に写像してから特徴ベクトルに落とし込むという点で、カーネルトリックに類似していますが、この論文で主張されているように、カーネルトリックに比べて直感的で良いかなと思います。また、弱学習器次第では他の分野にも応用できそうです。

A P300 BCI for the Masses: Prior Information Enables Instant Unsupervised Spelling

自分が学生時代やっていた研究とドンピシャだったので、思わず目を引いた論文。
脳波など脳活動を計測し、その信号でコンピュータを操作する技術をブレインコンピュータインタフェース（BCI）と言います。その1つとしてP300 Spellerと呼ばれる脳波を利用して文字入力を行うシステムが提唱されています。P300 Spellerでは、ユーザに行または列が１つずつパッと光るアルファベットや記号が羅列されたマトリックスを提示し、ユーザはその中の入力したい文字に着目します。すると、ユーザが入力したい文字が光った時、頭頂部付近から計測される脳波には、光ってから約300ms後にピーク成分が発生します。従って、計測した脳波にP300成分が発生しているか否かを識別すれば、ユーザがどの文字を入力したいか分かるという仕組みです。
ということで、この識別に機械学習の手法を利用した研究が過去多く提案されているのですが、その多くは教師あり学習の手法を利用するために、ユーザ本人から教師データとなる脳波データを予め採取する必要があります。この教師データ取得のための事前計測は、ユーザにとって負担になるため、P300 Spellerの1つの課題として挙がっています。じゃあ、他人の脳波予め採取しておいてそれを識別器に学習させればいいじゃないかと考えますが、個人差の関係上単純な方法ではうまくいきません。
そんな中、教師なし学習によるP300 Spellerの手法を今回取り上げる論文の著者らが提案しました*2。これは、ユーザが入力したい文字はどれであるかという変数を隠れ変数とみなし、EMアルゴリズムによってパラメータを学習する手法です。しかし、この手法ではユーザが使用していくうちに徐々に精度が良くなっていきますが、最初はほぼランダムな結果になってしまうという問題があります。そこで著者らは、他人の脳波データ及び言語モデル（n-gramモデル）を事前情報として利用することで、この問題を解消しようという手法を今回取り上げる論文で提案しています。
他人の脳波を利用する手法は転移学習の枠組みの1つといえます。この手法のモデルとして、著者らは以下のようなものを提案しています。

p({\bf \mu}_w) = {\cal N} (0 , \alpha_p {\bf I} ) , p({\bf w}_s |{\bf \mu}_w) = {\cal N} ({\bf w}_s |{\bf \mu}_w , \alpha_s {\bf I})
p(c_{s,t}) = \frac{1}{C} , p({\bf x}_{s , t , i} | c_{s,t} , {\bf w}_s , \beta_s) = {\cal N} ({\bf x}_{s , t , i} {\bf w}_s | y_{s , t , i}(c_{s,t}) , \beta_s)

ここで、c_{s , t}は被験者sがt文字目として入力したい文字を表します。また、{\bf x}_{s , t , i}は脳波の特徴ベクトル（行ベクトル）で、被験者s、t文字目においてi回目に取得した脳波データという意味になります。y_{s , t , i}(c_{s,t})は、そのi回目において、入力した文字c_{s , t}が光ったか否かを表すラベル変数です。更に、{\bf w}_sは今求めたい被験者sのパラメータベクトル（脳波の特徴ベクトルをP300が発生しているか否かに写像するベクトル）、\alpha_s及び\beta_sはそれぞれのガウス分布の精度でこれもEMアルゴリズムで逐次更新されます。{\bf \mu}_wは各被験者のパラメータベクトルの事前情報となる平均ベクトル、\alpha_pはそのガウス分布の精度です。
このようなモデルを元にc_{s , t}を隠れ変数と考えてEMアルゴリズムを適用して、{\bf w}_s, \alpha_s, \beta_sを求めます。大まかな学習の流れは以下の通り。

まず、最初のS人の被験者についての学習は、{\bf \mu}_w = {\bf 0} ,\alpha_p = 0として次の操作を収束するまで繰り返します。

  1. c_{s , t}の事後分布p(c_{s , t} | {\bf X}_{s ,t} , {\bf w}_s^{old} , \beta_s^{old})を求める
  2. \arg \displaystyle{\max_{{\bf w}_s ,\beta_s}} \sum_{c_s} p(c_{s} | {\bf X}_{s} , {\bf w}_s^{old} , \beta_s^{old}) \ln p( {\bf X}_{s} , c_s | {\bf w}_s , \beta_s) + \ln p({\bf w}_s |{\bf \mu}_w)を求めて更新する
  3. \alpha_sを更新する

次に被験者s=1 \cdots Sの{\bf w}_s及び\alpha_sを得たうえでの、S + 1人目の被験者についての学習ですが、これは

{\bf \mu}_w = \frac{1}{\alpha_p} \sum_{s = 1}^{S} \alpha_s {\bf w}_s , \alpha_p = \sum_{s = 1}^{S} \alpha_s
として上記同様の繰り返しを実行します。

具体的な更新式は書くのがしんどくなってきたので割愛します。論文に載っているのでそちらを参照してください。
上記のモデルで文字の出現確率p(c_{s,t})は\frac{1}{C}と一定値として扱っていましたが、これをn-gramの言語モデルに置き換えた手法もこの論文中で提案しています。
手法自体は素直にEMアルゴリズムを適用してるので、結構理解しやすかったです。後、P300 Spellerの学習アルゴリズムに関する論文が今でも結構出てるのが少し意外な気がしました。この論文に僕が知らなかったデータセットも載ってたので、久々になんかやってみようかなと思ったり思わなかったり・・・

ということで、前回のと合わせて、あきらめずに読んだNIPS 2012の論文うち4本紹介しました（紹介というより備忘録に近いですが）。ときどき、こうやって学生の頃やってた研究を思い出すのもいい刺激になります。

*1:通常、SIFTやSURFと言うと特徴点検出と特徴量の2段階をひっくるめた手法を指しますが、ここでは特徴量だけのことを指します

*2:P.-J. Kindermans, D. Verstraeten, and B. Schrauwen, A bayesian model for exploiting application constraints to enable unsupervised training of a P300-based BCI. PLoS ONE, 04 2012