NIPS2017論文メモ

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

Deep Subspace Clustering Networks

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

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

これらの方法では非線形なクラスタリングは行えない。そこで、オートエンコーダを利用しようというのがこの論文の提案。オートエンコーダはニューラルネットを入力と出力を同じデータとして学習させることで、中間層の出力が非線形な特徴抽出となるという手法。入力から特徴を抽出する中間層まで層をエンコーダ層、特徴を抽出する中間層から出力までの層をデコーダ層と呼ぶ。この論文では、このエンコーダ層とデコーダ層の間にとするSelf-Expressive層を新たに設けることを提案(はエンコーダ層で特徴抽出されたサンプルの集合)。これは、線形の全結合となる。このと他の層の重みを求めるため、下記の損失関数を最小化。

はニューラルネットの重み係数集合、はエンコーダ層の出力集合（バッチ）を表す。ノルム項はとの両方を実験。

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

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

Mixture-Rank Matrix Approximation for Collaborative Filtering

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

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

本論文では、このPMFをベースに複数のランクkを重み付きで混合させたモデルを提案。具体的には下記の確率モデルを構築している。

ここで、はランクkの行列のi行目（ユーザiの特徴ベクトル）、はランクkの行列のj行目（アイテムjの特徴ベクトル）を表す。また、はユーザiにおけるランクkモデルの重み、はアイテムjにおけるランクkモデルの重み。これにより、ランクkのモデルがユーザi、アイテムjにマッチするならは大きくなるはずであり、マッチしないならは小さくなるはず（期待的には0になる）。これらのをMAP推定により求めようというのが提案手法。

MAP推定のため事前分布を与える。、の事前分布はそれぞれ、平均、共分散、のガウス分布としている（行列の事前分布はこれらの総乗）。また、はそれぞれスパースなベクトルにしたいので、の事前分布はラプラス分布とする（それぞれ、平均、分散というパラメータを持つ）。これらを積算することによって導き出される事後分布のlogにマイナスをつけた損失関数を最小化する。実際には、そのままの損失関数を最小化するのは難しいので、イェンセンの不等式で導き出される下限を最小化する。

はを固定して確率的勾配法で解く。一方、　も学習したいので、今度はを固定して最尤推定する（こちらは、解析的に求まる）。これらの操作を収束するまで繰り返す。

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は次の式で、表されます。

(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を構築したり、色情報を利用したりと、様々な類似尺度を考えることができます。そこで、それぞれの類似尺度に対応するカーネル関数を構築して組み合わせれば、ハッピーになれるだろうという考え方です。
多くの場合、次のようにカーネル関数の線形和を考えます。
（はカーネル関数の数）
このカーネル関数を学習の目的関数に組み込んで、最適なを求めるのが、メジャーな手法といえるでしょう。必要なカーネル関数の重み（）だけ大きくなるので、どれが学習に寄与する類似尺度かもわかるというメリットもあります。
有名なのは、大規模データをマルチカーネル+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つのセンテンスとなります。
従って、文章に含まれるセンテンス数をとすると、GetLengthメソッドは、
このを返すことになります。
また、セグメントのスコアを計算する際、単語単位で計算しますが、
今回はストップワードを除いた名詞のみを抽出して利用します。
文章に含まれる名詞数（ストップワード除く）をとします。

セグメントのスコア

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

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

（ただし、）

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

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

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

ソースコード

単語の抽出

何はともあれ、単語を抽出しなければ始まりません。
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

行列の計算

行列は疎な行列なので、
scipyのスパース行列クラスの１つlil_matrixとして扱います。
まず、上記の2関数を用いて、単語のリストとセンテンスのリストを取得します。
その後、forループを用い、各々センテンスで単語を抽出、
そのセンテンスに含まれる単語のインデックスを取得、
行列におけるそのインデックスの行に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)関数を呼び、
行列と単語リスト、センテンスリストを得ます。
そして、行列を計算します。
後は、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)

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

所感

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

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

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

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

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

NIPS2012自棄読み その2

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

Learning Image Descriptors with the Boosting-Trick

