1 回帰とクラス分類の基礎

機械学習とは、\boldsymbol{y}=f(\boldsymbol{x})に従うデータ集合\left\{\boldsymbol{x}_n,\boldsymbol{y}_n\right\}に対し、関数fを推定するアルゴリズムの総称である。 第 1章では、機械学習の初歩的な事例として線型回帰最近傍法を実装し、機械学習の感覚的理解を目指す。

1.1 線型回帰

機械学習の中でも、標本\left\{\boldsymbol{x}_n,\boldsymbol{y}_n\right\}が入力\boldsymbol{x}に対する出力\boldsymbol{y}を明示的に指定する場合を教師あり学習と呼ぶ。 教師あり学習の中でも、出力すべき変数\boldsymbol{y}が連続値を取る場合を回帰と呼び、線型回帰はその初歩と言える。 線型回帰は、Eq. 1に示す通り、従属変数yを自由変数x_kに加重w_kを掛けた和で説明するモデルである。

y = f(\boldsymbol{x}) = w_0 + \displaystyle\sum_{k=1}^K w_k x_k = \begin{pmatrix} w_0 \\ w_1 \\ \vdots \\ w_K\end{pmatrix} \cdot \begin{pmatrix} 1 \\ x_1 \\ \vdots \\ x_K\end{pmatrix} = {}^t\!\boldsymbol{w}\ \boldsymbol{x}.\qquad(1)

また、何らかの非線型な写像\phi_k(\boldsymbol{x})の線型結合で回帰を行うEq. 2のモデルを線型基底関数モデルと呼ぶ。

y = f(\boldsymbol{x}) = \displaystyle\sum_{k=0}^K w_k \phi_k(\boldsymbol{x}) = {}^t\!\boldsymbol{w}\ \boldsymbol{\phi}(\boldsymbol{x}).\qquad(2)

Fig.  1はy\!\sim\!\mathcal{N}\left(y|x^3\!-\!3x,2500\right)に従う\left\{x_n,y_n\right\}\phi_k(x)\!=\!x^kの線型基底関数モデルを適用した例である。

Figure 1: linear basis function model.

線型基底関数モデルでは、他にも\phi_k(x)\!=\!\mathcal{N}\left(x|\mu_k,\sigma_k^2\right)を採用すれば、\mu_k近傍の局所的な曲線を表現できる。

