ハイパーデクノボウブログ

2014-09-26

【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{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})

(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波をそれぞれ適当なゲインをかけて足し合わせたものを使います。で、結果はこちら。
f:id:YamagenSakam:20140926120453p:image:w640

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

まとめ

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

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

2014-01-19

マルチ○○学習まとめ

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

マルチカーネル学習

複数の特徴ベクトル表現、または、類似の尺度が考えられるタスクに対し、複数のカーネル関数を組み合わせたカーネル学習を考えましょうという手法の枠組みです。例えば、画像データの場合、特徴ベクトルとしてbag-of-keypointsを構築したり、色情報を利用したりと、様々な類似尺度を考えることができます。そこで、それぞれの類似尺度に対応するカーネル関数を構築して組み合わせれば、ハッピーになれるだろうという考え方です。
多くの場合、次のようにカーネル関数の線形和を考えます。
k(¥bf{x} , ¥bf{y}) = ¥sum_{i=1}^{N} ¥beta_i k_i(¥bf{x_i , y_i})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つのジャンルに定められないことが多々あります。このような問題をマルチラベル学習では考えており、学習用サンプルに複数ラベル(カテゴリ)が割り当てられ、テスト用データの分類結果も複数ラベルが出力されるような手法が様々提案されています。
手法としては、大まかに分けると既存の学習手法をマルチラベルに拡張するアプローチと、マルチラベル問題をシングルラベルの問題に変換するアプローチの2つがあります。前者はSVMの拡張やAdaBoostの拡張、k近傍法の拡張など様々なものが提案されています。
後者の方もたくさんのアルゴリズムが提案されていています。例を挙げると、そのラベルが付いているか否かの2値分類問題に変換する手法(ジャンルがSFか否か、ミステリーか否かといった分類器をそれぞれ構築)、存在するラベルの組み合わせを1つのラベルとしてみなす手法(<SF>と言う1つのラベル、<SF+ミステリー>という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:マルチインスタンス学習の枠組みでは「インスタンス」と呼びます。

2013-12-30

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

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

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

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

テキストセグメンテーションとは、文章を意味的な段落に分割することを言います。
例えば、ある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 ¥time 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}  ¥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とすると、
そのセグメントに含まれるセンテンスの番号tt1 ¥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のスパース行列クラスの1つ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などの潜在トピックモデルを用いたセグメンテーション手法
いずれは試してみたいと思います。

2013-05-05

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 ¥(-l_i f(¥bf{x}_i , ¥bf{y}_i)¥)
を最小にする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}) = ¥(¥bf{B}^T ¥bf{h} ( ¥bf{x}) ¥)^T ¥bf{W} ¥(¥bf{B}^T ¥bf{h} ( ¥bf{y}) ¥)

というように、特徴ベクトル同士の符号付き内積の形で表せます。
従って、¥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では、ユーザに行または列が1つずつパッと光るアルファベットや記号が羅列されたマトリックスを提示し、ユーザはその中の入力したい文字に着目します。すると、ユーザが入力したい文字が光った時、頭頂部付近から計測される脳波には、光ってから約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}は被験者st文字目として入力したい文字を表します。また、¥bf{x}_{s , t , i}脳波の特徴ベクトル(行ベクトル)で、被験者st文字目において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¥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

2013-03-18

NIPS2012自棄読み その1

気が付いたら3月も中旬です。昨年末に27歳の誕生日を迎えたばかりだと思っていたら、もう3ヶ月経ってしまいました。このまま、あっという間に三十路を迎えるのでしょう。時の流れは早くて残酷です。生え際の後退もかなり進行してきました。悲しいです。昔の写真をみると涙が出そうになります。悲しすぎるので、僕の誕生日より少し前に開催されたNIP2012の論文を自棄読みしました。

Multi-Task Averaging

多くのマルチタスク学習は分類問題や回帰問題に対して提案された手法ですが、この論文は母平均をマルチタスク学習で推定するという珍しい(?)手法Multi-Task Averaging(MTA)を提案しています。異なるタスク同士でも、お互い情報を共有し合って母平均を推定することで推定精度を上げられると主張しています。早速、目的関数です。

¥{ {Y}_t^* ¥}_{t=1}^T = {¥rm{argmin}} ¥limits_{¥{ ¥hat{Y}_t ¥}_{t=1}^T} ¥frac{1}{T} ¥sum_{t=1}^{T} ¥sum_{i=1}^{N_t} ¥frac{¥(Y_{ti} - ¥hat{Y}_t ¥)^2}{¥sigma_t^2} + ¥frac{¥gamma}{T^2} ¥sum_{r=1}^{T} ¥sum_{s=1}^{T} A_{rs}(¥hat{Y}_r - ¥hat{Y}_s)