画像のローカルな領域の特徴量（特徴ベクトル）記述手法としてSIFTやSURFなど勾配ベースの手法があります*1。これらは、回転・スケール変化などに不変であるため、物体検出等によく利用されますが、不変と言ってもやはり限界はつきもので、これらの特徴量を使っても検出に失敗することも多々あります。
そこで、ラベル付きのデータを用いて、この特徴量自体を教師あり学習的に求めて精度を上げようという枠組みがあります。ここで紹介する論文は、その枠組みの中で、領域の勾配情報を基にした弱学習器 + AdaBoostによって高次元に写像し、そこから特徴量に落とし込むという手法を提案しています。この弱学習器については、論文を参照していただくこととして、ここでは特徴量抽出までの大まかな流れを紹介します。
及びを画像のパッチとします。また、を両画像パッチの類似度を示す関数とし、教師データのラベルは与えられた両パッチが類似しているか ()、類似していないか()を表すとします。そして、以下の関数

を最小にするは何ぞや、ということを考えます。
この論文の手法（Low-Dimensional Boosted Gradient Maps : L-BGM）は、個の勾配ベースの弱学習器を用意し、それを用いてを次のような関数として定義しています。

これを上記の目的関数に当てはめて、弱学習器とその重み行列を学習すればよいのですが、両方を同時に最適化するのは難しい。そこで、最初にAdaBoostの手法を使って各を学習し、その後に目的関数を最小にするようなを確率的勾配法で求めるという2stepの手法を提案しています。また、AdaBoostにより各学習器の重み係数が出てきますが、その多くは0になるようです。そのため、個の学習器の一部である個だけを使うことになります。ということで、実際にはの行列ではなく、使用する学習器に対応する要素で構築したの行列を求め、使用することになります（当然、）。
これらを求めると次は、画像パッチに対する各弱学習器の分類結果及びから特徴量を構築します。をランクの対称行列とすると、この行列は次のように分解できます。

はの行列、は要素が1 or -1の対角行列です。すると、類似度は

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

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つといえます。この手法のモデルとして、著者らは以下のようなものを提案しています。

,
,

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

まず、最初の人の被験者についての学習は、として次の操作を収束するまで繰り返します。

1. の事後分布を求める
2. を求めて更新する
3. を更新する

次に被験者の及びを得たうえでの、人目の被験者についての学習ですが、これは

,
として上記同様の繰り返しを実行します。

具体的な更新式は書くのがしんどくなってきたので割愛します。論文に載っているのでそちらを参照してください。
上記のモデルで文字の出現確率はと一定値として扱っていましたが、これを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

NIPS2012自棄読み その１

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

Multi-Task Averaging

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

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

はという要素を持つ対角行列、はのグラフラプラシアン、は各タスクのサンプル平均を縦に並べたベクトルです。
さて、この論文ではタスク数が2の場合について解析し、このタスク間の実際の母平均の差がそれぞれのタスクの分散に比べて小さい時、MTAによる推定結果の平均二乗誤差は、単純なサンプル平均による推定結果の平均二乗誤差より小さくなるという結論を得ています（つまり、タスク間の母平均が近いならMTAは有効）。また、2つのタスクの平均二乗誤差の和を最小にするは、

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

ただし、

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

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

Hamming Distance Metric Learning

入力ベクトルをの２値のベクトル（バイナリコードベクトル）として表現することで、近傍探索などにおいてメモリ使用量や計算速度の改善が期待できます。では、どのようにバイナリコードベクトルに変換するのが適切か？これをサンプル集合から学習するのがこの論文の目的です。
を入力ベクトル、を目的のバイナリコードベクトルとすると、その関係を

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

ただし、

、はベクトルの非ゼロ要素の数。
このを最小化するを求めることで、類似度の順序関係を保持したバイナリコード変換が得られることを期待できます。しかし、この目的関数は非凸かつで微分不可という面倒な性質を持っています。ということで、損失の上限式に対して最小化を行うことを考えます。の上限は以下のようになります。

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

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

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

【Python】動的計画法によるシーケンシャルデータのセグメンテーション

学生の頃、表題のアルゴリズムを研究に応用しようとしていましたが、すっかり忘れているので、思い出すために改めてPythonで実装してみました。当時は、matlabでfor文を回して実現していましたが、今回は、勉強がてらラムダ式やら再帰呼び出しやらを駆使して書いてみようと思います。

