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 \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などの潜在トピックモデルを用いたセグメンテーション手法も
いずれは試してみたいと思います。

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 \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$ として次の操作を収束するまで繰り返します。

$c_{s , t}$ の事後分布 $p(c_{s , t} | {\bf X}_{s ,t} , {\bf w}_s^{old} , \beta_s^{old})$ を求める
$\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)$ を求めて更新する
$\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自棄読み　その１

機械学習

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

Multi-Task Averaging

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

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

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

${\bf Y}^* = \left({\bf I} + \frac{\gamma}{T} \Sigma {\bf L}\right)^{-1} {\bf \bar{Y}}$

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

$a^* = \frac{2}{\Delta^2}$

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

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

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

$\displaystyle a^* = \frac{2}{ \frac{1}{T \left(T - 1 \right)} \sum_{r=1}^{T} \sum_{s=1}^{T} \left(\mu_r - \mu_s \right)^2}$

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

Hamming Distance Metric Learning

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

${\bf h} = b\left({\bf x};{\bf w}\right) = \rm{sign}\left(f\left({\bf x};{\bf w}\right)\right)$

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

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

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

$\displaystyle \max_{{\bf g} , {\bf g}^+ , {\bf g}^-} \left\{\ell_{{\rm triple}} \left({\bf g} , {\bf g}^+ , {\bf g}^-\right) + {\bf g}^T f\left({\bf x};{\bf w}\right) + {{\bf g}^+}^T f\left({\bf x}^+;{\bf w}\right) + {{\bf g}^-}^T f\left({\bf x}^-;{\bf w}\right)\right\} \\ \displaystyle \hspace{1em} - \max_{{\bf h}} {\bf h}^T f\left({\bf x};{\bf w}\right) - \max_{{\bf h}^+} {{\bf h}^+}^T f\left({\bf x};{\bf w}\right) - \max_{{\bf h}^-} {{\bf h}^-}^T f\left({\bf x};{\bf w}\right)$

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

$\displaystyle {\bf w}^{t+1} \leftarrow {\bf w}^{t} + \eta \left[\frac{\partial f\left({\bf x};{\bf w}\right)}{\partial {\bf w}} \left({\bf h} - {\bf g} \right) + \frac{\partial f\left({\bf x}^+;{\bf w}\right)}{\partial {\bf w}} \left({\bf h}^+ - {\bf g}^+ \right) + \frac{\partial f\left({\bf x}^-;{\bf w}\right)}{\partial {\bf w}} \left({\bf h}^- - {\bf g}^- \right) - \lambda {\bf w} \right]$

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

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

2013-01-27

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

アルゴリズム動的計画法

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

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

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

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

これを、上記のスコア最大化という文脈に当てはめると以下のような評価関数になります。
$\displaystyle J = \sum_{i=1}^{k} S_i = \sum_{i=1}^{k}\left( - \frac{1}{t - t_{i-1}}\sum_{t = t_{i - 1}}^{t_{i} - 1} \left(f(t) - \mu_i \right)^2 \right)$
ここでは評価関数の最大化という問題に当てはめたいので、誤差の二乗平均にマイナス付与したものをスコアとしています。ここでの時系列データは離散的にサンプリングされたデータであり、そのサンプル数は有限個（ $N$ 個）とするので、 $t_i$ は $0 \leq t_i < t_{i+1} \leq N , t_0 = 0 , t_k = N$ である整数です。この $J$ を最大化する $\{t_0 , t_1 , \cdots , t_k \}$ の組み合わせを求めるために、今回取り上げる動的計画法（Dynamic Programming, 以下DP）がよく利用されます。

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

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

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

作者: 金谷健一
出版社/メーカー: 共立出版
発売日: 2005/09/01
メディア: 単行本
購入: 29人クリック: 424回
この商品を含むブログ (42件) を見る

さて、各セグメントのスコア $S_i$ はセグメントを囲む境界位置 $t_{i-1},t_i$ の関数になります。そのため、 $S_i = h \left(t_{i-1} , t_i \right)$ とすると解きたい問題は以下のようになります。

$\displaystyle \hat{J} = \max_{t_1 , t_2 , \cdots , t_{k-1}}\left(h(0 , t_1) + h(t_1 , t_2) + h(t_2 , t_3) + \cdots + h(t_{k-1} , N) \right)$