Y_{ti}タスクti番目のサンプル、A_{rs}タスクrタスクsの類似度を表します。その他記号のNotationは論文の方をご覧ください。この式の第1項は各タスクの経験誤差最小化の作用、第2項にはタスク間の類似度が高いほどその推定平均値の差が小さくなるように制御する作用があります。
ここで、A_{rs}が全て0以上、かつ、A_{rs}を要素とする行列¥bf{A}が対称行列だとすると、上記の最小化問題は解析的に解を得ることができます。

¥bf{Y}^* = ¥(¥bf{I} + ¥frac{¥gamma}{T} ¥Sigma ¥bf{L}¥)^{-1} ¥bf{¥bar{Y}}

¥Sigma¥Sigma_{tt}={¥sigma_t^2}/{N_t}という要素を持つ対角行列、¥bf{L}¥bf{A}のグラフラプラシアン¥bf{¥bar{Y}}は各タスクのサンプル平均¥(¥frac{1}{N_t}¥sum_{i=1}^{N_t}Y_{it}¥)を縦に並べたベクトルです。
さて、この論文ではタスク数が2の場合について解析し、このタスク間の実際の母平均の差¥Deltaがそれぞれのタスク分散に比べて小さい時、MTAによる推定結果の平均二乗誤差は、単純なサンプル平均による推定結果の平均二乗誤差より小さくなるという結論を得ています(つまり、タスク間の母平均が近いならMTAは有効)。また、2つのタスクの平均二乗誤差の和を最小にするaは、

a^* = ¥frac{2}{¥Delta^2}

として得られます。
この結果を用いて、タスク間の類似度行列¥bf{A}を求める方法を2つ提案しています。ここでは、そのうちの1つであるConstant MTAについて少しだけ記述します。やりたいことは、以下の式のように、実際の各タスクの母平均¥bf{¥mu}MTAによる推定結果の平均二乗誤差を最小化することです。

R¥(¥bf{¥mu} , ¥bf{W}¥bf{¥bar{Y}} ¥) = E ¥[¥(¥bf{W} ¥bf{¥bar{Y}} - ¥bf{¥mu} ¥)^T ¥(¥bf{W} ¥bf{¥bar{Y}} - ¥bf{¥mu} ¥)¥] = ¥rm{tr} ¥(¥bf{W} ¥Sigma ¥bf{W} ¥) + ¥bf{¥mu}^T ¥(¥bf{I} - ¥bf{W} ¥)^T ¥(¥bf{I} - ¥bf{W} ¥) ¥bf{¥mu}
ただし、
 ¥bf{W}= ¥(¥bf{I} + ¥frac{¥gamma}{T} ¥Sigma ¥bf{L}¥)^{-1}

これを最小化する¥bf{A}を求めたいのですが、タスク数が2より大きい場合、求めたいパラメータの数が多くなり解析的に¥bf{A}を求めるのは、格段に難しくなります。そこで、行列の全要素はaであるという大胆な仮定を置きます。更に、¥gamma=1で全タスク分散は同じであるという仮定を置くと、損失を最小にするaは、

a^* = ¥frac{2}{ ¥frac{1}{T ¥(T - 1 ¥)} ¥sum_{r=1}^{T} ¥sum_{s=1}^{T} ¥(¥mu_r - ¥mu_s ¥)^2}

となります。ただし、¥muは実際の母平均の値で当然分からないので、結局のところサンプル平均の¥bar{y}を用います。とどのつまり、タスク間の平均差の二乗¥Delta^2タスクの組み合わせ分求めて平均したものを分母に用いるという方法で、タスク数が2の場合の拡張をしてることになります。結果をみると、やはりタスク間の類似度が高ければ、このConstant MTAは良く働くようです。

Hamming Distance Metric Learning

入力ベクトル¥{-1 , 1¥}の2値のベクトルバイナリコードベクトル)として表現することで、近傍探索などにおいてメモリ使用量や計算速度の改善が期待できます。では、どのようにバイナリコードベクトルに変換するのが適切か?これをサンプル集合から学習するのがこの論文の目的です。
¥bf{x} ¥in ¥bf{R}^pを入力ベクトル¥bf{h} ¥in ¥bf{R}^qを目的のバイナリコードベクトルとすると、その関係を

¥bf{h} = b¥(¥bf{x};¥bf{w}¥) = ¥rm{sign}¥(f¥(¥bf{x};¥bf{w}¥)¥)

という式で表すことができます。fは、¥bf{R}^p ¥rightarrow ¥bf{R}^qとする線形または非線形写像です。多くのバイナリコード変換の学習は、fを決定づけるパラメータ¥bf{w}を学習します。これまで、このfとして、線形写像を使うよとか多段ニューラルネット使いましょうなどの、様々な写像が提案されていますが、この論文ではfを決めず、それらを包括するような枠組みを提案しています。
まず、この論文の前提としてサンプル集合¥cal{D}は、¥bf{x}¥bf{x}^+の類似度は¥bf{x}¥bf{x}^-の類似度より高いという三つ組¥(¥bf{x} , ¥bf{x}^+ , ¥bf{x}^-¥)の集合として得られるものとします(入力インプットベクトル同士の類似度を計算できるなら、この三つ組は得ることができます)。そして、バイナリコードへの変換を学習するために以下のような損失関数を定義しています。