シーケンシャルデータのセグメンテーションとは

まず、今回取り上げる問題はどういうものか、等間隔でサンプリングされた離散的な時系列データを例にとって説明します。下図のように、時系列データを個のセグメント分割すると、それぞれのセグメントごとに何らかのスコアが得られるとします。さて、セグメントの境界をどこに置けば、スコアの合計値を最大化できるのでしょう？というのが今回解きたい問題です。

具体的な応用を1つ挙げると時系列データの近似があります。例えば、下図のようにセグメント内の平均値で時系列データを近似し、その二乗誤差平均が最小になるようなセグメンテーションを求めるというものがあります。

これを、上記のスコア最大化という文脈に当てはめると以下のような評価関数になります。

ここでは評価関数の最大化という問題に当てはめたいので、誤差の二乗平均にマイナス付与したものをスコアとしています。ここでの時系列データは離散的にサンプリングされたデータであり、そのサンプル数は有限個（個）とするので、はである整数です。このを最大化するの組み合わせを求めるために、今回取り上げる動的計画法（Dynamic Programming, 以下DP）がよく利用されます。

DPをどのように適用するか

DPは1950年代にBellmanによって考案されたアルゴリズムで、DPを応用した関数のセグメンテーション近似手法も同著者によって提案されています*1。DPの基本的な事項は下記の書籍が参考になります。僕もDPの元論文ではなく、この本を参考にセグメンテーションを実装しています。

これなら分かる最適化数学―基礎原理から計算手法まで

さて、各セグメントのスコアはセグメントを囲む境界位置の関数になります。そのため、とすると解きたい問題は以下のようになります。

これは次のような再帰的な部分問題に分割することができます。

ここで、を求めることが考えてみると、各の値に対して、を最大にするを探索しなければなりません。そのため、毎回の計算結果が必要になります。その際、を予め求めその結果を保持しておけば、計算時間の節約になります。同様にを求めるときも、予めの結果を保持しておけば、時間的に効率が良いです。このように途中までの計算結果を保持しておき、次のステップで利用するというのはまさに動的計画法の考え方です。
この操作を繰り返していくと、最終的には目的関数の最大値とその時のが求まります。が分かれば、を最大にするは既に計算して保持してあるはずなので、すぐにわかります。同様に遡っていくことで、最適なまで求めることができます。これが動的計画法を利用したシーケンスデータのセグメンテーションアルゴリズムです。

動的計画法を実現する再帰関数

ここから、書いたソースコードを紹介していきます。なお、いつものように

import numpy as np

としています。では、まずセグメンテーションの起点となる関数です。

def SequenceSegmentation(SequenceData , k):
    f_ini = lambda t : (SequenceData.h(0 , t) , 0)
    ResultList = []
    DPSegmentation(SequenceData , k , f_ini , 1 , ResultList)
    ResultList.reverse()
    T = np.zeros(k + 1)
    T[1:k] = np.array(ResultList)[: , 1]
    T[k] = SequenceData.GetLength()
    return T

シーケンシャルデータのオブジェクト（詳しくは後述）とセグメント数を与え、最適なセグメント境界位置を返します。実際に動的計画法でセグメンテーションを行う部分は、この関数で呼びだしているDPSegmentationという再帰関数です。f_iniは上述のにあたる関数オブジェクトです。ただし、上記の数式上ではその値を返す関数ですが、プログラムでは値とその時に選ばれたのタプル（つまり、とのタプル）を返す関数としています。f_iniも同様で、なのでタプルの第2要素を全て0とします。これを初期値としてDPSegmentationに与えると、ResultListに各セグメントの境界位置（実際には、それとの値のタプル）がリストで得られるという寸法です。
ということで、その動的計画法の再帰関数です。

def DPSegmentation(SequenceData , k , f_old , i , ResultList):
    if i == k - 1:
        N = SequenceData.GetLength()
        Result = GetMax((lambda t_old: f_old(t_old)[0] +
                            SequenceData.h(t_old , N )) , range(i , N ) )   #(4)
        ResultList.append(Result)
        return Result
    g = lambda t: lambda t_old: f_old(t_old)[0] + SequenceData.h(t_old , t) 
    f = lambda t: GetMax(g(t) , range(i , t))                               #(1)
    f = MemoizeFunc(f)                                                      #(2)
    Result = DPSegmentation(SequenceData , k , f , i + 1 , ResultList)      #(3)
    ResultList.append(f(Result[1]))
    return f(Result[1])