y = f(\boldsymbol{x}) = \left\{ \begin{aligned} &\displaystyle\sum_{k=0}^K w_k x^k& &\enspace\mathrm{where}\enspace \phi_k(x)=x^k,\\ &\displaystyle\sum_{k=0}^K \displaystyle\frac{w_k}{\sqrt{2\pi\sigma_k^2}}\exp\left\{-\displaystyle\frac{(x-\mu_k)^2}{2\sigma_k^2}\right\}& &\enspace\mathrm{where}\enspace \phi_k(x)=\mathcal{N}\left(x|\mu_k,\sigma_k^2\right). \end{aligned} \right.

実際の観測にはノイズ\varepsilon_nが重畳する。ノイズは釣鐘型に分布するので、正規分布を仮定してEq. 3を得る。

y = f(\boldsymbol{x}) + \varepsilon \sim p(y|\boldsymbol{x},\boldsymbol{w},\beta) = \mathcal{N}\left(y|f(\boldsymbol{x}),\beta^{-1}\right).\qquad(3)

なお、\betaは分散の逆数で、ノイズの精度である。観測された標本\left\{\boldsymbol{x}_n,y_n\right\}の出現確率は、Eq. 4で求まる。

p(\left\{y_n\right\}|\left\{\boldsymbol{x}_n\right\},\boldsymbol{w},\beta) = \displaystyle\prod_{n=1}^N \mathcal{N}\left(y_n|{}^t\!\boldsymbol{w}\ \boldsymbol{\phi}(\boldsymbol{x}_n), \beta^{-1}\right).\qquad(4)

Eq. 4を尤度と呼ぶ。尤度は、与えられた訓練データを良く表現する母数としての加重\boldsymbol{w}の妥当性を示す。 線型回帰の学習は、尤度の最大化を目指す。今回は、計算の容易性から、Eq. 5の対数尤度を最大化する。

\log p(\left\{y_n\right\}|\left\{\boldsymbol{x}_n\right\},\boldsymbol{w},\beta) = \displaystyle\frac{N}{2}\log\beta - \displaystyle\frac{N}{2}\log2\pi - \displaystyle\frac{\beta}{2}\displaystyle\sum_{n=1}^N\{y_n-{}^t\!\boldsymbol{w}\ \boldsymbol{\phi}(\boldsymbol{x}_n)\}^2.\qquad(5)

Eq. 5から定数項を取り除き、加重\boldsymbol{w}に依存する項を取り出すと、Eq. 6の2乗誤差関数E(\boldsymbol{w})を得る。 Eq. 6の最小化は、最小二乗法とも呼ばれる。さて、誤差E(\boldsymbol{w})を最小化する具体的な方法を検討しよう。

E(\boldsymbol{w}) = \displaystyle\frac{1}{2} \displaystyle\sum_{n=1}^N \left\{y_n - \displaystyle\sum_{k=0}^K w_k \phi_k(\boldsymbol{x}_n)\right\}^2.\qquad(6)

誤差E(\boldsymbol{w})連続関数なので、勾配\nabla E(\boldsymbol{w})が0の点で最小になる。ここで、Eq. 7の勾配法を導入する。 勾配法とは、加重\boldsymbol{w}を極値に向けて動かす操作を繰り返し、最終的に\nabla E(\boldsymbol{w})=0で収束させる手法である。

\boldsymbol{w}' = \boldsymbol{w} - \eta \nabla E(\boldsymbol{w}) = \boldsymbol{w} + \eta \displaystyle\sum_{n=1}^N \{y_n - {}^t\!\boldsymbol{w}\ \boldsymbol{\phi}(\boldsymbol{x}_n)\} \boldsymbol{\phi}(\boldsymbol{x}_n).\qquad(7)

係数\eta学習率と呼ばれ、加重\boldsymbol{w}の発散を防止するため、\eta\nabla E(\boldsymbol{w})\ll\boldsymbol{w}となるように調整する必要がある。 勾配法による線型回帰をRegressionクラスに実装する。引数eは学習率\etaで、XYは標本\left\{x_n,y_n\right\}を表す。

class Regression(e: Double, XY: Seq[(Double, Double)], p: Seq[Double=>Double]) {
    val w = Array.fill[Double](p.size)(0)

配列pは、写像\phi_k(x)を格納する長さKの列\left\{\phi_k\right\}である。下記のapplyメソッドは、Eq. 2を計算する。

    def apply(x: Double) = w.zip(p.map(_(x))).map{case (w,x)=>w*x}.sum

最後に、Eq. 7による加重\boldsymbol{w}の更新を実装する。反復回数を固定せずに、収束判定を行うと実用的である。

    for(n<-1 to 1000; (x,y) <- XY; k<-0 until p.size) w(k) += e * (y-this(x)) * p(k)(x)
}

完成したRegressionクラスにy=x^3+\varepsilonの標本を与えて、基底\left\{x^3,x^2,x,1\right\}で線型回帰を行う例を示す。

val pts = -10.to(10).map(x=>(x.toDouble, math.pow(x,3) + util.Random.nextGaussian))
val reg = new Regression(0.000001, pts, 0.to(3).map(k=>(x: Double)=>math.pow(x,k)))

第 5章で学ぶニューラルネットワークは、活性化関数を追加した線型基底関数モデルと本質的に等価である。

1.2 最近傍法

教師あり学習の中でも、従属変数\boldsymbol{y}が離散値を取る場合をクラス分類と呼び、最近傍法はその初歩と言える。 Fig.  2 (a)に示す最近傍法は、未知の点を分類する際に近傍のk点を参考にする。kの決め方が重要である。

a b

Figure 2: k nearest neighbor model.. a — k nearest neighbor diameters., b — k\!=\!10 region segmentation.

最近傍法は遅延学習とも呼ばれ、事前の学習処理が不要で、分類時に初めて標本を参照する点が特徴である。 下記のKNNクラスは、未分類の点と訓練データの各点との距離を計算し、最近傍のk個の点で多数決を行う。

class KNN[D,T](k: Int, data: Seq[(D,T)], d: (D,D)=>Double) {

引数dには距離関数dを指定する。例えば、平方ユークリッド距離を使うには、下記のquad関数を与える。

def quad(a: Seq[Double], b: Seq[Double]) = (a, b).zipped.map(_-_).map(x => x * x).sum

距離関数dの決め方は重要な課題である。典型的には平方ユークリッド距離やマンハッタン距離を採用する。 厳密には、Eq. 8の距離の公理を満たす必要があるが、距離の比較ができれば、厳密な定義は不要である。

\left\{ \begin{aligned} d(\boldsymbol{x}, \boldsymbol{y}) &\geq 0, \\ \boldsymbol{x} = \boldsymbol{y} &\Leftrightarrow d(\boldsymbol{x}, \boldsymbol{y}) = 0, \\ d(\boldsymbol{x}, \boldsymbol{y}) &= d(\boldsymbol{y}, \boldsymbol{x}), \\ d(\boldsymbol{x}, \boldsymbol{y}) &\leq d(\boldsymbol{x}, \boldsymbol{z}) + d(\boldsymbol{z}, \boldsymbol{y}). \end{aligned} \right.\qquad(8)

下記のapplyメソッドは、座標\boldsymbol{x}近傍のk点を引数dataで与えた標本\left\{\boldsymbol{x}_n,y_n\right\}から探し、多数決を採る。

    def apply(x: D) = data.sortBy(s=>d(x,s._1)).take(k).groupBy(_._2).maxBy(_._2.size)._1
}

Fig.  2 (b)は、混合正規分布に従う標本をKNNクラスで学習して、変数\boldsymbol{x}の空間を塗り分けた結果である。 実用的には、距離計算の計算量を抑えるため、R木などの空間データベースで探索空間を局限すべきである。

2 決定木の学習と汎化性能

Eq. 9は、気象条件\boldsymbol{x}に対して質問と条件分岐を繰り返し、海水浴の是非yを判断する決定木の例である。

y \approx f(\boldsymbol{x}) = \begin{cases} 0 & \text{if $\mathrm{wavy}(\boldsymbol{x}) = 1$} \\ \text{otherwise} \left\{ \begin{aligned} & 0 && \text{if $\mathrm{rain}(\boldsymbol{x}) = 1$} \\ & 1 && \text{if $\mathrm{rain}(\boldsymbol{x}) = 0$} \\ \end{aligned} \right\} & \text{if $\mathrm{wavy}(\boldsymbol{x}) = 0$} \\ \end{cases}\qquad(9)

理想的な決定木は、簡潔明瞭である。故に、決定木の学習は標本\left\{\boldsymbol{x}_n\right\}に対する質問回数の最小化を目指す。 質問回数の最小化は、第 2.1に述べる情報量の加法性を勘案すれば、質問の情報利得の最大化と同義である。

2.1 情報利得の最大化

ある確率分布に従う記号列\left\{y_n\right\}を出力する装置を情報源と呼び、特定の情報yの価値I(y)情報量と呼ぶ。 真に価値ある情報は、稀有な筈である。試みに、情報量I(y)を記号yの出現確率P(y)の反比例で定義する。

I(y) = \displaystyle\frac{1}{P(y)}.\qquad(10)

Eq. 10の定義は自然に思えるが、記号列\left\{y_n\right\}の情報量I(\left\{y_n\right\})を計算すると、Eq. 11の違和感に気付く。

I(\left\{y_n\right\}) = \displaystyle\prod_{n=1}^N \displaystyle\frac{1}{P(y_n)} = \displaystyle\prod_{n=1}^N I(y_n).\qquad(11)

情報量が情報の価値を表す量であるなら、情報量I(\left\{y_n\right\})は、情報量\left\{I(y_n)\right\}の和になって然るべきである。 この直観を情報量の加法性と呼ぶ。そこで、対数関数を利用して、情報量I(y)の定義をEq. 12で修正する。

I(y) = \log_2 \displaystyle\frac{1}{P(y)} = - \log_2 P(y) \geq 0.\qquad(12)

また、情報源Yに対し、情報量I(y)の期待値H(Y)を定義する。期待値H(Y)エントロピーと呼ばれる。

H(Y) = \displaystyle\sum_{y \in Y} P(y) I(y) = - \displaystyle\sum_{y \in Y} P(y) \log_2 P(y) \geq 0.\qquad(13)

なお、エントロピーH(Y)0になる状況は、情報源Yが常に同じ記号yのみを出力する場合に限られる。 記号yの集合Yを質問Qで分割し、K個の部分集合\left\{Y_k\right\}を得た場合、Eq. 14のG(Q)を情報利得と呼ぶ。

G(Q) = H(Y) - H(Y|Q) = H(Y) - \displaystyle\sum_{k=1}^K \displaystyle\frac{|Y_k|}{|Y|} H(Y_k) \geq 0.\qquad(14)

集合Y_kは、決定木の子ノードに該当する。Y_kは更なる質問Q_kで分割されて、再帰的に木構造を構築する。 決定木を辿り、質問を重ねる度に、エントロピーは平均的に減少し、0に収束した時点でyの値が確定する。

trait Node[T] {
    def apply(x: Seq[Int]): T
}

決定木の動作は煩雑なので、1回の質問を複数のNodeの実装クラスに分解することで、実装を簡素化する。 Questionクラスは、引数Yに説明変数と従属変数の標本\left\{\boldsymbol{x}_n,y_n\right\}を受け取り、決定木を再帰的に構築する。

case class Question[T](Y: Seq[(Seq[Int], T)]) extends Node[T] {
    lazy val freqs = Y.groupBy(_._2).map(_._2.size.toDouble / Y.size)
    lazy val ent = freqs.map(f => -f * math.log(f)).sum / math.log(2)
    lazy val major = Y.groupBy(_._2).maxBy(_._2.size)._1
    lazy val v = Y.head._1.indices.map(Variable(Y, _)).minBy(_.t.ent)
    def apply(x: Seq[Int]) = if(ent - v.t.ent < 1e-5) major else v(x)
}

freqsはP(y)の値を記憶し、entはH(Y)の値を記憶する。majorはP(y)を最大化するyの値を記憶する。 applyメソッドは分類を行う。学習した質問に照らして、変数\boldsymbol{x}に対応する変数yの値を再帰的に推論する。

case class Variable[T](Y: Seq[(Seq[Int], T)], axis: Int) extends Node[T] {
    val t = Y.map(_._1(axis)).distinct.map(Division(Y, axis, _)).minBy(_.ent)
    def apply(x: Seq[Int]) = t(x)
}

Variableクラスの引数axisは分岐する軸を表し、その閾値はDivisionクラスの引数valueで指示する。 引数axisや引数valueには、Eq. 14の条件付きエントロピーH(Y|Q)を最小化する軸と閾値が選ばれる。

case class Division[T](Y: Seq[(Seq[Int], T)], axis: Int, value: Int) extends Node[T] {
    val sn1 = Question(Y.filter(_._1(axis) >  value))
    val sn2 = Question(Y.filter(_._1(axis) <= value))
    val ent = (sn1.ent * sn1.Y.size + sn2.ent * sn2.Y.size) / Y.size
    def apply(x: Seq[Int]) = if(x(axis) >= value) sn1(x) else sn2(x)
}

entはH(Y|Q)を保持する。applyメソッドは、条件分岐を行って、子ノードのapplyメソッドを実行する。 Fig.  3は、混合正規分布に従う標本をQuestionクラスで学習し、変数\boldsymbol{x}の空間を塗り分けた結果である。

Figure 3: region segmentation by a decision tree.

学習結果は過剰に複雑な境界線を描き、標本には忠実だが、母集団から乖離している。これを過学習と呼ぶ。 過学習を防ぐには、決定木の枝刈りを行うか、第 2.2で学ぶアンサンブル学習により汎化性能の改善を図る。

2.2 汎化誤差の最小化

標本\left\{\boldsymbol{x}_n,y_n\right\}から得た関数\hat{f}(\boldsymbol{x})と真のf(\boldsymbol{x})の間には汎化誤差が生じ、期待値を分解するとEq. 15を得る。

\int_{\boldsymbol{x}}P(\boldsymbol{x})(y-\hat{f}(\boldsymbol{x}))^2d\boldsymbol{x} = \mathbf{V}\!\left[\,y-f(\boldsymbol{x})\,\right] + \left(\underset{\hat{f}(\boldsymbol{x})}{\mathbf{E}}\!\left[\,\,\right] - f(\boldsymbol{x})\right)^2 + \mathbf{V}\!\left[\,\hat{f}(\boldsymbol{x})\,\right].\qquad(15)

T個の分類器\left\{\hat{f}_t(\boldsymbol{x})\right\}弱学習器と称して糾合し、投票で誤差の抑制を図る手法をアンサンブル学習と呼ぶ。

\hat{f}(\boldsymbol{x}) = arg\,max_k \displaystyle\frac{1}{T} \displaystyle\sum_{t=1}^T \mathbb{I}\left(\hat{f}_t(\boldsymbol{x})=k\right) \enspace\mathrm{where}\enspace \mathbb{I}\left(\hat{f}_t(\boldsymbol{x})=k\right) = \begin{cases} 1 & \text{if $\hat{f}_t(\boldsymbol{x})=k$}, \\ 0 & \text{if $\hat{f}_t(\boldsymbol{x})\neq k$}. \end{cases}\qquad(16)

中でもバギングと呼ばれる手法では、Eq. 17に示す相関を抑制することで、Eq. 15の第3項の抑制を図る。

\mathbf{V}\!\left[\,\hat{f}(\boldsymbol{x})\,\right] = \displaystyle\frac{1}{T^2}\displaystyle\sum_{i=1}^T \displaystyle\sum_{j=1}^T \mathbf{C}\!\left[\,\mathbb{I}\left(f_i(\boldsymbol{x})=k\right),\mathbb{I}\left(f_j(\boldsymbol{x})=k\right)\,\right].\qquad(17)

下記のBaggingクラスは、要素の重複を許して濃度Nの標本をT通り抽出し、Eq. 17の相関を抑制する。

case class Bagging[T](Y: Seq[(Seq[Int], T)], T: Int, N: Int) extends Node[T] {
    val t = Seq.fill(T)(Question(Seq.fill(N)(Y(util.Random.nextInt(Y.size)))))
    def apply(x: Seq[Int]) = t.map(_(x)).groupBy(identity).maxBy(_._2.size)._1
}

他方、ブースティングと呼ばれる手法では、弱学習器\hat{f}_t(\boldsymbol{x})\hat{f}_{t-1}(\boldsymbol{x})が判断を誤る点を重点的に学習する。 \hat{f}_t(\boldsymbol{x})\hat{f}_{t-1}(\boldsymbol{x})は対等でなく、その信頼度に基づく加重\left\{w_t\right\}を付与され、Eq. 18に示す加重投票を行う。

\hat{f}(\boldsymbol{x}) = arg\,max_k \displaystyle\sum_{t=1}^T w_t \mathbb{I}\left(\hat{f}_t(\boldsymbol{x}_n)=k\right).\qquad(18)

学習の目標は、Eq. 19に示す指数誤差の期待値を最小化する弱学習器\left\{\hat{f}_t(\boldsymbol{x})\right\}と加重\left\{w_t\right\}の設定にある。

\underset{\displaystyle\displaystyle\sum_{n=1}^N \exp \left\{-\displaystyle\displaystyle\displaystyle\frac{}{}{}{}{1}{K} \displaystyle\displaystyle\displaystyle\sum_{k=1}^K \mathbb{I}_{k}\!\left({y_n}\right) \mathbb{I}_{k}\!\left({\hat{f}(\boldsymbol{x})}\right)\right\}}{\mathbf{E}}\!\left[\,\,\right] \enspace\mathrm{where}\enspace \mathbb{I}_{k}\!\left({y}\right) = \begin{cases} 1& \text{if $y=k$},\\ \displaystyle\frac{1}{1-K}& \text{if $y\neq k$}. \end{cases}\qquad(19)

最後の弱学習器\hat{f}_T(\boldsymbol{x})の学習に着目し、標本の分布P_T(\boldsymbol{x},y)を導入して、Eq. 19をEq. 20に変形する。

\underset{\displaystyle\displaystyle\sum_{n=1}^N P_T(\boldsymbol{x}_n,y_n) \exp \left\{-\displaystyle\displaystyle\displaystyle\frac{}{}{}{}{1}{K} \displaystyle\displaystyle\displaystyle\sum_{k=1}^K \mathbb{I}_{k}\!\left({y_n}\right) w_T \mathbb{I}_{k}\!\left({\hat{f}_T(\boldsymbol{x}_n)}\right)\right\}}{\mathbf{E}}\!\left[\,\,\right].\qquad(20)

Eq. 19の変形の過程で現れる分布P_T(\boldsymbol{x},y)は、弱学習器\hat{f}_{T-1}(\boldsymbol{x})が判断を誤る点を重点的に学習させる。

P_T(\boldsymbol{x}_n,y_n) = \exp \left\{-\displaystyle\displaystyle\frac{}{}{1}{K} \displaystyle\displaystyle\sum_{k=1}^K \mathbb{I}_{k}\!\left({y_n}\right) \displaystyle\displaystyle\sum_{t=1}^{T-1} w_t \mathbb{I}_{k}\!\left({\hat{f}_t(\boldsymbol{x}_n)}\right)\right\}.\qquad(21)

Eq. 21に従う標本を抽出する操作は、ノイマンの棄却法を利用して、下記のResampleクラスに実装する。

case class Resample[T](Y: Seq[(Seq[Int], T)], P: Seq[Double]) extends Node[T] {
    def reject(i: Int) = if(util.Random.nextDouble * P.max < P(i)) Y(i) else null
    val data = new collection.mutable.ArrayBuffer[(Seq[Int], T)]
    while(data.size < P.size) data += reject(util.Random.nextInt(P.size)) -= null
    val quest = Question(data)
    val error = Y.zip(P).map{case ((x, y), p) => if(quest(x) != y) p else 0}.sum
    def apply(x: Seq[Int]) = quest(x)
}

下記のAdaStageクラスはM個の候補\left\{\hat{f}_{tm}(\boldsymbol{x})\right\}から誤差E_tが最小の候補を選び、弱学習器\hat{f}_t(\boldsymbol{x})とする。

case class AdaStage[T](Y: Seq[(Seq[Int], T)], P: Seq[Double], M: Int) extends Node[T] {
    val best = List.fill(M)(Resample(Y, P.map(_ / P.sum))).minBy(_.error)
    val W = math.log((1 / best.error - 1) * (Y.map(_._2).toSet.size - 1))
    def isOK = best.error < 0.5
    def apply(x: Seq[Int]) = best(x)
    def apply(x: Seq[Int], y: T): Double = if(best(x) == y) W else 0
    val next = Y.zip(P).map{case ((x, y), p) => p * math.exp(W - this(x, y))}
}

WはEq. 20を最小化する加重W_Tであり、Eq. 22で計算する。ただし、Kはクラスの異なり数である。

\hat{w}_T = \left\{\log \left(\displaystyle\displaystyle\frac{}{}{1}{E_T} - 1\right) + \log (K-1)\right\} \enspace\mathrm{where}\enspace E_T = \displaystyle\sum_{n=1}^N P(\boldsymbol{x}_n,y_n) \mathbb{I}\left(\hat{f}_T(\boldsymbol{x}_n) \neq y_n\right).\qquad(22)

下記のAdaBoostクラスは、弱学習器\hat{f}_t(\boldsymbol{x})の誤差E_t0.5を上回るまで\hat{f}_t(\boldsymbol{x})を生成し、加重投票を行う。

case class AdaBoost[T](Y: Seq[(Seq[Int], T)], M: Int) extends Node[T] {
    val stages = Seq(AdaStage(Y, Y.map(_ => 1.0 / Y.size), M)).toBuffer
    while(stages.last.isOK) stages += AdaStage(Y, stages.last.next, M)
    def apply(x: Seq[Int], y: T): Double = stages.init.map(_(x, y)).sum
    def apply(x: Seq[Int]) = Y.map(_._2).distinct.maxBy(this(x, _))
}

Eq. 19を最小化する\hat{f}(\boldsymbol{x})はEq. 23を満たす。故に、指数誤差の最小化はEq. 15の第2項を抑制する。

\hat{f}(\boldsymbol{x}) = arg\,min_k (K-1) \left\{\log P(y=k|\boldsymbol{x}) - \displaystyle\displaystyle\frac{}{}{1}{K} \displaystyle\displaystyle\sum_{k=1}^K \log P(y=k|\boldsymbol{x})\right\} = arg\,min_k P(y=k|\boldsymbol{x}).\qquad(23)

比較のため、Fig.  3と同じ標本をBaggingクラスとAdaBoostクラスに学習させた結果をFig.  4に示す。

a b

Figure 4: region segmentation by ensemble learning.. a — Bagging class., b — AdaBoost class.

バイアスバリアンス理論では、Eq. 15の汎化誤差の第2項をバイアスと呼び、第3項をバリアンスと呼ぶ。 決定木は、バイアスが低くバリアンスが高いモデルなので、ブースティングよりもバギングが効果的である。

2.3 情報源符号化と可逆圧縮

第 2章や第 5章には、エントロピーやKullback–Leibler 情報量などの情報理論に所縁のある概念が登場した。 付録 2.3では、情報理論を情報源符号化の観点から理解するため、可逆圧縮ツールhuffをD言語で実装する。

$ ./huff encode alice.txt bob.rz
$ ./huff decode bob.rz alice.txt

2.3.1 情報源符号化

確率P(s)に従う記号sの系列\left\{s\right\}を出力する情報源Sに対し、記号sの集合A_Sを情報源S字母と呼ぶ。 情報源Sの字母A_Sから別の字母A_Cへの写像C情報源符号と呼び、符号Cの値c\!\in\!C符号語と呼ぶ。

Table 1: constant length.
s c P(s) H(s)
A 00 0.5 1.00
B 01 0.2 2.32
C 11 0.2 2.32
D 10 0.1 3.32
Table 2: variable length.
s c P(s) H(s)
A 0 0.5 1.00
B 10 0.2 2.32
C 110 0.2 2.32
D 111 0.1 3.32

Table  1,  2は、符号の例である。符号語cの長さを符号語長と呼び、特に2進表記の場合にビット数と呼ぶ。 特にTable  1,  2 1に挙げた固定長符号の場合、符号語長L(c)の下限は字母A_Sの濃度の対数で決定される。

L(C) \ge \lceil \log_2 |A_S| \rceil = \log_2 4 = 2.

記号sの情報量I(s)に応じて符号語長L(c)を変化させる可変長符号は、平均符号語長\overline{L(C)}を削減できる。

\overline{L(C)} = \displaystyle\sum_{s \in A_S} P(s) L(C(s)) = 0.5 \times 1 + 0.2 \times 2 + (0.2 + 0.1) \times 3 = 1.8

Eq. 24に示すシャノンの情報源符号化定理は、可変長符号Cが達成できる\overline{L(C)}の理論的な下限を与える。

\overline{L(C)} \geq H(P) = -\displaystyle\sum_{s \in A_S} P(s) \log P(s).\qquad(24)

交差エントロピーH(P,Q)は、符号Cが想定する分布Q(s)が真の分布P(s)と異なる場合の下限を与える。 その際の余分な符号語長がKullback–Leibler 情報量D\!\left(P\|Q\right)であり、分布PQの差を測る尺度でもある。

H(P,Q) = -\displaystyle\sum_{s \in A_S} P(s) \log Q(s) = -\displaystyle\sum_{s \in A_S} P(s) \left\{\log P(s) - \log\displaystyle\frac{}{}{P(S)}{Q(S)}\right\} = H(P) + D\!\left(P\|Q\right).\qquad(25)

現実には、符号語長L(c)は整数に限られ、情報源Sの真の分布P(s)に従う可変長符号の実現は困難である。 実際の圧縮アルゴリズムは、Kullback–Leibler 情報量の現実的な低減を目指して様々な亜種が派生している。

2.3.2 ハフマン符号

ハフマン符号とは、記号sに対して、L(c)\!\sim\!I(s)となる最短の符号語cを充てる代表的な符号化方式である。 通常、ハフマン符号はFig.  5に示すハフマン木で表す。左右の分岐が符号語cの2進数の1桁に相当する。

Figure 5: Huffman coding tree.

符号語は、ハフマン木を巡回することで得られる。下記は、円周率1,000桁の符号語列の冒頭部の例である。

01111110 11011100 10100010 00101101 00001100 00011001 11010001 10100110
01110011 01010110 11100011 01100101 10101110 10000011 11101000 10011100
10110011 10101110 11000111 00100011 11100001 01111110 00001010 11111100
11101100 10011001 10000010 00100111 11111110 00110111 01110011 11111010
10001110 10101111 10011001 00001110 10100011 11110111 10000101 00000111
10001010 11011110 11111110 11110100 00101010 11100001 11111001 11010001
01011010 00101001 11111111 01110111 00111011 11110001 10011100 11001101
11111100 00000011 11100000 10100100 11101111 00100000 11000100 11001111
10011010 10001110 00011011 01101111 01100000 11111010 00111001 01111110

但し、実際には符号語の前にハフマン木の復元に使う情報を埋め込む必要がある。以下、実装例を掲載する。

import core.bitop,std.algorithm,std.bitmanip,std.conv,std.file;

まず、ハフマン木を格納する方法を検討する。専用のクラスを定義しても良いが、単なる配列で十分である。 N種類の記号を符号化する場合、ハフマン木の枝の数はN\!-\!1本であるので、ノードの数は2N\!-\!1個となる。

size_t parent[511];
size_t child0[511];
size_t child1[511];
size_t weight[511];

countSymbols関数は、符号化前のバイト列を引数に取り、各記号の出現回数を数えて正規化した値を返す。

ubyte[256] countSymbols(ubyte[] source) {
    size_t count[256];
    foreach(sym; source) count[sym]++;
    size_t maxim = reduce!(max)(count);
    foreach(s, cnt; count) if (cnt > 0) {
        double deg = 255.0 * cnt / maxim;
        count[s] = to!ubyte(max(deg, 1));
    } else count[s] = 0;
    return to!(ubyte[256])(count);
}

findRareSymbol関数は、前掲の配列weightを走査して、出現頻度が最低の2個のノードの組を探索する。

size_t findRareSymbol(size_t exclusion) {
  size_t idx = exclusion;
    size_t min = size_t.max;
  foreach(i, c; weight) if(c) {
    if(i == exclusion) continue;
    if(c < min) min = c,idx = i;
  }
  return idx;
}

addHuffmanTreeNode関数は、findRareSymbol関数が見つけたノードの組を纏めたノードを配列に加える。

bool addHuffmanTreeNode(const size_t n) {
  size_t i = findRareSymbol(510);
  size_t j = findRareSymbol(i);
  if(i == j && i > 255) return true;
  weight[n] = weight[i] + weight[j];
  weight[i] = weight[j] = 0;
  parent[i] = parent[j] = n;
  child0[n] = i;
  child1[n] = j;
  return false;
}

N種類の記号を符号化する際は、この操作をN\!-\!1回だけ繰り返すことで、完全なハフマン木を構築できる。 ASCIIを符号化する際には、N\!-\!1\!=\!255回の反復である。下記のcreateHuffmanTree関数の通り実装した。

size_t createHuffmanTree(ubyte[] count) {
    weight[0..256] = to!(size_t[])(count);
    size_t root = 0;
    foreach(n; 256..511) {
        if(!addHuffmanTreeNode(n)) root = n;
    }
    return root;
}

以上で、ハフマン木を構築する処理が完成した。下記のsymbolToCodeWord関数で、記号1件を符号化する。 これは、ハフマン木を末端から根に向け巡回することで、記号sに対応する符号語cの2進表現を取得する。

bool[] symbolToCodeWord(const size_t s) {
    bool[] code = new bool[0];
    for(long n=s; parent[n]; n=parent[n]) {
        code ~= child1[parent[n]] == n;
    }
    return code.reverse;
}

下記のencode関数は、復元時のバイト数と、各記号の出現頻度を出力したのち、符号語を順番に出力する。 なお、符号化に使ったハフマン木は復号時にも必要なので、countSymbols関数の返り値も冒頭に出力する。

void encode(string from, string target) {
    auto plain = cast(ubyte[]) read(from);
    createHuffmanTree(countSymbols(plain));
    ubyte[] coded = new ubyte[8];
    coded.append!(size_t)(plain.length);
    coded ~= countSymbols(plain);
    ubyte buf, i;
    foreach(s; plain) {
        foreach(b; symbolToCodeWord(s)) {
            if(b) buf |= 1; else buf &= ~1;
            if(++i % 8 != 0) buf <<= 1;
            else coded ~= buf, buf = 0;
        }
    }
    if(i % 8) coded ~= buf <<= 7 - i % 8;
    write(target, coded);
}

下記のdecode関数は、符号語列を1ビットずつ走査してハフマン木を辿り、符号化前の記号列を復元する。

void decode(string from, string target) {
    auto coded = cast(ubyte[]) read(from);
    auto size = coded[0..8].peek!size_t();
    auto freq = coded[8..264];
    size_t root = createHuffmanTree(freq);
    size_t s = root;
    ubyte[] plain;
    foreach(size_t bits; coded[264..$]) {
        foreach_reverse(i; 0..8) {
            const int bit = bt(&bits, i);
            s = (bit? child1: child0)[s];
            if(s >= 256) continue;
            plain ~= to!ubyte(s), s = root;
            if(--size == 0) break;
        }
    }
    write(target, plain);
}

最後に、main関数を定義する。huffの最初の引数がencodeの場合は圧縮して、decodeの場合は解凍する。

void main(string[] args) {
    if(args[1] == "encode") encode(args[2], args[3]);
    if(args[1] == "decode") decode(args[2], args[3]);
}

以上で、可逆圧縮ツールhuffが完成した。dmdを利用する場合、下記のコマンドによりコンパイルできる。

$ dmd huff -O -release -inline

なお、符号化したファイルの内容を確認する際には、xxd等の2進ダンプのツールを利用すると便利である。

$ xxd -b -g 1 -c 8 bob.rz

3 混合正規分布と最尤推定

第 3章では、観測データ\boldsymbol{x}が従う何らかの確率分布P(\boldsymbol{x})を仮定し、分布P(\boldsymbol{x})の母数を求める手順を学ぶ。 単峰型の分布としては正規分布が代表的だが、Fig.  6に示す混合正規分布なら多峰型の分布を表現できる。

Figure 6: Gaussian mixture model.

混合正規分布はK個の正規分布の線型和である。変数\boldsymbol{x}\!\in\!\mathbb{R}^Dは、確率w_kで正規分布\mathcal{N}\left(\boldsymbol{x}|\boldsymbol{\mu}_k,S_k\right)に従う。

P(\boldsymbol{x}) = \displaystyle\sum_{k=1}^K w_k \mathcal{N}\left(\boldsymbol{x}|\boldsymbol{\mu}_k,S_k\right) \enspace\mathrm{where}\enspace\displaystyle\displaystyle\sum_{k=1}^K w_k = 1.\qquad(26)

\boldsymbol{\mu}_kS_kk番目の正規分布の平均と分散共分散行列である。正規分布\mathcal{N}\left(\boldsymbol{x}|\boldsymbol{\mu},S\right)はEq. 27で定義される。

\mathcal{N}\left(\boldsymbol{x}|\boldsymbol{\mu},S\right) = \displaystyle\frac{1}{(\sqrt{2\pi})^D\sqrt{|S|}} \exp \left\{-\displaystyle\frac{1}{2} {}^t\!(\boldsymbol{x}-\boldsymbol{\mu})\ S^{-1} (\boldsymbol{x}-\boldsymbol{\mu})\right\}.\qquad(27)

母数\left\{w_k,\boldsymbol{\mu}_k,S_k\right\}の値を推定する前に、標本\left\{\boldsymbol{x}_n\right\}K個の部分集合に分類する問題を第 3.1節で検討する。

3.1 クラスタリング

相互に近接した2点が同じ部分集合に配属されるように標本\left\{\boldsymbol{x}_n\right\}を分割する操作をクラスタリングと呼ぶ。 直感的には、最適な部分集合C_kは、重心\boldsymbol{\mu}_kと内部の座標\forall\boldsymbol{x}\!\in\!C_kとの距離の総和を最小化する筈である。

\mathcal{D} = \displaystyle\sum_{n=1}^N \displaystyle\sum_{k=1}^K z_{nk} \|\boldsymbol{x}_n-\boldsymbol{\mu}_k\|^2.\qquad(28)

ただし、z_{nk}は点\boldsymbol{x}_nの所属を表す変数で、Eq. 29に定義する。z_{nk}は観測\left\{\boldsymbol{x}_n\right\}の埒外の潜在変数である。

\hat{z}_{nk} = \begin{cases} 1& \text{if $\boldsymbol{x}_n \in C_k$}\\ 0& \text{if $\boldsymbol{x}_n \not\in C_k$}. \end{cases}\qquad(29)

最適な\left\{z_{nk}\right\}\left\{\boldsymbol{\mu}_k\right\}は、反復法で求める。まず\left\{\boldsymbol{\mu}_k\right\}を乱数で初期化し、次にEq. 30で\left\{z_{nk}\right\}を更新する。

\hat{z}_{nk} = \begin{cases} 1& \text{if $k= arg\,min_{k}\|\boldsymbol{x}_n-\boldsymbol{\mu}_k\|^2$}\\ 0& \text{if $k\neqarg\,min_{k}\|\boldsymbol{x}_n-\boldsymbol{\mu}_k\|^2$}. \end{cases}\qquad(30)

次に、Eq. 31で\left\{\boldsymbol{\mu}_k\right\}を更新する。以後、Eq. 30の操作とEq. 31の操作を交互に反復し、収束解を得る。

\hat{\boldsymbol{\mu}}_k = \displaystyle\frac{1}{N_k} \displaystyle\sum_{n=1}^N z_{nk} \boldsymbol{x}_n \enspace\mathrm{where}\enspace N_k = \displaystyle\sum_{n=1}^N z_{nk}\; \leftarrow\displaystyle\frac{\partial }{\partial }{\mathcal{D}}{\boldsymbol{\mu}_k} = 2\displaystyle\sum_{n=1}^N z_{nk} (\boldsymbol{x}_n-\boldsymbol{\mu}_k) = 0.\qquad(31)

Eq. 30Eq. 31の反復で\left\{C_k\right\}を求める手法をk-meansと呼ぶ。下記のKmeansクラスはk-meansを実装する。

class Kmeans(X: Seq[Seq[Double]], K: Int, d: (Seq[Double],Seq[Double])=>Double) {
    val M = Array.fill(K, X.map(_.size).min)(Math.random)

標本\left\{\boldsymbol{x}_n\right\}は引数Xに渡す。引数dは距離関数であり、通常はEq. 28に従って、下記の2乗関数を与える。

    def quad(a: Seq[Double], b: Seq[Double]) = (a,b).zipped.map(_-_).map(d=>d*d).sum

下記のapplyメソッドは、指定した座標\boldsymbol{x}に対し、Eq. 30に従って、至近の部分集合C_kの番号kを返す。

    def apply(x: Seq[Double]) = M.map(d(_,x)).zipWithIndex.minBy(_._1)._2

下記のestepメソッドは、標本\left\{\boldsymbol{x}_n\right\}を部分集合\left\{C_k\right\}に分配し、Eq. 31に従って、重心\left\{\boldsymbol{\mu}_k\right\}を計算する。

    def estep = X.groupBy(apply).values.map(c=>c.transpose.map(_.sum / c.size))

最後に、Eq. 30Eq. 31の反復を実装する。回数は固定せずに、距離\mathcal{D}の収束を終了条件にすると完璧である。

    for(step<-1 to 100) (estep, M).zipped.foreach(_.copyToArray(_))
}

Fig.  7は、K\!=\!2の混合正規分布をKmeansクラスで分割した結果で、特に黒縁の2点は重心\left\{\boldsymbol{\mu}_k\right\}を表す。

Figure 7: k-means clustering on Gaussian mixture model.

概してクラスタリングは、標本\left\{\boldsymbol{x}_n\right\}が入力\boldsymbol{x}に対する出力\boldsymbol{y}を含まないので、教師なし学習に分類される。

3.2 期待値最大化法

後述の通り、k-meansは第 3.2節で学ぶ期待値最大化法の特殊例であり、母数の初期値を求める際に役立つ。 第 3.2節では、変数\left\{z_{nk}\right\}確率変数と見なし、重心\boldsymbol{\mu}_kに加えて加重w_kと分散S_kを推定する手順を学ぶ。

P(z_{nk}=1|\boldsymbol{x}_n) = \displaystyle\frac{ w_k \mathcal{N}\left(\boldsymbol{x}_n|\boldsymbol{\mu}_k,S_k\right)}{\displaystyle\displaystyle\sum_{k=1}^{K} w_k \mathcal{N}\left(\boldsymbol{x}_n|\boldsymbol{\mu}_k,S_k\right)} = \gamma_{nk}.\qquad(32)

Eq. 32に定義した\gamma_{nk}は、観測\boldsymbol{x}_nを得た後の潜在変数z_{nk}の確率分布を表すため、事後確率と呼ばれる。 以後、母数\left\{w_k,\boldsymbol{\mu}_k,S_k\right\}の最尤推定の式を導出する。まず、母数の値の妥当性を担保する尤度\mathcal{L}を定義する。

\mathcal{L}\left(\left\{w_k,\boldsymbol{\mu}_k,S_k\right\}\right) = P(\left\{\boldsymbol{x}_n\right\}|\left\{w_k,\boldsymbol{\mu}_k,S_k\right\}) = \displaystyle\prod_{n=1}^N P(\boldsymbol{x}_n|\left\{w_k,\boldsymbol{\mu}_k,S_k\right\}).\qquad(33)

Eq. 33にEq. 26を代入して、Eq. 34を得る。対数は、微小値の積による指数部のアンダーフローを防ぐ。

\log\mathcal{L}\left(\right){\left\{w_k,\boldsymbol{\mu}_k,S_k\right\}} = \log\displaystyle\prod_{n=1}^N \displaystyle\sum_{k=1}^K w_k \mathcal{N}\left(\boldsymbol{x}_n|\boldsymbol{\mu}_k,S_k\right) = \displaystyle\sum_{n=1}^N \log \displaystyle\sum_{k=1}^K w_k \mathcal{N}\left(\boldsymbol{x}_n|\boldsymbol{\mu}_k,S_k\right).\qquad(34)

以後、尤度\mathcal{L}の最大化を目指す。しかし、解析的な求解は困難である。Eq. 35に\boldsymbol{\mu}_kの偏微分の例を示す。

\displaystyle\frac{\partial }{\partial \boldsymbol{\mu}_k}\log\mathcal{L} = \displaystyle\frac{\partial }{\partial \boldsymbol{\mu}_k}\displaystyle\sum_{n=1}^N \log\displaystyle\sum_{k=1}^K w_k \mathcal{N}\left(\boldsymbol{x}_n|\boldsymbol{\mu}_k,S_k\right) = \displaystyle\sum_{n=1}^N \gamma_{nk} S_k^{-1} (\boldsymbol{x}_n - \boldsymbol{\mu}_k).\qquad(35)

Eq. 35の偏微分を0にする\boldsymbol{\mu}_kが、平均\boldsymbol{\mu}_kの推定値\hat{\boldsymbol{\mu}_k}となる。分散S_kの推定値\hat{S_k}も同様に計算する。

\left\{ \begin{aligned} \hat{\boldsymbol{\mu}_k} &= \displaystyle\frac{1}{N_k} \displaystyle\sum_{n=1}^N \gamma_{nk} \boldsymbol{x}_n,\\ \hat{S_k} &= \displaystyle\frac{1}{N_k} \displaystyle\sum_{n=1}^N \gamma_{nk} (\boldsymbol{x}_n - \boldsymbol{\mu}_k)\;{}^t\!(\boldsymbol{x}_n - \boldsymbol{\mu}_k)\ . \end{aligned} \right\} \enspace\mathrm{where}\enspace N_k = \displaystyle\sum_{n=1}^N \gamma_{nk}.\qquad(36)

分散\hat{S}_kは、変数\boldsymbol{x}の各次元の独立性を仮定する場合は、Eq. 37で代用できる。\circ要素ごとの積を表す。

\hat{S_k} = \displaystyle\frac{1}{N_k} \left(\displaystyle\sum_{n=1}^N \gamma_{nk} \boldsymbol{x}_n^{\circ 2}\right) - \boldsymbol{\mu}_k^{\circ 2}.\qquad(37)

加重w_kの推定値\hat{w_k}は、Eq. 26の制約条件を満たす必要があるので、ラグランジュの未定乗数法で求める。

\hat{w_k} = \displaystyle\frac{N_k}{N}.\qquad(38)

事後確率\left\{\gamma_{nk}\right\}が求まれば母数も求まるが、\left\{\gamma_{nk}\right\}の計算に\left\{w_k,\boldsymbol{\mu}_k,S_k\right\}が必要であり、求解は困難である。 そこで補助関数法を導入し、尤度\mathcal{L}の下限を与える補助関数Qの最大化により、間接的に\mathcal{L}を最大化する。

\log\mathcal{L}\left(\right){\boldsymbol{\theta}} = \max_{\boldsymbol{\gamma}} Q(\boldsymbol{\gamma}, \boldsymbol{\theta})\; \enspace\mathrm{where}\enspace \boldsymbol{\gamma} = \left\{\gamma_{nk}\right\},\;\boldsymbol{\theta}=\left\{w_k,\boldsymbol{\mu}_k,S_k\right\}.

補助関数Qは、Eq. 39の各等式の更新を交互に反復することで最大化でき、最終的には有限な値に収束する。

\begin{aligned} \hat{\boldsymbol{\gamma}}^{t+1} &= arg\,max_{\boldsymbol{\gamma}} Q(\boldsymbol{\gamma}, \boldsymbol{\theta}), \\ \hat{\boldsymbol{\theta}}^{t+1} &= arg\,max_{\boldsymbol{\theta}} Q(\boldsymbol{\gamma}, \boldsymbol{\theta}). \end{aligned}\qquad(39)

最尤推定に補助関数法を併用し、Eq. 39の反復により最尤推定を行う手法を期待値最大化法と呼ぶ。 ここでf(x)を凸関数、\left\{\gamma_n\right\}を正の実数列とすると、Eq. 40が成立し、これをイェンゼンの不等式と呼ぶ。

\displaystyle\sum_{n=1}^N \gamma_n f(x_n) \geq f\left(\displaystyle\sum_{n=1}^N \gamma_n x_n\right) \enspace\mathrm{where}\enspace \displaystyle\sum_{n=1}^N \gamma_n=1.\qquad(40)

\logが凹関数である点に注意して、Eq. 40にEq. 34を当て嵌めると、Eq. 41に示す補助関数Qを得る。 Eq. 39に関連して、補助関数Qを最大化する\hat{\boldsymbol{\gamma}}\hat{\boldsymbol{\theta}}を求めると、Eq. 32とEq. 36Eq. 38を得る。

\log\mathcal{L}\left(\right){\boldsymbol{\theta}} \geq \displaystyle\sum_{n=1}^N \displaystyle\sum_{k=1}^K \gamma_{nk} \log \displaystyle\frac{w_k \mathcal{N}\left(\boldsymbol{x}_n|\boldsymbol{\mu}_k,S_k\right)}{\gamma_{nk}} = Q(\boldsymbol{\gamma}, \boldsymbol{\theta}).\qquad(41)

期待値最大化法には別の解釈もある。仮に変数\left\{z_{nk}\right\}が観測可能ならば、対数尤度はEq. 42で定義される。 実際には、変数\left\{z_{nk}\right\}の真の値は観測不可能なので、変数\left\{z_{nk}\right\}の事後分布に関して尤度の期待値を考える。

\log P(\left\{\boldsymbol{x}_n,z_{nk}\right\}|\boldsymbol{\theta}) = \displaystyle\sum_{n=1}^N \displaystyle\sum_{k=1}^K z_{nk} \log \left\{w_k \mathcal{N}\left(\boldsymbol{x}_n|\boldsymbol{\mu}_k,S_k\right)\right\}.\qquad(42)

Eq. 39やEq. 32で\boldsymbol{\gamma}を求める操作は、Eq. 43に示す期待対数尤度の最新値を求める操作と等価であり、 Eq. 39やEq. 36Eq. 38で\boldsymbol{\theta}を求める操作は、その期待対数尤度を極大点まで近付ける操作と等価である。

\underset{\boldsymbol{z}}{\mathbf{E}}\!\left[\,\log P(\left\{\boldsymbol{x}_n,z_{nk}\right\}|\boldsymbol{\theta})\,\right] = \displaystyle\sum_{n=1}^N \displaystyle\sum_{k=1}^K \gamma_{nk} \log \left\{w_k \mathcal{N}\left(\boldsymbol{x}_n|\boldsymbol{\mu}_k,S_k\right)\right\}.\qquad(43)

これが期待値最大化法の名の由来である。 余談だが、分散をスカラー行列\varepsilon Eとすれば\varepsilon\!\to\!\inftyの極限でEq. 44が成立し、更に\gamma_{nk}\!\to\!z_{nk}も成立する。

\lim_{\varepsilon\to0} \varepsilon\log \left\{w \mathcal{N}\left(\boldsymbol{x}|\boldsymbol{\mu},\varepsilon E\right)\right\} = \lim_{\varepsilon\to0} \left\{\varepsilon\log w - \varepsilon\displaystyle\displaystyle\frac{}{}{}{}{D}{2} \log (2\pi\varepsilon) - \displaystyle\displaystyle\frac{}{}{1}{2} \|\boldsymbol{x}-\boldsymbol{\mu}\|^2\right\} = -\displaystyle\frac{1}{2} \|\boldsymbol{x}-\boldsymbol{\mu}\|^2.\qquad(44)

従って、Eq. 45が成立し、更に、期待値最大化法によるEq. 45の最大化はEq. 28の最小化に帰結する。

\varepsilon\underset{}{\mathbf{E}}\!\left[\,\,\right][\boldsymbol{z}]{\log P(\left\{\boldsymbol{x}_n,z_{nk}\right\}|\boldsymbol{\theta})} \to -\displaystyle\frac{1}{2} \displaystyle\sum_{n=1}^N \displaystyle\sum_{k=1}^K z_{nk} \|\boldsymbol{x}_n-\boldsymbol{\mu}_k\|^2 + C.\qquad(45)

期待値最大化法を実装する前にK個の正規分布\left\{\mathcal{N}\left(\boldsymbol{x}|\boldsymbol{\mu}_k,S_k\right)\right\}の混合分布を下記のGMMクラスに実装する。

class GMM(val D: Int, val K: Int) {
    val W = Array.fill(K)(1.0 / K)
    val M = Array.fill(K, D)(math.random)
    val S = Array.fill(K, D)(math.random)

配列WとMとSは、加重\left\{w_k\right\}と平均\left\{\boldsymbol{\mu}_k\right\}と分散\left\{S_k\right\}である。次のapplyメソッドはEq. 27を計算する。

    def apply(x: Seq[Double]) = for((w,m,s)<-(W,M,S).zipped) yield {
        val p = math.exp(-(x,m,s).zipped.map((x,m,s)=>(x-m)*(x-m)/s).sum/2)
        w * p / (math.pow(2 * math.Pi, .5 * m.size) * math.sqrt(s.product))
    }
}

下記のEMクラスは、期待値最大化法を実装する。標本\left\{\boldsymbol{x}_n\right\}は引数Xに、正規分布の個数は引数Kに渡す。

class EM(X: Seq[Seq[Double]], K: Int) {
    val gmm = new GMM(X.map(_.size).min, K)

下記のapplyメソッドは、事後確率\gamma_kの最大値を与えるkを返す。第 3.1節のapplyメソッドに相当する。

    def apply(x: Seq[Double]) = gmm(x).toSeq.zipWithIndex.maxBy(_._1)._2

次にM-射影を実装する。Eq. 38Eq. 36Eq. 37に従って母数\left\{w_k,\boldsymbol{\mu}_k,S_k\right\}を更新し、Eq. 34の尤度を返す。

    def mstep(P: Seq[Seq[Double]]): Double = {
        for(k<-0 until K) gmm.W(k) = P(k).sum / X.size
        val m1 = P.map((_,X).zipped.map((p,x)=>x.map(x=>p*x*1)).transpose)
        val m2 = P.map((_,X).zipped.map((p,x)=>x.map(x=>p*x*x)).transpose)
        for(k<-0 until K; d<-0 until gmm.D) {
            gmm.M(k)(d) = m1(k)(d).sum / P(k).sum
            gmm.S(k)(d) = m2(k)(d).sum / P(k).sum - math.pow(gmm.M(k)(d), 2)
        }
        X.map(x=>math.log(gmm(x).sum)).sum
    }

mstepメソッドの引数は、標本\left\{\boldsymbol{x}_n\right\}の事後確率\left\{\gamma_{nk}\right\}である。これは、Eq. 32に従って求める。 下記は、期待値最大化を交互に反復する。回数を固定せずに、尤度の収束を目安にすると実用的である。

    for(step <- 1 to 100) mstep(X.map(gmm(_)).map(p=>p.map(_/p.sum)).transpose)
}

Fig.  8は、Fig.  7と同じ標本をEMクラスで学習した結果で、等高線は、混合正規分布の密度関数を示す。 特に、 8 (a)は密度関数のヒートマップを表す。また、 8 (b)は混合正規分布によるクラスタリングの結果である。

a b

Figure 8: expectation maximization on a Gaussian mixture model.. a — density map., b — cluster map.

Fig.  7に比べ、正規分布が接する谷間に位置する標本は、正規分布の等高線を正しく反映した分類になる。

4 潜在的ディリクレ配分法

第 4章では、観測データの生成過程を確率変数の因果関係などの連鎖で表現するグラフィカルモデルを学ぶ。 自然言語処理に焦点を当て、百科事典の記事を訓練データに使う。まず、XMLのダンプデータを入手する。

$ wget https://dumps.wikimedia.org/jawiki/latest/jawiki-latest-abstract.xml

次に、を利用して、記事の題名と本文を抽出し、適当な配列に格納する。

val docs = scala.xml.XML.loadFile("jawiki-latest-abstract.xml") \\ "doc"
val data = Map(docs.map(d => (d\"title").text->(d\"abstract").text): _*)

記事には、文章を単語列に分解する形態素解析を予め施す。形態素解析器は、が便利だろう。

val gosen = net.java.sen.SenFactory.getStringTagger(null)
val words = gosen.analyze(data("Wikipedia: Scala"), new ArrayList[Token]())

助詞や助動詞や接続詞など機能語は、機械学習の際には無用のため、形態素解析の段階での排除を推奨する。 Eq. 46に示すtf-idfを指標に、特定の記事dで高頻度だが、他の記事で低頻度な単語t_dを残す方法もある。

\mathrm{tf\,idf}(t_d \in d,d \in D) = \mathrm{tf}(t_d,d) \mathrm{idf}(t_d,D)\; \enspace\mathrm{where}\enspace\left\{ \begin{aligned} \mathrm{tf}(t_d,d) &= \log \biggl(1+\displaystyle\frac{N_{t_dd}}{\displaystyle\displaystyle\sum_t N_{td}}\biggr) \\ \mathrm{idf}(t,D) &= \log \displaystyle\frac{|D|}{|\left\{d: t \in d\right\}|} \end{aligned} \right.\qquad(46)

ただし、N_{td}は文書dに単語tが出現する回数である。また、Eq. 46の対数は、\mathrm{tf}\mathrm{idf}の変動を抑制する。

4.1 単純ベイズ分類器

第 4.1節で学ぶ単純ベイズ分類器は、単語w_nの列である文書dが、話題cから生成された確率を推定する。

P(d|c) = P(w_1,\dots,w_{N_d}|c) = \displaystyle\prod_{n=1}^{N_d} P(w_n|c).\qquad(47)

Eq. 47のP(d|c)は、観測変数たる文書dが潜在変数たる話題cから生成された場合の、仮説の尤度を表す。 単語の共起を考慮するとEq. 48になるが、寧ろ独立性に拘るEq. 47は、この分類器の単純たる理由である。

P(w_1,\dots,w_{N_d}|c) = \displaystyle\prod_{n=1}^{N_d} P(w_n|c,w_1,...,w_{n-1}).\qquad(48)

文書dを生成した話題\hat{c}_dを推測する方法を検討する。直感的には、話題\hat{c}_dは確率P(c|d)の最大値を与える。

\hat{c}_d = arg\,max_{c \in C} P(c|d).\qquad(49)

確率P(d)が話題cに対し独立である点に留意して、Eq. 49にベイズの定理を適用すると、Eq. 50を得る。

\hat{c}_d = arg\,max_{c \in C} \displaystyle\frac{P(d|c)P(c)}{P(d)} = arg\,max_c P(d|c)P(c).\qquad(50)

独立性を仮定したEq. 47をEq. 50に代入すると、Eq. 51を得る。話題を判定する際は、Eq. 51を使う。

\hat{c}_d = arg\,max_{c \in C} P(c) \displaystyle\prod_{n=1}^{N_d} P(w_n|c).\qquad(51)

なお、確率P(w|c)は未知の単語wに対して0になる。対策として、Eq. 52のラプラス平滑化で補正する。 Eq. 52の確率P(w)は未知語w事前確率を、確率P(w|c)は訓練データを学習した後の事後確率を表す。

P(w|c) = \displaystyle\frac{N_{wc}+1}{\displaystyle\displaystyle\sum_{w\in V} (N_{wc}+1)}\; \leftarrow w \sim P(w) = \mathrm{Uniform}(V) = \displaystyle\frac{1}{|V|}.\qquad(52)

下記のNaiveBayesクラスは、観測dを得て潜在変数cの事後確率P(d|c)を最大化する話題\hat{c}_dを推定する。 型引数WとCは単語と話題を表す。引数docsとcatsには、文書\left\{d_i\right\}と話題\left\{c_i\right\}を同じ順番に並べて渡す。

class NaiveBayes[W,C](docs: Seq[Seq[W]], cats: Seq[C]) {
    val N = scala.collection.mutable.Map[(W,C),Double]().withDefaultValue(0)

下記の配列Pは、話題cの事前確率P(c)を保持する。単に、訓練データで話題\left\{c\right\}が出現した頻度である。

    val P = cats.groupBy(c=>c).map{case (c,s)=>c->s.size.toDouble/docs.size}

次に、単語wと話題cの共起の回数N_{wc}を配列Nに格納する。配列Nは未知語wに対してN_{wc}\!=\!0を返す。

    for((d,c)<-(docs,cats).zipped; w<-d) N(w,c) += 1

Pwcメソッドは、単語wを観測した場合の、その原因たる話題cの尤度P(w|c)をEq. 52に従って計算する。

    def Pwc(w: W, c: C) = (N(w,c)+1) / docs.flatten.distinct.map(N(_,c)+1).sum

Pcdメソッドは、文書dを観測した際の、その原因たる話題cの事後確率P(c|d)をEq. 51により計算する。

    def Pcd(c: C, d: Seq[W]) = math.log(P(c)) + d.map(w=>math.log(Pwc(w,c))).sum

最後にapplyメソッドを定義する。未知の文書dに対し、事後確率P(c|d)を最大化する話題\hat{c}_dを推定する。

    def apply(d: Seq[W]) = cats.distinct.maxBy(Pcd(_, d))
}

Fig.  9 (a)は、固有名詞を特徴量として東日本と西日本の記事を学習し、46都道府県を分類した結果である。 同様に、北海道と東北、関東と中部、近畿と中国、四国に九州の8地方に分類した結果をFig.  9 (b)に示す。

a b

Figure 9: Japanese map division into regions based on classification of Wikipedia pages.. a — 2-regional division., b — 8-regional division.

4.2 単語の話題の推定

第 4.2節で学ぶ潜在的ディリクレ配分法は、話題に応じて単語の意味が変化する複雑な言語モデルを扱える。 実際の自然言語の文では、短文でも複数の話題に言及し得るので、単語と意味を区別して考える必要がある。

I wrote a program drinking a cup of coffee.

Fig.  10に示す通り、複雑な言語モデルは単語wや話題zを始め、複数の確率変数の依存関係で表現される。

Figure 10: latent Dirichlet allocation model.

文書d\!\in\!Dに出現するn番目の単語w_{dn}の話題を変数z_{dn}で表す。両変数はEq. 53の多項分布に従う。 ただし、Vは語彙量を、Kは話題の数を表し、N_v\in\left\{\right\}{0,1}N_k\in\left\{\right\}{0,1}は単語vと話題kの出現を表す。

\begin{aligned} w_{dn} &\sim P(w|z) = \displaystyle\prod_{v=1}^V \phi_{zv}^{N_v} = \mathrm{Mul}\left(\boldsymbol{\phi}_z\right) \enspace\mathrm{where}\enspace \displaystyle\sum_{v=1}^V N_v = 1,\\ z_{dn} &\sim P(z|d) = \displaystyle\prod_{k=1}^K \theta_{dk}^{N_k} = \mathrm{Mul}\left(\boldsymbol{\theta}_d\right) \enspace\mathrm{where}\enspace \displaystyle\sum_{k=1}^K N_k = 1. \end{aligned}\qquad(53)

\phi_{zv}は話題zが単語vを生起する確率P(v|z)を表し、\theta_{dk}は文書dが話題kを生起する確率P(k|d)を表す。 単語wと文書dの結合確率はEq. 54で与える。話題z潜在変数で観測の埒外にあるので、周辺化する。

P(w,d) = P(d) \displaystyle\sum_{z=1}^K P(w|z) P(z|d) = \int_z \phi_z(w) \theta_d(z) dz.\qquad(54)

観測された文書dを的確に表す母数\phi_zや母数\theta_dを推定する場合は、Eq. 55の尤度関数の極値を探索する。

\mathcal{L}\left(\phi_z,\theta_d\right) = P(\boldsymbol{w}_d|\phi_z,\theta_d) = \displaystyle\prod_{n=1}^{N_d} \int_z \phi_z(w_{dn}) \theta_d(z) dz.\qquad(55)

ただし、母数\theta_dの個数は文書の量に比例し、膨大な数になるため、単なる最尤推定では過学習を頻発する。 対策としてベイズ推定を導入し、文書dを観測した後の母数\theta_dや母数\phi_zの分布を事後分布として計算する。

P(\phi_z,\theta_d|\boldsymbol{w}_d) = \displaystyle\frac{P(\boldsymbol{w}_d|\phi_z,\theta_d)P(\phi_z,\theta_d)}{P(\boldsymbol{w}_d)} \propto P(\boldsymbol{w}_d|\phi_z,\theta_d)P(\phi_z,\theta_d).\qquad(56)

Eq. 56のP(\phi_z,\theta_d)事前分布と呼び、文書dを観測する前に予想された母数\phi_zや母数\theta_dの分布を表し、 過学習を抑える働きを持つ。事後分布は尤度と事前分布の積に比例し、最適な母数は事後確率を最大化する。

\begin{aligned} \hat{\phi}_z &= arg\,max_{\phi_z} P(\phi_z,\theta_d|\boldsymbol{w}_d),\\ \hat{\theta}_d &= arg\,max_{\theta_d} P(\phi_z,\theta_d|\boldsymbol{w}_d). \end{aligned}

なお、Eq. 56は、1件の文書を観測する場合を想定した。文書が複数ある場合には、Eq. 57を適用する。

P(\phi_z,\theta_1,...,\theta_D|\boldsymbol{w}) \propto P(\phi_z,\theta_1,...,\theta_D) \displaystyle\prod_{d=1}^D P(\boldsymbol{w}_d|\phi_z,\theta_d).\qquad(57)

Eq. 57は、文書d_tを観測した後の事後分布を事前分布として次の文書d_{t+1}を学習する連鎖の様子を表す。 連鎖を容易にする目的から、事前分布はEq. 58のディリクレ分布で定義する。\Gamma\left(z\right)ガンマ関数である。

\mathrm{Dir}\left(\boldsymbol{\theta}|\boldsymbol{\alpha}\right) = \displaystyle\frac{\Gamma\left(\displaystyle\displaystyle\displaystyle\sum_{k=1}^K \alpha_k\right)}{\displaystyle\displaystyle\prod_{k=1}^K \Gamma\left(\alpha_k\right)} \displaystyle\prod_{k=1}^K \theta_k^{\alpha_k-1}\; \enspace\mathrm{where}\enspace \Gamma\left(z\right) = \int_0^\infty t^{z-1} e^{-t} dt.\qquad(58)

Eq. 53より、Eq. 55の尤度は多項分布を描くが、Eq. 59より、事後確率もディリクレ分布に従う。

\mathrm{Dir}\left(\boldsymbol{\theta}|\boldsymbol{n}+\boldsymbol{\alpha}\right) \propto \mathrm{Mul}\left(\boldsymbol{n}|\boldsymbol{\theta}\right) \mathrm{Dir}\left(\boldsymbol{\theta}|\boldsymbol{\alpha}\right).\qquad(59)

Eq. 59の性質により、事前確率は事後確率と共通の確率密度関数で定義でき、これを共役事前分布と呼ぶ。

\begin{aligned} \phi_k & \sim \mathrm{Dir}\left(\phi|\boldsymbol{\beta}\right),\\ \theta_d & \sim \mathrm{Dir}\left(\theta|\boldsymbol{\alpha}\right). \end{aligned}\qquad(60)

母数\phi_kと母数\theta_dの事前分布をEq. 60に設定し、Eq. 56のベイズ推定に組み込むと、Eq. 61を得る。

P(\boldsymbol{\phi},\boldsymbol{\theta}|\boldsymbol{w},\boldsymbol{\alpha},\boldsymbol{\beta}) \propto \displaystyle\prod_{k=1}^K \mathrm{Dir}\left(\phi_k|\boldsymbol{\beta}\right) \displaystyle\prod_{d=1}^D \mathrm{Dir}\left(\theta_d|\boldsymbol{\alpha}\right) \displaystyle\prod_{n=1}^{N_d} \int_z \phi_z(w_{dn}) \theta_d(z) dz.\qquad(61)

Eq. 61の事後確率を最大化する点を微分により求めれば、母数\phi_kや母数\theta_dの具体的な数値が得られるが、 微分は困難である。そこで、話題z_{dn}の標本を近似的に生成し、Eq. 62で母数を計算する方針を検討する。

\begin{aligned} \hat{\phi}_{kv} &= \displaystyle\frac{N_{kv} + \beta_v}{\displaystyle\displaystyle\sum_{v=1}^V \left(N_{kv} + \beta_v\right)},\\ \hat{\theta}_{dk} &= \displaystyle\frac{N_{dk} + \alpha_k}{\displaystyle\displaystyle\sum_{k=1}^K \left(N_{dk} + \alpha_k\right)}. \end{aligned}\qquad(62)

ただし、N_{kv}は話題kで単語vが出現した回数を表し、N_{dk}は文書dで話題kの単語が出現した回数を表す。 時刻tの標本\left\{z_{dn}^t\right\}は、Eq. 63の提案分布に従う乱数により生成する。これをギブスサンプリングと呼ぶ。

\forall z_{dn}^t \sim P(z_{dn}|z_{11}^t,...,z_{dn-1}^t,z_{dn+1}^{t-1},...,z_{DN_D}^{t-1}).\qquad(63)

変数\phi_zと変数\theta_dも潜在変数なので、本来は標本に含むべきであるが、Eq. 64の周辺化により除去できる。

P(\boldsymbol{z},\boldsymbol{w}|\boldsymbol{\alpha},\boldsymbol{\beta}) = \left(\displaystyle\prod_{k=1}^K \int \mathrm{Dir}\left(\phi_k|\boldsymbol{\beta}\right) \displaystyle\prod_{d=1}^D \displaystyle\prod_{n=1}^{N_d} \phi_{z_{dn}}\!(w_{dn})\;d\phi_k\right)\! \left(\displaystyle\prod_{d=1}^D \int \mathrm{Dir}\left(\theta_d|\boldsymbol{\alpha}\right) \displaystyle\prod_{n=1}^{N_d} \theta_d(z_{dn})\;d\theta_d\right).\qquad(64)

詳細は省略するが、ベータ関数の積分が出現する点に注目しつつ、Eq. 64を変形すると、Eq. 65を得る。

P(\boldsymbol{z},\boldsymbol{w}|\boldsymbol{\alpha},\boldsymbol{\beta}) = \displaystyle\prod_{k=1}^K \left\{ \displaystyle\frac{\Gamma\left(\displaystyle\displaystyle\displaystyle\sum_{v=1}^V\beta_v\right)}{\displaystyle\displaystyle\prod_{v=1}^V\Gamma\left(\beta_v\right)} \displaystyle\frac{\displaystyle\displaystyle\prod_{v=1}^V\Gamma\left(N_{kv}+\beta_v\right)}{\Gamma\left(N_k+\displaystyle\displaystyle\displaystyle\sum_{v=1}^V\beta_v\right)} \right\} \displaystyle\prod_{d=1}^D \left\{ \displaystyle\frac{\Gamma\left(\displaystyle\displaystyle\displaystyle\sum_{k=1}^K\alpha_k\right)}{\displaystyle\displaystyle\prod_{k=1}^K\Gamma\left(\alpha_k\right)} \displaystyle\frac{\displaystyle\displaystyle\prod_{k=1}^K\Gamma\left(N_{dk}+\alpha_k\right)}{\Gamma\left(N_d+\displaystyle\displaystyle\displaystyle\sum_{k=1}^K\alpha_k\right)} \right\}.\qquad(65)

Eq. 65で単語w_{dn}と話題z_{dn}に依存する部分をEq. 63の条件付き分布に代入すると、Eq. 66を得る。 ただし、Eq. 66で、\boldsymbol{w}^{\cancel{dn}}は単語w_{dn}を除外した単語の列であり、\boldsymbol{z}^{\cancel{dn}}\boldsymbol{w}^{\cancel{dn}}に対応する話題の列である。

P(z_{dn}|\boldsymbol{z}^{\cancel{dn}}) \propto P(z_{dn}=k,w_{dn}=v,\boldsymbol{z}^{\cancel{dn}},\boldsymbol{w}^{\cancel{dn}}|\boldsymbol{\alpha},\boldsymbol{\beta}) \propto \displaystyle\frac{N_{kv}^{\cancel{dn}} + \beta_v}{N_k^{\cancel{dn}} + \displaystyle\displaystyle\sum_{v=1}^V \beta_v} \displaystyle\frac{N_{dk}^{\cancel{dn}} + \alpha_k}{N_d^{\cancel{dn}} + \displaystyle\displaystyle\sum_{k=1}^K \alpha_k}.\qquad(66)

Eq. 66の提案分布に従う標本生成を何度も反復すれば、最終的に、標本\boldsymbol{z}は真の分布に収束する筈である。 下記のLDAクラスは、潜在的ディリクレ配分法を実装する。引数docsに文書の題名と単語列の配列を渡す。

class LDA[W,D](docs: Map[D, Seq[W]], K: Int, a: Double = 0.1, b: Double = 0.01) {

引数Kは話題の個数で、引数aとbは\forall a_k\!=\!a,\forall b_v\!=\!bを仮定した対称ディリクレ分布の母数\alpha\betaである。 LDAクラス内に、単語w_{dn}を包むWordクラスを実装する。乱数で初期化される変数zは、話題z_{dn}を表す。

    case class Word(v: W, var z: Int = util.Random.nextInt(K))

配列Wは集合\left\{w_{dn},z_{dn}\right\}に相当する。文書dと位置nを索引としてWordクラスのインスタンスを管理する。

    val W = docs.mapValues(_.map(Word(_)).toArray)

配列Vは語彙Vの要素v\!\in\!Vを索引として集合\left\{w_{dn},z_{dn}\right\}を管理する。変数N_{kv}の値を計算する際に使う。

    val V = W.flatMap(_._2).groupBy(_.v)

下記のNdkメソッドはEq. 62の変数N_{dk}に相当し、文書dに含まれる話題kの単語の出現回数を数える。

    def Ndk(d: D) = 0.until(K).map(k=>W(d).count(_.z==k) + a)

NkvメソッドはEq. 62の変数N_{kv}に相当し、単語の列\left\{w_{dn}\right\}の部分集合\left\{w_{dn}\!=\!v,z_{dn}\!=\!k\right\}の濃度を返す。

    def Nkv(v: W) = 0.until(K).map(k=>V(v).count(_.z==k) + b)

配列NkVはN_{kv}V全体に渡って合計した値である。実行効率の都合から、関数ではなく配列で実装した。

    val NkV = V.keys.map(Nkv).reduce((_,_).zipped.map(_+_)).toArray

次に、Eq. 63に従う乱数で標本\left\{z_{dn}^t\right\}を更新する処理を記述する。この処理は、収束まで何度も繰り返す。

    for(step<-1 to 1000; d<-docs.keys; (w,n)<-docs(d).zipWithIndex) {
        NkV(W(d)(n).z) -= 1
        W(d)(n).z = -1
        val S = (Nkv(w), Ndk(d), NkV).zipped.map(_*_/_).scan(.0)(_+_)
        val r = util.Random.nextDouble * S.last
        W(d)(n).z = S.tail.indexWhere(_ >= r)
        NkV(W(d)(n).z) += 1
    }

applyメソッドは、文書dで支配的な話題を探す。未知の文書には適用できず、専らクラスタリングに使う。

    def apply(d: D) = 0.until(K).maxBy(Ndk(d)(_))
}

固有名詞を特徴量として46都道府県の記事をLDAクラスで学習し、色を塗り分けた結果をFig.  11に示す。

a b

Figure 11: Japanese map division into regions based on clustering of Wikipedia pages.. a — 3-regional division., b — 5-regional division.

教師なし学習であるため第 4.1節と比べて分類結果の評価は難しく、変数N_{kv}を詳細に解析する必要がある。

5 ニューラルネットワーク

ニューラルネットワークは神経系を模倣した計算模型であり、神経細胞に相当する単位をニューロンと呼ぶ。

Figure 12: a neuron.

細胞はD本の樹状突起で信号\boldsymbol{x}を受容し、加重\boldsymbol{w}で合計した後、活性化関数fの値を軸索末端で放出する。

y = f(z) = f(\boldsymbol{w}\cdot\boldsymbol{x}) = f\left(\displaystyle\sum_{d=1}^D w_d x_d\right).\qquad(67)

関数fは、神経細胞が興奮する閾値に相当する。fが恒等関数f_\mathrm{id}(z)\!=\!zの場合、Eq. 67は線型回帰になる。 他方、関数fを階段関数やEq. 68のシグモイド関数で実装すると、Eq. 67はy\!\in\!\left\{0,1\right\}の2値分類になる。

f_\mathrm{sigm}(z) = \displaystyle\frac{1}{1 + e^{-z}} = \displaystyle\frac{1}{2}\tanh\displaystyle\frac{}{}{z}{2} + \displaystyle\frac{1}{2}.\qquad(68)

他にもFig.  13 (a)に示す\mathrm{tanh}関数やReLU関数が存在し、回帰やクラス分類など用途に応じて選択できる。 Fig.  13 (b)は、論理和をf=f_\mathrm{sigm}の分類器で学習した結果である。直線f(\boldsymbol{x})\!=\!0.5はクラスの境界を表す。

a b

Figure 13: neuron mechanism.. a — activation functions., b — OR learned by neuron.

分類器としてのニューロンは線型分類器と呼ばれ、境界が超平面になる線型分離可能な問題のみ学習できる。

5.1 誤差逆伝播法

Fig.  14は、ニューロンの多層化により非線型な回帰問題や分類問題に対応した多層パーセプトロンである。 左端で信号\boldsymbol{x}_1を受容する層を入力層、右端で信号\boldsymbol{x}_3を出力する層を{}と呼び、その間を{}と呼ぶ。

Figure 14: multi-layer perceptron.

信号\boldsymbol{x}は層を経る度に非線型に変換されるため、3層あれば任意の滑らかな関数を任意の精度で近似できる。 入力層と隠れ層の活性化関数をfとし、第m層の加重をW_mと表記する。出力\boldsymbol{y}はEq. 69で定義される。

\boldsymbol{y} = \boldsymbol{x}_3 = f(\boldsymbol{z}_2) = f(W_2\boldsymbol{x}_2) = f(W_2f(\boldsymbol{z}_1)) = f(W_2f(W_1\boldsymbol{x}_1)).\qquad(69)

Eq. 69の学習は、出力\left\{\boldsymbol{y}\right\}と正解\left\{\boldsymbol{t}\right\}の誤差を表す損失関数Eが最小になる加重Wの探索を通じて行う。 関数Eは、出力層の直前の関数fに応じて選ぶべきだが、Eq. 70の2乗誤差は、特定のfに拠らず使える。

E_\mathrm{sq}(\boldsymbol{y}\sim q(\boldsymbol{y}),\boldsymbol{t}\sim p(\boldsymbol{t})) = \displaystyle\frac{1}{2} \|\boldsymbol{y}-\boldsymbol{t}\|^2\;\enspace\mathrm{where}\enspace \|\boldsymbol{a}\| = \sqrt{|a_1|^2+\cdots+|a_d|^2}.\qquad(70)

Eq. 71の交差エントロピーE_\mathrm{CE}は分布pqの差を表すKullback–Leibler 情報量D\!\left(p\|q\right)の上限を与える。 誤差E_\mathrm{CE}の最小化はD\!\left(p\|q\right)の直接的な最小化に比べ、誤差逆伝播法に馴染む。詳細は第 5.3節で解説する。

E_\mathrm{CE}(p,q) = -\int p(\boldsymbol{y}) \log q(\boldsymbol{y}) d\boldsymbol{y} = -\int p(\boldsymbol{y}) \left\{\log p(\boldsymbol{y}) - \log\displaystyle\frac{}{}{p(\boldsymbol{y})}{q(\boldsymbol{y})}\right\} d\boldsymbol{y} = H(p) + D\!\left(p\|q\right).\qquad(71)

m層の第i成分から第m\!+\!1層の第j成分への経路i\!\to\!jの加重w_m^{ij}は、Eq. 72の勾配法で最適化できる。 係数\etaは学習率である。また、変数x_m^iは第m層の入力\boldsymbol{x}_mの第i成分でz_m^jW_m\boldsymbol{x}_mの第j成分である。

w_m'^{ij} = w_m^{ij}-\eta\displaystyle\frac{\partial }{\partial }{E}{w_m^{ij}} = w_m^{ij}-\eta \displaystyle\frac{\partial z_m^j}{\partial w_m^{ij}} \displaystyle\frac{\partial x_{m+1}^j}{\partial z_m^j} \displaystyle\frac{\partial E}{\partial x_{m+1}^j} = w_m^{ij}-\eta x_m^i \displaystyle\frac{\partial f}{\partial z_m^j}(z_m^j) \displaystyle\frac{\partial E}{\partial x_{m+1}^j}.\qquad(72)

各層で損失関数Eの導関数が必要だが、Eq. 73の漸化式を通じ、出力層から入力層に向けて逆伝播できる。

\displaystyle\frac{\partial E}{\partial x_m^i} = \displaystyle\sum_{j=1}^J \displaystyle\frac{\partial z_m^j}{\partial x_m^i} \displaystyle\frac{\partial x_{m+1}^j}{\partial z_m^j} \displaystyle\frac{\partial E}{\partial x_{m+1}^j} = \displaystyle\sum_{j=1}^J w_m^{ij} \displaystyle\frac{\partial f}{\partial z_m^j}(z_m^j) \displaystyle\frac{\partial E}{\partial x_{m+1}^j}.\qquad(73)

なお、活性化関数fはEq. 72Eq. 73の適用を助ける目的で、ReLU関数を除き微分可能な関数で設計される。

\displaystyle\frac{\partial f_\mathrm{sigm}}{\partial z_m^j}(z_m^j) = \displaystyle\frac{e^{-z_m^j}}{(1+e^{-z_m^j})^2} = x_{m+1}^j(1-x_{m+1}^j).\qquad(74)

損失関数Eの逆伝播による学習を誤差逆伝播法と呼ぶ。出力層の直前のE_\mathrm{sq}の導関数は、Eq. 75で求まる。

\displaystyle\frac{\partial E_\mathrm{sq}}{\partial x_3^j}(\boldsymbol{x}_3,\boldsymbol{t}) = \displaystyle\frac{\partial }{\partial x_3^j} \displaystyle\frac{1}{2} \|\boldsymbol{x}_3-\boldsymbol{t}\|^2 = x_3^j - t^j.\qquad(75)

以上の議論に基づき、ニューラルネットワークと誤差逆伝播法を実装する。まず、活性化関数fを定義する。

trait Active {
    def fp(z: Seq[Double]): Seq[Double]
    def bp(y: Seq[Double]): Seq[Double]
}

Activeを継承したSigmoidクラスは関数f_\mathrm{sigm}を表す。Eq. 68Eq. 74をfpとbpの各メソッドで実装する。

class Sigmoid extends Active {
    def fp(z: Seq[Double]) = z.map(z=>1/(1+Math.exp(-z)))
    def bp(z: Seq[Double]) = fp(z).map(y=>y*(1-y))
}

Fig.  14の入力層と隠れ層と出力層に相当し、Eq. 69Eq. 73の順伝播と逆伝播を定義するNeuralを宣言する。

trait Neural {
    val dim: Int
    def apply(x: Seq[Double]): Seq[Double]
    def apply(x: Seq[Double], t: Seq[Double]): Seq[Double]
}

2重定義のapplyメソッドで伝播を実行する。次に実装するOutputクラスは、出力層を具体的に定義する。

class Output(val dim: Int = 1, loss: (Double,Double)=>Double = _-_) extends Neural {
    def apply(x: Seq[Double]) = x
    def apply(x: Seq[Double], t: Seq[Double]) = (x,t).zipped.map(loss)
}

同様にNeuralを継承してNeuronクラスを実装する。Neuronクラスは、入力層ないし隠れ層の働きをする。

class Neuron(val dim: Int, act: Active, next: Neural, sgd: ()=>SGD) extends Neural {
    val W = Seq.fill(next.dim, dim)(sgd())
    def apply(x: Seq[Double]) = next(act.fp(z(x)))
    def apply(x: Seq[Double], t: Seq[Double]): Seq[Double] = {
        val xE = next(act.fp(z(x)), t)
        val zE = (xE, act.bp(z(x))).zipped.map(_*_)
        for((w,ze)<-W zip zE; (w,x)<-w zip x) w.update(x*ze)
        return W.transpose.map((_,zE).zipped.map(_.w*_).sum)
    }
    def z(x: Seq[Double]) = W.map((_,x).zipped.map(_.w*_).sum)
}

続くOffsetクラスはNeuronクラスを隠蔽し、層の入力\boldsymbol{x}_mに定数1を追加する。定数項の役割を果たす。

class Offset(val dim: Int, act: Active, next: Neural, sgd: ()=>SGD) extends Neural {
    val hidden = new Neuron(dim + 1, act, next, sgd)
    def offset = hidden.W.map(_.last.w)
    def apply(x: Seq[Double]) = hidden(x:+1d)
    def apply(x: Seq[Double], t: Seq[Double]) = hidden(x:+1d, t).init
}

最後にSGDクラスを定義する。これは、加重w_m^{ij}を保持すると同時に、勾配法によるw_m^{ij}の更新も実装する。

abstract class SGD(var w: Double = math.random) {
    def update(e: Double): Unit
}

下記のPlainSGDクラスは、Eq. 72に示した通りの基礎的な勾配法を実装する。学習率\etaは引数eに渡す。

class PlainSGD(e: Double = 0.01) extends SGD {
    def update(e: Double) = this.w -= e * this.e
}

誤差逆伝播法の実装は以上である。活性化関数f_\mathrm{sigm}と2乗誤差E_\mathrm{sq}の3層パーセプトロンの実装例を示す。

val model3 = new Output(1, _-_)
val model2 = new Offset(3, new Sigmoid, model3, ()=>new PlainSGD)
val model1 = new Offset(2, new Sigmoid, model2, ()=>new PlainSGD)

境界線が超平面では表現不可能な非線型分離の分類問題の例として、排他的論理和の学習を実験してみよう。

for(n<-1 to 1000000; x<-0 to 1; y<-0 to 1) model1(Seq(x,y), Seq(x^y))

Fig.  15 (a)に定数項なしの、 15 (b)に定数項を含む場合の学習結果を示す。定数項の有無で境界線が変化する。

a b

Figure 15: exclusive OR learned by a three-layer perceptron.. a — Neuron + Neuron., b — Offset + Offset.

同じ境界線でも、加重の初期値を変えて何度も試行すれば、収束に要する時間が変化する様子を観察できる。 余談だが、Fig.  13 (b)に示した論理和の境界線を再現するには、線型分類器を下記の通りに構築し学習する。

val model = new Offset(2, new Sigmoid, new Output, ()=>new PlainSGD)
for(n<-1 to 1000000; x<-0 to 1; y<-0 to 1) model(Seq(x,y), Seq(x|y))

5.2 鞍点と学習率

Eq. 76に与える関数Eに対し、原点\boldsymbol{O}は勾配が\nabla E\!=\!\boldsymbol{0}であるが、極小点ではない。これを鞍点と呼ぶ。

\Delta E = \displaystyle\frac{\partial f}{\partial x} \Delta x + \displaystyle\frac{\partial f}{\partial y} \Delta y = 2x \Delta x - 2y \Delta y\;\enspace\mathrm{where}\enspace E(x,y) = x^2 - y^2.\qquad(76)

Fig.  16 (a)に示す通り、鞍点の近傍ではEq. 72の勾配法は停滞し、運悪く鞍点に嵌れば、そこで収束する。 Fig.  16 (b)は、5個の初期値で排他的論理和を学習した際の誤差E_\mathrm{sq}の推移だが、学習の停滞が観察できる。

a b

Figure 16: comparison of PlainSGD and AdaDelta.. a — saddle point avoidance mechanism., b — loss diminution through training.

対策として、最適解の近傍で学習率を小さく、鞍点の近傍では大きく設定する適応的勾配法が効果的である。 Eq. 77に示すAdaGradは、学習率を勾配\nabla E_{mt}^{ij}の期待値に反比例させつつ、時刻tに伴って減衰させる。

\Delta w_{mt}^{ij} = -\displaystyle\frac{\eta}{t\sqrt{\underset{(\nabla E_m^{ij})^2}{\mathbf{E}}\!\left[\,\,\right]_t}}\; \enspace\mathrm{where}\enspace \underset{(\nabla E_m^{ij})^2}{\mathbf{E}}\!\left[\,\,\right]_t = \displaystyle\frac{1}{t} \displaystyle\sum_{\tau=0}^t (\nabla E_{mt}^{ij})^2,\; \underset{(\nabla E_m^{ij})^2}{\mathbf{E}}\!\left[\,\,\right]_0 = \varepsilon.\qquad(77)

Eq. 78に示すAdaDeltaは、直近の勾配を重視する。加重w_m^{ij}と勾配\nabla E_m^{ij}スケール変換の効果も持つ。

\Delta w_{mt}^{ij} = -\displaystyle\frac{\sqrt{\underset{(\Delta w_m^{ij})^2}{\mathbf{E}}\!\left[\,\,\right]_t+\varepsilon}}{\sqrt{\underset{(\nabla E_m^{ij})^2}{\mathbf{E}}\!\left[\,\,\right]_t+\varepsilon}} \nabla E_{mt}^{ij}\; \enspace\mathrm{where}\enspace \underset{x}{\mathbf{E}}\!\left[\,\,\right]_t = \rho \underset{x}{\mathbf{E}}\!\left[\,\,\right]_{t-1} + (1-\rho) x_t,\; \underset{x}{\mathbf{E}}\!\left[\,\,\right]_0 = 0.\qquad(78)

下記のAdaDeltaクラスに実装する。引数rは係数\rhoを表す。また、引数eはゼロ除算を防ぐ係数\varepsilonである。

class AdaDelta(r: Double = 0.95, e: Double = 1e-8) extends SGD {
    var eW,eE = 0.0
    def update(E: Double) = {
        this.eE = r*eE + (1-r) * math.pow(1*E, 2)
        val n = math.sqrt(eW+e) / math.sqrt(eE+e)
        this.eW = r*eW + (1-r) * math.pow(n*E, 2)
        this.w -= n * E
    }
}

PlainSGDクラスでFig.  16 (a)の軌跡を描くには419回の更新を要したが、AdaDeltaクラスは66回である。

5.3 多クラス分類

出力y\!\in\!\left\{1,...,K\right\}の分類器の活性化関数は、関数f_\mathrm{sigm}Kクラスに拡張したソフトマックス関数を使う。

y \sim q(y) = f_\mathrm{smax}(\boldsymbol{z}) = \displaystyle\frac{f_\mathrm{exp}(\boldsymbol{z})}{\|f_\mathrm{exp}(\boldsymbol{z})\|}\; \enspace\mathrm{where}\enspace f_\mathrm{exp}(\boldsymbol{z}) = {}^t\!e^{z_1} \cdots e^{z_K}\ .\qquad(79)

直感的には、最適な事後分布q(y)は正解p(t)との差を表すKullback–Leibler 情報量D\!\left(p\|q\right)を最小化する。 Kullback–Leibler 情報量は非負で、Eq. 80に示すギブスの不等式よりq\!=\!pに限りD\!\left(p\|q\right)\!=\!0が成立する。

D\!\left(p\|q\right) = \int_K p(y) \log \displaystyle\frac{p(y)}{q(y)} dy \geq \int_K p(y) \left(1-\displaystyle\frac{q(y)}{p(y)}\right) dy = 0\;\leftarrow \log x \leq x-1.\qquad(80)

通常はD\!\left(p\|q\right)の上限を与えるEq. 71の交差エントロピーE_\mathrm{CE}を通じて、間接的にD\!\left(p\|q\right)を最小化する。 理由は、出力層の直前の勾配\nablaE_\mathrm{CE}(\boldsymbol{z})がEq. 81で容易に求まり、Eq. 71でH(p)が定数である点にある。

\displaystyle\frac{\partial E_\mathrm{CE}}{\partial z_k} = - \displaystyle\frac{\partial }{\partial z_k} \displaystyle\sum_{i=1}^K t_i \left(\log e^{z_i} - \log \displaystyle\sum_{j=1}^K e^{z_j}\right) = - t_k + \displaystyle\sum_{i=1}^K t_i y_k = -t_k + y_k.\qquad(81)

下記のSoftmaxクラスにEq. 79を実装する。出力層の直前での利用を前提に、逆伝播は恒等関数にした。

class Softmax extends Active {
    def fp(z: Seq[Double]) = z.map(math.exp(_)/z.map(math.exp).sum)
    def bp(z: Seq[Double]) = z.map(_=>1.0)
}

Softmaxクラスは、出力層の直前での利用に限定される点を除き、他の活性化関数と同じ要領で利用できる。

val model3 = new Output(4, _-_)
val model2 = new Offset(3, new Softmax, model3, ()=>new PlainSGD)
val model1 = new Offset(2, new Sigmoid, model2, ()=>new PlainSGD)

Fig.  17は、xy軸上の4点を標本に与えて、3層パーセプトロンで国際信号旗のzuluを学習した結果である。

Figure 17: maritime signal flag zulu learned by a three-layer perceptron.

原理上は、第 5.1節と同様に2乗誤差E_\mathrm{sq}で学習できるが、勾配\nablaE_\mathrm{sq}が緩慢なため、収束に時間がかかる。

6 サポートベクターマシン

サポートベクターマシンとは、Fig.  18に示す判別境界を学習し、標本\left\{\boldsymbol{x}_i\right\}を正負に分割する分類器である。 両側の集団からの最短距離dが最大になる直線を学習することで、未知の\boldsymbol{x}を誤分類する可能性を抑制する。

Figure 18: a support vector machine.

判別境界\boldsymbol{w}\! \cdot \!\boldsymbol{x}\!+\!c\!=\!0の学習は制約付き最大化問題であり、標本\left\{\boldsymbol{x}_i,t_i\right\}は条件Eq. 82を満たす必要がある。

t_i(\boldsymbol{w} \cdot \boldsymbol{x}_i + c) \geq 1\; \enspace\mathrm{where}\enspace t_i= \begin{cases} +1 & \text{if $\boldsymbol{x}_i\in\mathrm{positive\;group}$} \\ -1 & \text{if $\boldsymbol{x}_i\in\mathrm{negative\;group}$} \end{cases}\qquad(82)

ただし、t_iは点\boldsymbol{x}_iが帰属する集団を示す変数である。集団と判別境界との距離dは、Eq. 83で計算できる。

d(\left\{\boldsymbol{x}_i\right\}) = \min\displaystyle\frac{}{}{|\boldsymbol{w} \cdot \boldsymbol{x}_i + c|}{\|\boldsymbol{w}\|} = \displaystyle\frac{1}{\|\boldsymbol{w}\|}.\qquad(83)

現実には、正負の境界が曖昧な場合には、Eq. 82の遵守が困難なため、ソフトマージンによる緩和を図る。 即ち、制約条件をEq. 84に変更し、若干の誤分類を見逃す代わりに、誤分類点\boldsymbol{x}_iにヒンジ損失\xi_iを課す。

t_i (\boldsymbol{w} \cdot \boldsymbol{x}_i + c) \geq 1 - \xi_i\; \enspace\mathrm{where}\enspace \xi_i = \begin{cases} 0 & \text{if $t_i(\boldsymbol{w} \cdot \boldsymbol{x}_i + c) > 1$} \\ |t_i - (\boldsymbol{w} \cdot \boldsymbol{x}_i + c)| & \text{if $t_i(\boldsymbol{w} \cdot \boldsymbol{x}_i + c) \leq 1$}. \end{cases}\qquad(84)

誤分類の抑制と距離dの最大化は、Eq. 85の最小化で達成する。ここで、\|\boldsymbol{w}\|^2はL2正則化の役割を担う。

f(\left\{\boldsymbol{x}_i,t_i\right\}, \boldsymbol{w}, c) = C \displaystyle\sum_{i=1}^N \xi_i + \displaystyle\frac{1}{2}\|\boldsymbol{w}\|^2 \enspace\mathrm{where}\enspace C>0.\qquad(85)

誤分類点の集合E(\left\{\boldsymbol{x}_i\right\})に対し、Eq. 86が成立する。故に、定数Cは誤分類率の限度を間接的に決定する。

\displaystyle\sum_{i=1}^N \xi_i > |E(\left\{\boldsymbol{x}_i\right\})| \leftarrow \forall\boldsymbol{x}_i \in E(\left\{\boldsymbol{x}_i\right\}),\;\xi_i>1.\qquad(86)

6.1 凸二次計画問題

Eq. 85の最小化は、未定乗数法とKarush–Kuhn–Tucker 条件の適用により、凸二次計画問題で達成できる。

L(\boldsymbol{w},c,\xi,\lambda,\mu) = \displaystyle\frac{1}{2}\|\boldsymbol{w}\|^2 +C \displaystyle\sum_{i=1}^N \xi_i - \displaystyle\sum_{i=1}^N \lambda_i \{t_i(\boldsymbol{w} \cdot \boldsymbol{x}_i + c) - 1\} - \displaystyle\sum_{i=1}^N \mu_i \xi_i.\qquad(87)

\lambda_i\mu_iは座標\boldsymbol{x}_iに関する未定乗数である。係数\boldsymbol{w}cで偏微分して、Lが最小になる条件Eq. 88を得る。

\tilde{L}(\lambda) = \min_{\boldsymbol{w},c} L(\boldsymbol{w},c, \lambda) = \displaystyle\sum_{i=1}^N \lambda_i \left\{1 - \displaystyle\displaystyle\frac{}{}{1}{2} \displaystyle\displaystyle\sum_{j=1}^N \lambda_j t_i t_j (\boldsymbol{x}_i \cdot \boldsymbol{x}_j)\right\} \enspace\mathrm{where}\enspace\left\{ \begin{aligned} 0 &= \displaystyle\sum_{i=1}^N \lambda_i t_i \\ \lambda_i &= C - \mu_i\\ \boldsymbol{w} &= \displaystyle\sum_{i=1}^N \lambda_i t_i \boldsymbol{x}_i \end{aligned} \right..\qquad(88)

Eq. 84は不等式ゆえ、Karush–Kuhn–Tucker 条件より、Eq. 88が最適解を与える場合、Eq. 89を満たす。 誤分類点の存在を無視すると、Eq. 88Eq. 89より、判別境界から距離dの点のみが係数\boldsymbol{w}の決定に寄与する。

\lambda_i \left\{t_i\left(\displaystyle\displaystyle\sum_{j=1}^N \lambda_j t_j \boldsymbol{x}_j \cdot \boldsymbol{x}_i + c\right) + \xi_i - 1\right\} = 0\; \enspace\mathrm{where}\enspace\left\{ \begin{aligned} \lambda_i &\geq 0 \\ \mu_i &\geq 0 \\ \mu_i\xi_i&=0 \end{aligned} \right..\qquad(89)

これをサポートベクトルと呼ぶ。より遠方の点は、\lambda_i\!=\!0を満たす必要から、係数\boldsymbol{w}の決定に無関係である。 以後、双対問題\tilde{L}(\lambda)を最大化する。今回は、実装が容易で計算が高速な逐次最小問題最適化法を利用する。

t_i \delta_i + t_j \delta_j = 0\; \enspace\mathrm{where}\enspace \left\{ \begin{aligned} \delta_i &= \hat{\lambda}_i - \lambda_i \\ \delta_j &= \hat{\lambda}_j - \lambda_j \end{aligned} \right\}\; \leftarrow 0 = \displaystyle\sum_{i=1}^N \lambda_i t_i.\qquad(90)

これは、Eq. 89を破る点\boldsymbol{x}_iが存在する限り、任意の点\boldsymbol{x}_jを選び、Eq. 90を満たす局所的な最適化を施す。 Eq. 88の\tilde{L}(\lambda)から、1回の最適化による\lambda_i,\lambda_jの増分\delta_i,\delta_jを含む部分式を抜き出すと、Eq. 91を得る。

d\tilde{L}(\delta_i, \delta_j) = \delta_i +\delta_j - \displaystyle\frac{1}{2} |\delta_i t_i \boldsymbol{x}_i + \delta_j t_j \boldsymbol{x}_j|^2 - \displaystyle\sum_{n=1}^N \lambda_n t_n \boldsymbol{x}_n \cdot (\delta_i t_i \boldsymbol{x}_i + \delta_j t_j \boldsymbol{x}_j).\qquad(91)

Eq. 91にEq. 90の制約を加えると、Eq. 92を得る。なお、Eq. 92は変数t_iの符号に拘らず成立する。

d\tilde{L}(\delta_i) = t_i (t_i - t_j) \delta_i - \displaystyle\frac{1}{2} \delta_i^2 |\boldsymbol{x}_i - \boldsymbol{x}_j|^2 - t_i \delta_i \displaystyle\sum_{n=1}^N \lambda_n t_n \boldsymbol{x}_n \cdot (\boldsymbol{x}_i - \boldsymbol{x}_j).\qquad(92)

双対問題\tilde{L}(\lambda)の最大化は、d\tilde{L}(\delta_i,\delta_j)の極大点の探索と同義である。\delta_iでの偏微分により、Eq. 93を得る。 なお、更新後の値\hat{\lambda}_i\hat{\lambda}_jがEq. 89を逸脱する場合は、クリッピングにより強制的に区間[0,C]に収める。

\delta_i = - \displaystyle\frac{t_i}{|\boldsymbol{x}_i - \boldsymbol{x}_j|^2} \left\{\displaystyle\displaystyle\sum_{n=1}^N \lambda_n t_n \boldsymbol{x}_n \cdot (\boldsymbol{x}_i - \boldsymbol{x}_j) - t_i + t_j\right\}.\qquad(93)

座標\boldsymbol{x}_iに対し、条件Eq. 89を確認する際は、係数\boldsymbol{w}cの値が必要になるが、\boldsymbol{w}の値はEq. 88で求まる。 定数項cの値は、t_i(\boldsymbol{w}\! \cdot \!\boldsymbol{x}_i)を最小化する座標\boldsymbol{x}_iがサポートベクトルとなる点に着目し、Eq. 94で求める。

c = -\displaystyle\frac{1}{2} \left\{ \min_{i|t_i=+1} \displaystyle\displaystyle\sum_{j=1}^N \lambda_j t_j \boldsymbol{x}_i \cdot \boldsymbol{x}_j + \max_{j|t_j=-1} \displaystyle\displaystyle\sum_{i=1}^N \lambda_i t_i \boldsymbol{x}_j \cdot \boldsymbol{x}_j\right\}.\qquad(94)

下記のSVMクラスはサポートベクターマシンを実装する。標本\left\{\boldsymbol{x}_i\right\}は引数Xに渡す。引数kは内積を表す。

class SVM(X: Seq[(Seq[Double],Int)], C: Double, k: (Seq[Double],Seq[Double])=>Double) {
    val data = X.map(Data.tupled)
    var c = 0.0

SVMクラス内にDataクラスを実装する。これは、標本のi番目の点\boldsymbol{x}_iと正解t_iと未定乗数\lambda_iを記憶する。

    case class Data(x: Seq[Double], t: Int) {
        var l = 0.0
    }

判別境界を表す係数\boldsymbol{w}は変数に記憶せず、内積\boldsymbol{w}\! \cdot \!\boldsymbol{x}の形で実装する。これは第 6.2節に向けた布石である。

    def wx(x: Data) = data.map(d => d.l * d.t * k(x.x,d.x)).sum

kktメソッドは、Eq. 89のKarush–Kuhn–Tucker 条件を検査する。実数の比較では丸め誤差に注意を払う。

    def kkt(x: Data) = (x.t * this(x.x) - 1) match {
        case e if e < -1e-10 => x.l >= C
        case e if e > +1e-10 => x.l <= 0
        case _ => true
    }

clipメソッドは、Eq. 93の最適化を行う際に、未定定数の差分\delta_i\delta_jの閉区間[0,C]からの逸脱を防ぐ。

    def clip(xi: Data, xj: Data, d: Double): Double = {
        if(d < lower(xi, xj)) return lower(xi, xj)
        if(d > upper(xi, xj)) return upper(xi, xj)
        return if(d.isNaN) 0 else d
    }

lowerメソッドは、差分\delta_i\delta_jの下限を計算する。下限はEq. 90の\delta_iの定義とEq. 88Eq. 89から導ける。

    def lower(xi: Data, xj: Data) = xi.t * xj.t match {
        case  1 => math.max(-xi.l, +xj.l - C)
        case -1 => math.max(-xi.l, -xj.l)
    }

同様に、差分\delta_i\delta_jの上限もEq. 90の\delta_iの定義とEq. 88Eq. 89から導出し、upperメソッドに実装する。

    def upper(xi: Data, xj: Data) = xi.t * xj.t match {
        case  1 => math.min(+xj.l, -xi.l + C)
        case -1 => math.min(-xi.l, -xj.l) + C
    }

次に実装するposとnegは、正負のサポートベクトルの内積\boldsymbol{w}\! \cdot \!\boldsymbol{x}を求める。Eq. 94を計算する際に使う。

    def pos = data.filter(_.t == +1).map(wx(_)).min
    def neg = data.filter(_.t == -1).map(wx(_)).max

下記は、逐次最小問題最適化法の反復処理である。Eq. 89を破る点\boldsymbol{x}_iを探し、Eq. 93に従って更新する。

    while(data.count(!kkt(_)) >= 2) {
        val a = data.find(!kkt(_)).get
        val b = data(util.Random.nextInt(X.size))
        val sub = wx(Data((a.x,b.x).zipped.map(_-_), 0))
        val den = k(a.x,a.x) - 2*k(a.x,b.x) + k(b.x,b.x)
        val del = clip(a, b, -a.t * (sub-a.t+b.t) / den)
        a.l += del
        b.l -= del * a.t * b.t
        this.c = -0.5 * (pos+neg)
    }

下記のapplyメソッドは、未知の点\boldsymbol{x}に対して\boldsymbol{w}\! \cdot \!\boldsymbol{x}\!+\!cを計算する。返り値の正負は、帰属する集団を表す。

    def apply(x: Seq[Double]) = wx(Data(x, 0)) + c
}

完成したSVMクラスを使う際は、標本\boldsymbol{x}_iに加え、ソフトマージンの定数Cと内積の定義を引数に指定する。

val svm = new SVM(data, 1e-10, (a, b) => (a, b).zipped.map(_*_).sum)

Fig.  19は、線型分離可能な標本をSVMクラスで学習した結果である。ただし 19 (b)は2個の誤分類点を含む。 黒の点線は判別境界を表し、赤と青の点線は\boldsymbol{w}\! \cdot \!\boldsymbol{x}\!+\!c\!=\!\pm1を表す。赤と青の濃淡は\boldsymbol{w}\! \cdot \!\boldsymbol{x}\!+\!cの相対値を表す。

a b

Figure 19: decision surface learned by a linear SVM.. a — sample data separable by a line., b — sample data with outlier points.

 19 (a)は、両集団を直線で完璧に分離できる場合、判別境界がサポートベクトルの中間に位置する様子を表す。  19 (b)は、判別境界が正の領域にある負の誤分類点に誘引され、かつサポートベクトルが消滅する様子を表す。

6.2 特徴空間の変換

第 6.2節では、サポートベクターマシン等の線型分類器を非線型分離な問題に対応させるカーネル法を学ぶ。 根底には、標本を高次元空間に写像することで、像の分布を疎にした上で、線型分離を達成する着想がある。

\Phi: \boldsymbol{x} \mapsto \Phi_{\boldsymbol{x}}\;\; \mathrm{e.g.}\; \begin{pmatrix}x\\y\end{pmatrix} \mapsto z=\displaystyle\frac{1}{2\pi} \exp \left(-\displaystyle\frac{x^2 + y^2}{2}\right).\qquad(95)

Fig.  20は、同心円に並ぶ正負の集団をEq. 95でz軸上に投影し、z軸に垂直な分離平面を得る例である。

Figure 20: conversion of linearly inseparable problem to separable.

Eq. 82–Eq. 94の座標\boldsymbol{x}\Phi_{\boldsymbol{x}}に変更し、第 6.2節の議論に流用する。例えば、係数\boldsymbol{w}はEq. 96で求める。 基本的に、像\Phi_{\boldsymbol{x}}の次元Dが高ければ分離の機会は増えるが、Eq. 82–Eq. 94の計算は\mathcal{O}(D)の時間を要する。

\boldsymbol{w} = \displaystyle\sum_{i=1}^N \lambda_i t_i \Phi_{\boldsymbol{x}_{i}}.\qquad(96)

ここで、Eq. 89Eq. 93Eq. 94を注意深く観察すると、内積\Phi_{\boldsymbol{x}_{i}}\! \cdot \!\Phi_{\boldsymbol{x}_{j}}が求まれば\Phi_{\boldsymbol{x}}の値は不要な点に気付く。 実は、像の計算を経ず内積を計算する手段が存在する。まず、Eq. 97に示す2変数の対称関数を定義する。

k: \mathbb{M} \times \mathbb{M} \to \mathbb{R}\;\enspace\mathrm{where}\enspace k(\boldsymbol{x}_i,\boldsymbol{x}_j) = k(\boldsymbol{x}_j,\boldsymbol{x}_i).\qquad(97)

記号\mathbb{M}可測空間を指す。関数kは、Eq. 98に示す正定値性を満たす場合、正定値カーネルと呼ばれる。

\forall\boldsymbol{x}\in\mathbb{R}^n,\; {}^t\!\boldsymbol{x}\ \begin{pmatrix} k(\boldsymbol{x}_1,\boldsymbol{x}_1) & \cdots & k(\boldsymbol{x}_1,\boldsymbol{x}_n) \\ \vdots & \ddots & \vdots \\ k(\boldsymbol{x}_n,\boldsymbol{x}_1) & \cdots & k(\boldsymbol{x}_n,\boldsymbol{x}_n) \\ \end{pmatrix} \boldsymbol{x} > 0.\qquad(98)

Eq. 98の関数kに対し、Eq. 99の像\Phi_{\boldsymbol{x}_{i}}の線型結合からなる空間H_k再生核ヒルベルト空間と呼ぶ。

H_k = \left\{f(\boldsymbol{x}) = \displaystyle\displaystyle\sum_{n=1}^N a_n \Phi_{\boldsymbol{x}_{n}} = \displaystyle\displaystyle\sum_{n=1}^N a_n k(\boldsymbol{x}, \boldsymbol{x}_n)\right\} \enspace\mathrm{where}\enspace \Phi: \boldsymbol{x} \mapsto \Phi_{\boldsymbol{x}_{n}} = k(\boldsymbol{x}, \boldsymbol{x}_n).\qquad(99)

空間H_kの元fgに対し、内積\left\langle f,g\right\rangleは非退化かつ正定値な対称双線型形式であり、Eq. 100で定義する。

\left\langle f,g\right\rangle = \displaystyle\sum_{i=1}^N \displaystyle\sum_{j=1}^N a_i b_j k(\boldsymbol{x}_i, \boldsymbol{x}_j)\; \enspace\mathrm{where}\enspace\left\{ \begin{aligned} f(\boldsymbol{x}) &= \displaystyle\sum_{i=1}^N a_i k(\boldsymbol{x}, \boldsymbol{x}_i) \\ g(\boldsymbol{x}) &= \displaystyle\sum_{j=1}^N b_j k(\boldsymbol{x}, \boldsymbol{x}_j) \\ \end{aligned} \right\}\; a_i,b_j \in \mathbb{R}.\qquad(100)

余談ながら、内積と距離を備え、かつ完備な空間をヒルベルト空間と呼び、Eq. 101の性質を再生性と呼ぶ。

\left\langle f,k(\boldsymbol{x},\boldsymbol{x}_j)\right\rangle = \displaystyle\sum_{i=1}^N a_i k(\boldsymbol{x}_i, \boldsymbol{x}_j) = f(\boldsymbol{x}_j).\qquad(101)

Eq. 101の性質は、像\Phi_{\boldsymbol{x}_{i}}\Phi_{\boldsymbol{x}_{j}}の明示的な計算を経ずに、内積\Phi_{\boldsymbol{x}_{i}}\! \cdot \!\Phi_{\boldsymbol{x}_{j}}の値が求まる可能性を示唆する。 この技法はカーネルトリックと通称される。関数kの定義次第では、無限次元の空間への変換も可能である。

\Phi_{\boldsymbol{x}_{i}} \cdot \Phi_{\boldsymbol{x}_{j}} = k(\boldsymbol{x}_i,\boldsymbol{x}_j) = \left\{ \begin{aligned} & \displaystyle\sum_{d=1}^D x_{id} x_{jd} && \text{$\cdots$ linear kernel} \\ & \exp \left(-\displaystyle\frac{|\boldsymbol{x}_i - \boldsymbol{x}_j|^2}{2\sigma^2}\right) && \text{$\cdots$ Gaussian kernel} \end{aligned} \right.\qquad(102)

Eq. 102のうち線型カーネルは単なる恒等変換に過ぎないが、ガウシアンカーネルには興味深い特性がある。

\exp \displaystyle\frac{\boldsymbol{x}_i \cdot \boldsymbol{x}_j}{\sigma^2} = \exp \displaystyle\sum_{d=1}^D \displaystyle\frac{x_{id} x_{jd}}{\sigma^2} = \displaystyle\prod_{d=1}^D \exp \displaystyle\frac{x_{id} x_{jd}}{\sigma^2} = \displaystyle\prod_{d=1}^D \displaystyle\sum_{n=0}^\infty \displaystyle\frac{1}{\sqrt{n!}} \left(\displaystyle\frac{x_{id}}{\sigma}\right)^n \displaystyle\frac{1}{\sqrt{n!}} \left(\displaystyle\frac{x_{jd}}{\sigma}\right)^n.\qquad(103)

Eq. 103は、Eq. 104のような無限次元のベクトルの内積と同義であり、事実上、無限次元への写像を表す。

\displaystyle\sum_{n=0}^\infty \displaystyle\frac{1}{\sqrt{n!}} \left(\displaystyle\frac{x_{id}}{\sigma}\right)^n \displaystyle\frac{1}{\sqrt{n!}} \left(\displaystyle\frac{x_{jd}}{\sigma}\right)^n = \begin{pmatrix} \displaystyle\frac{1}{\sqrt{0!}} \left(\displaystyle\frac{x_{id}}{\sigma}\right)^0 \\ \vdots \\ \displaystyle\frac{1}{\sqrt{n!}} \left(\displaystyle\frac{x_{id}}{\sigma}\right)^n \\ \vdots \end{pmatrix} \cdot \begin{pmatrix} \displaystyle\frac{1}{\sqrt{0!}} \left(\displaystyle\frac{x_{jd}}{\sigma}\right)^0 \\ \vdots \\ \displaystyle\frac{1}{\sqrt{n!}} \left(\displaystyle\frac{x_{jd}}{\sigma}\right)^n \\ \vdots \end{pmatrix}\qquad(104)

他の著名なカーネルの例では、Eq. 105のシグモイドカーネルは3層パーセプトロンに似た挙動で知られる。

\mathrm{sign}\left(\boldsymbol{w} \cdot \Phi_{\boldsymbol{x}} + c\right) = \mathrm{sign}\left(\displaystyle\sum_{i=1}^N \lambda_it_i k_s(\boldsymbol{x},\boldsymbol{x}_i) + c\right)\; \enspace\mathrm{where}\enspace k(\boldsymbol{x}_i, \boldsymbol{x}_j) = \tanh \left(c \displaystyle\sum_{d=1}^D x_{id} x_{jd} + \theta\right).\qquad(105)

Fig.  21は、SVMクラスの引数\texttt{k}にガウシアンカーネルを与えて、非線型分離な標本を学習した結果である。

a b

Figure 21: decision surface learned by a Gaussian SVM.. a — diamond-shaped samples., b — two concentric circles.

黒の点線は判別境界を表し、赤と青の濃淡は\boldsymbol{w}\! \cdot \!\boldsymbol{x}\!+\!cの相対値を表す。カーネル法の利点がよく理解できる。