${\rm s.t.} \hspace{1em} t_i \in {\mathbb N} \hspace{1em} , \hspace{1em} t_{i - 1} < t_i$

これは次のような再帰的な部分問題に分割することができます。
$f_1(t_1) = h(0 , t_1)$
$\displaystyle f_2(t_2) = \max_{t_1} \left(f_1(t_1) + h(t_1 , t_2)\right) \hspace{1cm} {\rm s.t.} \hspace{1em} 1 \leq t_1 < t_2$
$\vdots$
$\displaystyle f_i(t_i) = \max_{t_{i-1}} \left(f_i(t_{i-1}) + h(t_{i-1} , t_{i})\right) \hspace{1cm} {\rm s.t.} \hspace{1em} i - 1 \leq t_{i - 1} < t_i$
$\vdots$
$\displaystyle \hat{J} = \max_{t_{k - 1}} \left(f_{k - 1}(t_{k - 1}) + h(t_{k-1} , N)\right) \hspace{1cm} {\rm s.t.} \hspace{1em} k - 2 \leq t_{k - 1} < N$

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

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

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

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

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

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

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

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

この実装では関数のメモ化というテクニックを利用して、動的計画法の肝となる途中までの計算結果の保持を各 $f_i$ にやらせています。メモ化とは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)再帰呼び出し

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

(4)再帰の終着点

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

時系列データのクラス

上述の通り、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と交差する付近に各セグメントの境界が来ています。流石にセグメンテーション数 $k=5$ なので、近似という意味ではいまいちかもしれませんが、時系列的な推移の概略は分かるものになっていると思います。

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

先ほどの例では、セグメント内の分散にマイナスをつけたものをセグメント内のスコアとしましたが、別の時系列データクラスを定義し、その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で直線補間値を得て、誤差の二乗平均を計算します。プログラム上、最後の境界位置を $t_k = N$ としていますが、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]

2013-01-01

ICML2012の論文をいくつか3行程度で紹介する

機械学習

あけましておめでとうございます。
新年早々初詣にも行かず、4ヶ月滞ってたブログを更新するのが僕です。
ということで、今更ながらピックアップしてたICML2012の論文を読んでみました。
タイトルの通り3行で概要を書いていこうと思います。

Multiple Kernel Learning from Noisy Labels by Stochastic Programming

学習用サンプルのラベルに間違いが含まれることを考慮したマルチカーネル学習の手法
サンプル $i$ のラベルが合っているか否かを表すランダムな変数 $\xi_i$ を組み込んだ最適化問題として定式化
予め仮定したノイズレベルの範囲内で $\xi_i$ によって起こる最悪のケース（損失関数が最大になるケース）を考え、ランダム変数のない最適化問題へ変換

Dimensionality Reduction by Local Discriminative Gaussians

2次分類器のクロスバリデーション誤差最小化の近似を目的関数とした線形教師あり次元削減手法
その分類器のクラス条件付き確率密度は、平均及び共分散行列がサンプルのk近傍で計算されるローカルなガウス分布
目的関数を操作し、求めたい次元削減行列の直交制約を導入することで、固有値分解で解ける問題に変換

Adaptive Regularization for Weight Matrices*1

ベクトル $\bf{q}$ と $\bf{p}$ の類似度を ${\bf q}^{T}{\bf W}{\bf p}$ とし、 $\left({\bf q} , {\bf p}^{+} , {\bf p}^{-}\right)$ という三つ組（ ${\bf p}^{+}$ の方が ${\bf p}^{-}$ より ${\bf q}$ との類似度が高い）が与えられる時、 ${\bf W}$ をオンラインで学習する手法を提案
この論文の著者らが2009年に提案した、求めたいパラメータを平均とするガウス分布をオンラインで推定していくAROWという手法の行列版
${\bf W}$ をベクトルに展開することでAROWを適用できるが、そのままだと推定する共分散行列が馬鹿でかくなるので、対角行列など共分散行列に制限を設けて適用したアルゴリズムを提供

Discriminative Probabilistic Prototype Learning

1サンプルに対し複数の特徴ベクトルが与えられるデータセットにBoW*2を適用し1つの特徴ベクトルに再構築する際、再構築後の特徴ベクトルが分類しやすくなるように各Wordのプロトタイプを決める手法
各特徴ベクトルがどのWordに属するかをハードに割り当てるのではなく、ソフトマックス関数を用いて確率で表すことにより微分可能にして、尤度が最大になるプロトタイプを勾配法で推定
各教師サンプルのラベルも、どのクラスにどれくらいの確率で属するかという数値で与えられる