第1引数はシーケンシャルデータクラスから生成したオブジェクトです。見て分かる通りシーケンスの長さ（）を返すメソッドGetLength()と2つの引数（セグメントの境界位置）を受け取りセグメントのスコアを返すh(t1 , t2)が定義されている必要があります。第2引数はセグメント数、第3引数のf_oldは前ステップで作成した上述の関数にあたる関数オブジェクトを入れます。第4引数iはそのf_oldの添え字番号を表し、第5引数のリストに計算結果を詰めていきます。
ソースコードにコメントで番号をふったので、以下その番号に沿って説明していきます。

(1)関数オブジェクトfの作成

ここでは、として機能する関数オブジェクトfを作成しています。一行ですぱっと書く方法が分からなかったので、まずを定義しています。ソースコード上でgはカリー化された関数として定義しており、g(a)としたときには、t=aと固定し1つの引数t_oldのみ受け取る新たな関数が返ってくることに注意してください。そして、次にとするため、g(t)として返ってきた関数をGetMaxに渡しています。
関数GetMaxのソースコードは以下の通りです。これは関数と探索する値のリストを引数に与え、その探索範囲で関数の最大値を見つけるものです。初期化のところでも述べたように、はとのタプルを返すようにしたいので、GetMax最大値と最大にする変数の値のタプルを返すようにしています。

def GetMax(func , SearchVars):
    max_var = SearchVars[0]
    max_value = func(SearchVars[0])
    for var in SearchVars:
        if func(var) >  max_value:
            max_value =  func(var)
            max_var = var
    return max_value , max_var

(2)関数のメモ化

この実装では関数のメモ化というテクニックを利用して、動的計画法の肝となる途中までの計算結果の保持を各にやらせています。メモ化とは1度引数として与えて計算した値は保持しておいて、再度同様の引数が与えられたときに再計算する必要をなくすような関数にすることをいいます。メモ化に関してはM.Hiroiさんのホームページを参考にしました。
本実装では引数で与えた関数をメモ化してくれる関数をMemoizeFuncと定義しています。そのソースコードを以下に示しますが、ほとんど先ほど挙げた参考ページのまるパクリなのでそちらをみた方がいいと思います。

def MemoizeFunc(func):
    table = {}
    print table
    def memoized_func(*args):
        if not args in table:
            table[args] = func(*args)
        return table[args]
    return memoized_func

(3)再帰呼び出し

にあたる関数オブジェクトも作成したので、それを引数に与えてDPSegmentationを再帰呼び出しして次のステップに移ります。当然、fの添え字番号にあたるiもインクリメントして与えます。

(4)再帰の終着点

i==k-1ということは のところまで来たということになります。ここにきて初めて関数GetMaxが呼ばれます。その時、i==k-1のときのfの戻り値が必要なので、そのfが呼び出されますが、また、GetMax関数があり更にその中でi==k-2のときのfが呼び出され・・・という動作が繰り返され、最終的にはi==1のとき、つまりf_iniで定義した関数が呼び出されるという流れになります。従って、実質的な計算はこの段階に来て初めて行われ、それまでの過程は計算のための準備（関数オブジェクトの作成、メモ化）といえます。
後は、ResultListに各の値を詰めていけば終了です。

時系列データのクラス

上述の通り、SequenceSegmentationに渡すSequenceDataはGetLength()とh(t1,t2)のメソッドが定義されているオブジェクトである必要があります。また、今回は実数値の時系列データを取り扱うので、__init__()でnumpy.array型の実数値配列データを受け取り、これを時系列データの実際の数値とします。ということで時系列データクラスの定義は以下の通りです。

class TimeSeriesData:
    def __init__(self , Data):
        self.Data = Data

    def GetLength(self):
        return self.Data.shape[0]

    def h(self , t1 , t2):
        S = self.Data[t1:t2] - np.mean(self.Data[t1:t2])
        return  -np.sum(S*S)/(t2 - t1)

