甲斐性なしのブログ

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

【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の元論文ではなく、この本を参考にセグメンテーションを実装しています。

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

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

さて、各セグメントのスコア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]

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ともいう）

【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次元にマッピングした結果が、
僕の結果と類似しているものもあるので、とりあえずよしとします。

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

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

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

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

NIPS2011斜め読み part2

機械学習

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

Heavy-tailed Distances for Gradient Based Image Descriptors

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

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

Efficient Methods for Overlapping Group Lasso

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

Transfer Learning by Borrowing Examples for Multiclass Object Detection

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

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

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

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

NIPS2011斜め読み

機械学習

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

Sparse Filtering

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

Learning a Distance Metric from a Network

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

Metric Learning with Multiple Kernels

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

Trace Lasso: a trace norm regularization for correlated designs

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

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

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

などなど

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

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

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

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

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

で、手ごたえですが

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

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

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

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

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

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

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