Learning Task Grouping and Overlap in Multi-task Learning

各タスクのパラメータ $\bf{w}_t$ がタスク数より少ない基底ベクトルの線形和（ ${\bf w}_t = \bf{L}{\bf s}_t$ ： ${\bf L}$ は基底ベクトルを列ベクトルにとった行列）で表されるとしたマルチタスク学習手法
各 ${\bf s}_t$ の $l1$ ノルム正則化項と ${\bf L}$ のフロベニウスノルム正則化項を損失関数に加えた目的関数を解くことにより実現
${\bf s}_t$ の $l1$ ノルム正則化項の作用によって、 ${\bf s}_t$ はスパースになり無関係なタスク同士は共通する基底ベクトルを持たなくなる

以上、5本の論文を紹介しました。
3行でまとめるのは難しいですね。
かなり無理やりやってるので、1行が長文になってしまいました。
ということで、今年もよろしくお願いします。

*1:ICML2012のホームページではAdaptive Regularization for Similarity Measuresとなってる

*2:Bag of Words（Bag of Keypointsともいう）

2012-09-07

【Python】Heat Kernel行列をfor文を使わずに求める

ふと、過去の記事を眺めていたら、ESFSの記事でこんな記述があるのを発見。

高速化のためには、for文で計算している部分を行列計算になおすなど、
繰り返し構文をなるべく使わない工夫が必要になります。
実際にどのように実装したかは、気が向いたら書こうと思います。

Pythonのようなスクリプト言語の場合、繰り返し構文で計算するよりも、
numpyなどによる行列演算を用いて、一括で計算した方が速くなることが多々あります。
これはMATLABでも同じことが言えますし（最近のはJIT-Acceleratorがあるので
そこまで変わらないかも）、Rなんかも同様のようです。
ということで、今回はPythonでESFSやLE、LPP等で必要になるHeat Kernel行列を
（一部強引に）繰り返し構文を使わず実装したコードを紹介します。
まず、ここでいうHeat Kernel行列は以下のようなものを指します。

$N$ 個の $d$ 次元ベクトル ${\bf x}_1 , {\bf x}_2 , \cdots , {\bf x}_N$ があり、 ${\bf x}_i$ が ${\bf x}_j$ のn近傍、または、 ${\bf x}_j$ が ${\bf x}_i$ のn近傍の場合、
$W_{ij} = e^{-\frac{\parallel {\bf x}_i - {\bf x}_j \parallel_2^2}{t}}$
それ以外の場合は、 $W_{ij} = 0$ となる $W_{ij}$ を各要素として持つ $N \times N$ の対称行列 ${\bf W}$ 。

ということで、この行列をつくるためにやらければならないことは、

各データ間の距離の2乗を計算
その距離に基づいて、各データのn近傍のサンプルを導出
近傍関係から接続行列（n近傍なら接続）を導出
データ間の距離と接続行列から、Heat Kernel行列を計算

この4つです。
これらの計算にfor文を一切使わないことで、高速化に努めました。

前提

下記のソースコードは、以下のようにインポートしていることを前提とします。

import numpy as np
import scipy as sci
import scipy.linalg

各データ間の距離の2乗を計算