スコアh(t1 , t2)は、セグメント内の分散にマイナスをつけたものとしました。これで、このクラスから生成したオブジェクトをSequenceSegmentationに突っ込めば、時系列データをセグメント内の平均値で近似した時に誤差が最も小さくなるようなセグメンテーションが得られます。

実行と結果

以下のように、インタプリタに打ち込んで実行しました。

>>> data = np.sin(np.arange(0, np.pi * 4 ,(np.pi * 4)/150.0)) + 0.25 * np.random.randn(150)
>>> Sequence = TimeSeriesData(data)
>>> T = SequenceSegmentation(Sequence , 5)

データとしては2周期分の正弦波にテキトーなノイズをのせ、150点サンプリングしたものにしました。それを5つのセグメントに分割するように実行しています。ということで結果は以下の通り。元の時系列データの各点を青で示し、各セグメントを平均値で近似した曲線は赤で示しています。黒の線はセグメントの境界です。

大体、元々の正弦波が0と交差する付近に各セグメントの境界が来ています。流石にセグメンテーション数なので、近似という意味ではいまいちかもしれませんが、時系列的な推移の概略は分かるものになっていると思います。

スコアの計算方法をかえるだけで様々なセグメンテーションが可能

先ほどの例では、セグメント内の分散にマイナスをつけたものをセグメント内のスコアとしましたが、別の時系列データクラスを定義し、そのh(t1 , t2)を新たに定義することで、また別のセグメンテーションが得られます。
例えば、以下のコードはセグメントの境界位置にあたる点同士を直線で結び、その直線と実際の点の誤差の二乗平均にマイナスつけたものをスコアにしています。こうすることで、セグメントの境界の点同士を直線補間して得られる曲線が最も元の時系列データを近似するようなセグメントの境界位置を得ることができます。

class TimeSeriesLerp(TimeSeriesData):

    def AproximateSubSequence(self , t1 , t2):
        if t2 < self.GetLength():
            x1 = t1
            x2 = t2
        else:
            x1 = t1
            x2 = t2 - 1
        a = (self.Data[x2] - self.Data[x1]) / (x2 - x1)
        b = (x2 * self.Data[x1] - x1 * self.Data[x2])/(x2 - x1)
        return (a * np.arange(t1 , t2) + b)

    def h(self , t1 , t2):
        S = self.Data[t1:t2] - self.AproximateSubSequence(t1 , t2)
        return  -np.sum(S*S)/(t2 - t1)

面倒だったので、先ほどのTimeSeriesDataクラスを継承してhをオーバーライドしています。AproximateSubSequenceで直線補間値を得て、誤差の二乗平均を計算します。プログラム上、最後の境界位置をとしていますが、Data[N]はIndexErrorになるので、最後のセグメントだけ、セグメントの始点とN-1の点を結ぶようにするためif文で分岐させています。そしてこれを下記のようにしてセグメンテーションします。

>>> Sequence = TimeSeriesLerp(data)
>>> T = SequenceSegmentation(Sequence , 5)

dataは上述と同じ正弦波にノイズを乗せたものです。そしてセグメンテーションの結果は以下の通り。

平均で近似した時と違い、振幅の最大値・最小値付近に各セグメントの境界が来ています。近似と言う意味では、こちらの線形補間近似の方がより良く近似できていると思います。実際、得られたスコアもこちらの方が高い（つまり、誤差が小さい）です。

まとめ

ラムダ式やメモ化、再帰を利用すれば、かなりすっきりした形でシーケンシャルデータのセグメンテーションを実装できるなと感じました。ただ、これらに対する理解がまだまだ浅いので、無駄な事をしている箇所やもっと効率化できる部分があるかもしれません。そういった点があればご指摘いただければ幸いです。
また、今回は時系列データの近似のみを考えましたが、セグメント毎にスコアが得られるシーケンシャルデータであれば同様のアルゴリズムが適用できるので、なんか面白そうな題材があれば次は、それを試してみたいと思います。

*1:On the Approximation of Curves by Line Segments Using Dynamic Programming [R.Bellman , 1961]