¥cal{L}¥( ¥bf{w} ¥) = ¥sum_{¥(¥bf{x} , ¥bf{x}^+ , ¥bf{x}^-¥) ¥in ¥cal{D}} ¥ell_{¥rm{triple}} ¥(b¥(¥bf{x};¥bf{w}¥) , b¥(¥bf{x}^+;¥bf{w}¥) , b¥(¥bf{x}^-;¥bf{w}¥) ¥) + ¥frac{¥lambda}{2} ¥parallel ¥bf{w} ¥parallel_2^2
ただし、
¥ell_{¥rm{triple}} ¥(¥bf{h} , ¥bf{h}^+ , ¥bf{h}^-¥) = ¥[ ¥parallel ¥bf{h} - ¥bf{h}^+ ¥parallel_H - ¥parallel ¥bf{h} - ¥bf{h}^- ¥parallel_H + 1 ¥]_+

¥[¥alpha ¥]_+ = ¥max ¥(¥alpha , 0 ¥)¥parallel ¥cdot ¥parallel_Hベクトルの非ゼロ要素の数。
この¥cal{L}¥( ¥bf{w} ¥)を最小化する¥bf{w}を求めることで、類似度の順序関係を保持したバイナリコード変換が得られることを期待できます。しかし、この目的関数は非凸かつ¥bf{w}微分不可という面倒な性質を持っています。ということで、損失の上限式に対して最小化を行うことを考えます。¥ell_{¥rm{triple}} ¥(b¥(¥bf{x};¥bf{w}¥) , b¥(¥bf{x}^+;¥bf{w}¥) , b¥(¥bf{x}^-;¥bf{w}¥) ¥) の上限は以下のようになります。

¥max_{¥bf{g} , ¥bf{g}^+ , ¥bf{g}^-} ¥{¥ell_{¥rm{triple}} ¥(¥bf{g} , ¥bf{g}^+ , ¥bf{g}^-¥) + ¥bf{g}^T f¥(¥bf{x};¥bf{w}¥) +  ¥bf{g}^+^T f¥(¥bf{x}^+;¥bf{w}¥) +  ¥bf{g}^-^T f¥(¥bf{x}^-;¥bf{w}¥)¥} ¥¥  ¥hspace{20cm}  - ¥max_{¥bf{h}} ¥hspace{2cm} ¥bf{h}^T f¥(¥bf{x};¥bf{w}¥) - ¥max_{¥bf{h}^+} ¥hspace{2cm} ¥bf{h}^+^T f¥(¥bf{x};¥bf{w}¥) - ¥max_{¥bf{h}^-} ¥hspace{2cm} ¥bf{h}^-^T f¥(¥bf{x};¥bf{w}¥)

こうすることで¥bf{w}微分可能になります。後は確率的勾配法を用いてこれを解いていきます。具体的には、¥bf{w}に初期値を与えた後、サンプル集合¥cal{D}から三つ組¥(¥bf{x} , ¥bf{x}^+ , ¥bf{x}^-¥)をランダムに1つ抽出し、それに対し上の式を満たす¥bf{g} , ¥bf{g}^+ , ¥bf{g}^- , ¥bf{h} , ¥bf{h}^+ , ¥bf{h}^-を求めてから、次のように¥bf{w}を更新します。

¥bf{w}^{t+1} ¥leftarrow ¥bf{w}^{t} + ¥eta ¥[¥frac{¥partial f¥(¥bf{x};¥bf{w}¥)}{¥partial ¥bf{w}} ¥(¥bf{h} - ¥bf{g} ¥) + ¥frac{¥partial f¥(¥bf{x}^+;¥bf{w}¥)}{¥partial ¥bf{w}} ¥(¥bf{h}^+ - ¥bf{g}^+ ¥) + ¥frac{¥partial f¥(¥bf{x}^-;¥bf{w}¥)}{¥partial ¥bf{w}} ¥(¥bf{h}^- - ¥bf{g}^- ¥)  - ¥lambda ¥bf{w} ¥]

この手順を繰り返すして、パラメータ¥bf{w}を推定します。
ところが、この式にある¥bf{g} , ¥bf{g}^+ , ¥bf{g}^-を求めるのにも一工夫が必要です。この¥bf{g} , ¥bf{g}^+ , ¥bf{g}^-を求める過程も、このベクトル-11のどちらかしかとらないことを利用し動的計画法で解くなど、なかなかテクニカルで面白いですが、この部分の詳細は論文を参照してください。
あと、論文にある式(8)は¥ell ’ ¥(¥alpha ¥) = ¥[ ¥alpha - 1 ¥]_+となっていますが、¥ell ’ ¥(¥alpha ¥) = ¥[ ¥alpha + 1 ¥]_+のような気がします。どうなんでしょう・・・

とりあえず、今回はここまでにします。近いうちに第2弾を書ければなぁと思います。

Connection: close