ここで求めるのは次のような距離の2乗を要素として持つ行列です。
${\bf {\rm DistMat}} = \begin{equation} \begin{pmatrix} \parallel {\bf x}_1 - {\bf x}_1 \parallel_2^2 & \parallel {\bf x}_1 - {\bf x}_2 \parallel_2^2 & \cdots & \parallel {\bf x}_1 - {\bf x}_N \parallel_2^2 \\ \parallel {\bf x}_2 - {\bf x}_1 \parallel_2^2 & \parallel {\bf x}_2 - {\bf x}_2 \parallel_2^2 & \cdots & \parallel {\bf x}_2 - {\bf x}_N \parallel_2^2 \\ \vdots & \vdots & \ddots & \vdots \\ \parallel {\bf x}_N - {\bf x}_1 \parallel_2^2 & \parallel {\bf x}_N - {\bf x}_2 \parallel_2^2 & \cdots & \parallel {\bf x}_N - {\bf x}_N \parallel_2^2 \end{pmatrix} \end{equation}$
単純にループを回して、それぞれで距離を計算すれば、各データ間の距離は得られますが、
その実装方法だと $O(N^2)$ 回の繰り返しなので、 $N$ が大きくなればかなり遅くなります。
そのため、ここがHeat Kernel行列の計算の実装において、最も工夫すべき箇所だといえます。
さて、上記の行列ですが、次のように分解できます。
$\begin{eqnarray} {\bf {\rm DistMat}} &=& \begin{pmatrix} \parallel {\bf x}_1 \parallel_2^2 - 2 {\bf x}_1^T {\bf x}_1 + \parallel {\bf x}_1 \parallel_2^2 & \cdots & \parallel {\bf x}_1 \parallel_2^2 - 2 {\bf x}_1^T {\bf x}_N + \parallel {\bf x}_N \parallel_2^2 \\ \parallel {\bf x}_2 \parallel_2^2 - 2 {\bf x}_2^T {\bf x}_1 + \parallel {\bf x}_1 \parallel_2^2 & \cdots & \parallel {\bf x}_2 \parallel_2^2 - 2 {\bf x}_2^T {\bf x}_N + \parallel {\bf x}_N \parallel_2^2 \\ \vdots & \ddots & \vdots \\ \parallel {\bf x}_N \parallel_2^2 - 2 {\bf x}_N^T {\bf x}_1 + \parallel {\bf x}_1 \parallel_2^2 & \cdots & \parallel {\bf x}_N \parallel_2^2 - 2 {\bf x}_N^T {\bf x}_N + \parallel {\bf x}_N \parallel_2^2 \end{pmatrix} \\ & = & \begin{pmatrix}\parallel {\bf x}_1 \parallel_2^2 & \cdots & \parallel {\bf x}_1 \parallel_2^2 \\ \parallel {\bf x}_2 \parallel_2^2 & \cdots & \parallel {\bf x}_2 \parallel_2^2 \\ \vdots & \ddots & \vdots \\ \parallel {\bf x}_N \parallel_2^2 & \cdots & \parallel {\bf x}_N \parallel_2^2 \end{pmatrix} - 2 \begin{pmatrix} {\bf x}_1^T {\bf x}_1 & \cdots & {\bf x}_1^T {\bf x}_N \\ {\bf x}_2^T {\bf x}_1 & \cdots & {\bf x}_2^T {\bf x}_N \\ \vdots & \ddots & \vdots \\ {\bf x}_N^T {\bf x}_1 & \cdots & {\bf x}_N^T {\bf x}_N \end{pmatrix} + \begin{pmatrix}\parallel {\bf x}_1 \parallel_2^2 & \cdots & \parallel {\bf x}_N \parallel_2^2 \\ \parallel {\bf x}_1 \parallel_2^2 & \cdots & \parallel {\bf x}_N \parallel_2^2 \\ \vdots & \ddots & \vdots \\ \parallel {\bf x}_1 \parallel_2^2 & \cdots & \parallel {\bf x}_N \parallel_2^2 \end{pmatrix} \\ &=& {\bf H}^T -2{\bf X}^T {\bf X} + {\bf H} \end{eqnarray}$
ただし、
$\begin{eqnarray}{\bf H} = \begin{pmatrix}\parallel {\bf x}_1 \parallel_2^2 & \cdots & \parallel {\bf x}_N \parallel_2^2 \\ \parallel {\bf x}_1 \parallel_2^2 & \cdots & \parallel {\bf x}_N \parallel_2^2 \\ \vdots & \ddots & \vdots \\ \parallel {\bf x}_1 \parallel_2^2 & \cdots & \parallel {\bf x}_N \parallel_2^2 \end{pmatrix} & , & {\bf X} =\begin{pmatrix} {\bf x}_1 {\bf x}_2 \cdots {\bf x}_N \end{pmatrix} \end{eqnarray}$

後は、numpyのnp.diagやnp.tile（MATLABならrepmat）を利用して、
行列 $\bf{H}$ を生成すれば、繰り返し構文なしで
距離の2乗の行列を得ることができます。
サンプルの行列 $\bf{X}$ を引数にとり、
距離の2乗の行列を返す関数のソースコードは以下の通りです。

def GetDistanceMatrix(X):
    H = np.tile(np.diag(np.dot(X.T , X)) , (np.size(X , 1) , 1))
    G = np.dot(X.T , X)
    DistMat = H - 2 * G + H.T
    return DistMat

この関数に3次元1000個のデータ（3行1000列の行列）を入力したところ、
0.05秒程度で計算できるのに対し、forループを回して愚直に計算する関数の場合、
13秒程度かかります。
やはり、行列演算を用いて一括で計算した方が、格段に速いです。

距離行列に基づいて、各データのn近傍のサンプルを導出

np.argsortを用いて、各データにおける距離が小さいサンプル上位n件の
インデックスを求めるだけです。早速ソースコードです。

def GetnNearests(X , n):
    DistMat = GetDistanceMatrix(X)
    NearestsIdx = np.argsort(DistMat , axis=1)
    return NearestsIdx[: , 1:(n + 1)] , DistMat

近傍関係から接続行列（n近傍なら接続）を導出

上記のGetnNearests関数で得られるn近傍のインデックスから、
近傍関係なら1（接続）、そうでないなら0という要素で構成される行列を求めます。
ただし、ここでは対称行列がほしいので、どちらか一方から見て近傍関係に当たるなら、
このペアは接続ということにします。

def GetAdjacencyMatrix(X , n):
    (Idx , DistMat) = GetKNearests(X , n)
    AdjacencyMatrix = np.zeros((X.shape[1] , X.shape[1]) , 'bool');
    AdjacencyMatrix[np.tile(np.arange(Idx.shape[0]) ,
            (Idx.shape[1] , 1)).flatten('F') , Idx.flatten('C')] = True;
    AdjacencyMatrix |= AdjacencyMatrix.T #転置とorをとって対称行列にする
    return AdjacencyMatrix , DistMat

for文を使わない実装にこだわったため、実はここが一番苦労しました。
どうもスマートな方法が思いつかず、np.tileやflattenを駆使して、
1としたい要素のインデックスを無理やり並べて指定しています。
苦労した割にfor文を使用した実装と比べて、サンプル数1000程度では、
さほど速くなるわけではないので、ここは素直にfor文使った方が良いと思います。
もし、もっと良い方法があれば、教えてください。

データ間の距離と接続行列から、Heat Kernel行列を計算

ここまでくれば後は特別なことはしません。素直に実装します。

def GetHeatKernel(X , n , t):
    (AdjMat , DistMat) = GetAdjacencyMatrix(X , n)
    W = np.exp(- (DistMat / t) ) * AdjMat
    return W

np.expにnp.ndarrayのインスタンスを入力すると、
各要素のexpを返してくれます。
また、np.ndarrayクラス同士を*演算子で計算すると、
行列の積でなく、行列の各要素同士を掛け算した行列が返ってきます。

ついでだからLEを実装

Laplacian Eigenmapsです。LEについての解説は論文に譲ります。
Heat Kernel行列さえ求めれば、あとはグラフラプラシアン行列 $\bf{L}$ と
重みの対角行列 $\bf{D}$ 求めて、一般化固有値問題を解くだけです。

def LaplacianEigenmaps(X , n , t , m):
    W = GetHeatKernel(X , n , t)
    D = np.diag(np.sum(W , 0)) #重みの対角行列
    L = D - W #グラフラプラシアン
    (Lambda , V) = sci.linalg.eigh(L , D)
    Y = V[: , 1:(m + 1)].T
    return Y

おなじみの3次元スイスロールデータセット（サンプル数1000個）を2次元に次元削減してみました。

元データ

次元削減後（n = 10 , t = 5）

どうも、元論文のようにはならない・・・
とはいえ、他の論文をみるとスイスロールデータをLEで2次元にマッピングした結果が、
僕の結果と類似しているものもあるので、とりあえずよしとします。

どうやったら、ループ使わずに実装できるかと頭ひねるのも、なかなか楽しいですね。

2012-06-17

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

以前の記事でエンベを受験してきたことを報告しましたが、その合格発表が6月15日にありました。

結果は合格ということで一安心。
ただ、午後1，2ともに、自己採点の点数より思いのほか低く、合格ラインぎりぎりだったので少し焦りました。
たぶん、自分が「模範解答と多少違うけど、意図してることは同じだから○だろう」や
「少し、記述が足りないが大筋はあってるから、部分点くらいはもらえるだろう」などと勘ぐってたところが
実際は、×にされたのかなと思います。
何はともあれ、受かっていたのでえがったえがった