Scala で実装する機械学習入門
Let’s implement machine learning by Scala
無線部開発班
平成 28 年 2 月 22 日
目次
第 1 章 最近傍探索
2
第 2 章 ニューラルネットワーク
2.1
単純パーセプトロン
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
4
2.2
多層パーセプトロン
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
第 3 章 サポートベクターマシン
3.1
3.2
10
ハードマージンの場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ソフトマージンの場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
14
第 4 章 カーネルトリック
4.1
4.2
4.3
4.4
17
高次元への写像 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
ヒルベルト空間 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
再生核カーネル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
分類器への適用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
第 5 章 単純ベイズ分類器
18
19
21
23
第 6 章 最尤推定法
6.1
6.2
6.3
26
混合分布モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
最尤推定の理論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
最尤推定の実装 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
26
28
第 1 章 最近傍探索
k 近傍法は、分類対象のデータの近傍にある k 個のデータに多数決を取ってもらい、そのデータのクラスを
分類する手法である。Fig. 1.1 は 6 個の赤い正のサンプル点と 6 個の青い負のサンプル点に対し k 近傍法を
2.0
1.5
1.5
1.0
1.0
0.5
0.5
0
0.0
y
2.0
−0.5
−1.0
−1.0
0.5
0
0.0
0.5
y
行い、空間を正と負の領域に区分する例である。k = 2 とした場合と k = 3 とした場合で挙動が異なる。
−0.5
−0.5
0.0
0.5
x
1.0
1.5
2.0
−1.0
−1.0
−0.5
(a) k = 2
0.0
0.5
x
1.0
1.5
2.0
(b) k = 3
Fig. 1.1: k nearest neighbors.
実は、k 近傍法を使うにあたり、この k の適切な値の判断基準が問題になるのだが、本稿の議論の対象外と
する。また、距離関数をどのように定義するかも重要な問題である。典型的な例では、ユークリッド距離や
マンハッタン距離があるが、距離の公理を満たす関数ならば仔細を問わない。k 近傍法は、分類器としては
単純なアルゴリズムである。事前の学習が不要で、データを分類する際に初めて訓練データを参照する点が
特徴である。このような挙動から、遅延学習とも呼ばれる。では、実装してみよう。
knn.scala
case class Data(point: Seq[Double], answer: Any = 0)
Data クラスは 1 件のサンプルを表現するケースクラスで、空間内の座標と正負の正解ラベルをフィールドに
持つ。例えば、Fig. 1.1 の訓練データは下記のように定義した。ただし、1 が正のラベルで 0 が負のラベルで
ある。matplotlib で描画する都合からラベルを数値としたが、Any 以下の如何なるクラスでも構わない。
knn.scala
val samples = Seq(
Data(Seq(-0.10, +0.50),
Data(Seq(+0.60, +0.75),
Data(Seq(+0.75, +0.70),
Data(Seq(+0.20, +0.30),
Data(Seq(+0.40, +0.55),
Data(Seq(+0.80, +0.55),
1),
1),
1),
0),
0),
0),
Data(Seq(-0.45,
Data(Seq(+0.30,
Data(Seq(+1.20,
Data(Seq(+1.60,
Data(Seq(+0.60,
Data(Seq(+1.25,
2
+1.30),
+1.25),
+0.55),
+0.60),
+0.40),
+0.20),
1),
1),
1),
0),
0),
0))
無線部
Scala で実装する機械学習入門
開発班
学習用の訓練データも準備できたので、さっそく k 近傍法の本体を実装してみよう。
knn.scala
class KNN(k: Int, samples: Seq[Data], distance: (Data, Data) => Double) {
def predict(target: Data): Any = {
val knears = samples.sortWith((t1, t2) => {
val d1 = distance(t1, target)
val d2 = distance(t2, target)
d1 < d2
}).take(k).map(_.answer)
val answers = knears.toSet
val verdict = answers.map(ans => ans -> knears.count(_ == ans))
verdict.maxBy(tuple => tuple._2)._1
}
}
KNN クラスが k 近傍法による分類器である。以下のように k の値と訓練データと距離関数を与えて使用する。
knn.scala
val euclid = (a: Data, b: Data) => Math.sqrt(
Math.pow(a.point(0) - b.point(0), 2) +
Math.pow(a.point(1) - b.point(1), 2)
)
val knn = new KNN(3, samples , euclid) // k = 3
println(knn.predict(Data(Seq(0.0, 1.0))))
println(knn.predict(Data(Seq(1.0, 0.0))))
分類対象のデータを引数に predict メソッドを呼び出すと、訓練データとの距離を計算し、最近傍の k 個の
訓練データを取り出してラベルの多数決を行う。実際には、距離に応じた発言力の重み付けが必要な場合も
ある。Fig. 1.2 に示すように、少数派の集団が多数派の集団により囲まれている場合に、中心部を少数派に
分類したいからである。
2.0
1.5
1.0
0.0
0.50
y
0.5
−0.5
−1.0
−1.5
−2.0
−2.0 −1.5 −1.0 −0.5
0.0
x
0.5
1.0
1.5
2.0
Fig. 1.2: concentric circular neighbors.
また、分類のたびに全ての訓練データとの距離を求めるのは効率の点で不利である。単純だが用途に注意を
要する分類器と言えよう。なお、今回は平面上の点を分類したが、対象のデータは画像や文字列など様々な
形態が考えられる。これらを扱う際は、何らかの方法で距離空間への写像を行うことになる。写像後の値が
定義されるベクトル空間を特徴空間ないし素性空間と呼ぶ。
3
第 2 章 ニューラルネットワーク
ニューラルネットワークは脳神経の仕組みを模倣した数学モデルであり、最も歴史の長い機械学習モデルの
ひとつである。今回は、信号が出力方向にのみ伝搬するフィードフォワード型のニューラルネットワークを
実装する。教科書を読むと難しく思えるが、非常に単純なモデルである。手を動かして感覚を掴もう。
2.1
単純パーセプトロン
単純パーセプトロンとは、単体の神経細胞を模した、最も単純なニューラルネットワークである。D 次元の
入力 x と重み w との内積に、式 (2.1) のように活性化関数 f を適用するだけである。Fig. 2.1 に図示する。
(
y = f (w · x) = f
w0 +
D
∑
)
wd xd
.
(2.1)
d=1
Fig. 2.1: single layer perceptron.
活性化関数は、実数値の入力を受け取り、閾値以上なら正、未満なら負に分類する関数である。内積 w · x は
実数であり、活性化関数で真偽値に変換する。グラフを描けば瞭然だが、パーセプトロンとは、座標空間を
正と負の領域に分離する決定境界を用いて、データを弁別する分類器である。例えば、論理和を学習すると
Fig. 2.2 の直線を決定境界として習得し、直線より上側は正、下側は負に分類する。
Fig. 2.2: decision surface of OR operator.
4
無線部
Scala で実装する機械学習入門
開発班
直線上の点の分類は真でも負でも構わないが、本稿では真とする。決定境界を表す w は学習により求める。
学習の仕組みは後述するとして、分類器を実装する。Perceptron クラスが、パーセプトロンの本体である。
perceptron.scala
class Perceptron(samples: Seq[(Seq[Double], Boolean)], dim: Int) {
val weights = new Array[Double](dim + 1)
def predict(in: Seq[Double]) = ((1.0 +: in) zip weights).map {
case (i, w) => i * w
}.sum >= 0
}
weights は重み w に相当する。dim は入力 x の階数である。w の階数を dim+1 に設定するのは、式 (2.1) の
定数項 w0 の分である。predict メソッドは、式 (2.1) を実装したものである。以上で、分類器が完成した。
学習の仕組みを考えてみよう。n 番目のデータ xn の正解を tn = ± 1 とし、出力誤差を式 (2.2) で定義する。
t w · x
(youtput ̸= yanswer ),
n
n
En =
(2.2)
0
(y
=y
).
output
answer
学習の目的は誤差 En の総和を最小化することなので、勾配法、特に最急降下法を用いる。
w′ = w − η
N
∂ ∑
En .
∂w n=1
(2.3)
各次元の独立を仮定して d 次元目に着目すると、式 (2.3) は式 (2.4) に変形できる。
wd′
N
N
∑
∂ ∑
= wd − η
En = w d − η
tn xnd .
∂wd n=1
(2.4)
n;En ̸=0
式 (2.4) の計算を全次元で繰り返し実行することにより、徐々に重み w が改善される。では、実装してみよう。
perceptron.scala
val eta = 0.1
for(step <- 1 to 10000) samples.foreach{case (in, out) => predict(in) match {
case true if !out =>
for(k <- 0 until dim) weights(k + 1) -= eta * in(k)
weights(0) -= eta * 1
case false if out =>
for(k <- 0 until dim) weights(k + 1) += eta * in(k)
weights(0) += eta * 1
case _ =>
}}
以上を、先ほどの Perceptron クラスのコンストラクタとして追記すればよい。論理和を学習させてみよう。
perceptron.scala
val samples = Seq(
(Seq(0.0, 0.0), false),
(Seq(0.0, 1.0), true),
(Seq(1.0, 0.0), true),
(Seq(1.0, 1.0), true))
val p = new Perceptron(samples , 2);
println(p.predict(Seq(0.0, 0.0))) //
println(p.predict(Seq(0.0, 1.0))) //
println(p.predict(Seq(1.0, 0.0))) //
println(p.predict(Seq(1.0, 1.0))) //
false
true
true
true
5
無線部
Scala で実装する機械学習入門
開発班
本当に Fig. 2.2 に示した通りの直線を学習しているかを確かめるには、重み w の内容を出力してみれば良い。
perceptron.scala
println(p.weights.mkString(",")) // -0.1,0.1,0.1
論理和と同じ要領で、論理積を学習させることも可能である。しかし、排他的論理和の学習は必ず失敗する。
perceptron.scala
val samples = Seq(
(Seq(0.0, 0.0), false),
(Seq(0.0, 1.0), true),
(Seq(1.0, 0.0), true),
(Seq(1.0, 1.0), false))
val p = new Perceptron(samples , 2);
println(p.predict(Seq(0.0, 1.0))) //
println(p.predict(Seq(0.0, 1.0))) //
println(p.predict(Seq(0.0, 1.0))) //
println(p.predict(Seq(0.0, 1.0))) //
true
true
false
false
重み w は 0.0,-0.1,0.0 となり、y 軸そのものを学習しているが、これは至極当然の話で、排他的論理和は
直線では分離できない。決定境界が直線ないし超平面になる問題を線型分離可能と呼ぶ。論理和や論理積は
線型分離可能であるが、排他的論理和は線型分離不可能であり、単純パーセプトロンでは習得できない。
2.2
多層パーセプトロン
単純パーセプトロンは線型分離不可能な問題を扱うことができないが、パーセプトロンを何層も積み重ねる
ことで、任意の連続関数を近似できることが知られている。Fig. 2.3 に示す多層パーセプトロンである。
Fig. 2.3: multi layer perceptron.
パーセプトロンをいくつか並列に束ねることで出力 y をベクトルとし、そうしたパーセプトロンを何層にも
積み重ねることで多層パーセプトロンを構成する。左端で入力ベクトル x1 を受け取る層を入力層、右端で
出力ベクトル yN を出力する層を出力層と呼び、入力層と出力層の中間の層を隠れ層と呼ぶ。活性化関数を
f とすると、N 層パーセプトロンの最終段の出力 yN は式 (2.5) で計算される。
yN = f (WN yN −1 ) = f (WN f (WN −1 f (WN −2 ...f (W2 f (W1 x))...)).
(2.5)
ただし、便宜的に第 n 段の入力を xn 、出力を yn と表記する。従って、式 (2.6) が成立する。
xn+1 = yn .
6
(2.6)
無線部
Scala で実装する機械学習入門
開発班
活性化関数 f には式 (2.7) で定義される標準シグモイド関数を使用する。
fs (x) =
1
.
1 + e−x
(2.7)
シグモイド関数は階段関数を連続関数で近似したような 0 から 1 へのゆるやかな段差を描く (Fig. 2.4)。
Fig. 2.4: sigmoid function.
シグモイド関数を用いることで、パーセプトロンの学習に勾配法を適用する際、シグモイド関数の導関数が
シグモイド関数そのもので表現できるという強みが活かせる。
dfs
e−x
(x) =
= fs (x){1 − fs (x)}.
dx
(1 + e−x )2
(2.8)
今回は隠れ層が 1 層だけの 3 層パーセプトロンを実装しよう。手始めに Data クラスを定義する。
mlp.scala
case class Data(input: Seq[Double], answer: Seq[Double])
例えば、排他的論理和の訓練データは以下のように定義する。
mlp.scala
val samples = Seq(
Data(Seq(0.0, 0.0),
Data(Seq(0.0, 1.0),
Data(Seq(1.0, 0.0),
Data(Seq(1.0, 1.0),
Seq(0.0)),
Seq(1.0)),
Seq(1.0)),
Seq(0.0)))
この訓練データをコンストラクタの引数に取る MLP クラスを定義する。
mlp.scala
import scala.collection.mutable.Map
class MLP(samples: Seq[Data], I: Int, H: Int, O: Int) {
val w1 = Map[(Int, Int), Double]()
val w2 = Map[(Int, Int), Double]()
}
1
)で
I は入力層、H は隠れ層、O は出力層の階数である。入力層から隠れ層への伝搬は I 行 H 列の行列 (wih
2
重み付けされ、隠れ層から出力層への伝搬は H 行 O 列の行列 (who
) で重み付けされる。学習により重みが
変わることを考慮し、mutable.Map で重みを保持する。隠れ層と出力層の出力を式 (2.5) に基づき実装する。
7
無線部
Scala で実装する機械学習入門
開発班
mlp.scala
def hidden(args: Seq[Double]): Seq[Double] = {
val hid = new Array[Double](H)
for(i <- 0 until I; h <- 0 until H) hid(h) += w1((i, h)) * args(i)
for(h <- 0 until H) hid(h) = sigmoid(hid(h))
hid
}
def output(args: Seq[Double]): Seq[Double] = {
val out = new Array[Double](O)
for(h <- 0 until H; o <- 0 until O) out(o) += w2((h, o)) * args(h)
for(o <- 0 until O) out(o) = sigmoid(out(o))
out
}
hidden メソッドと output メソッドは、第 n 層の入力ベクトル xn を受け取って、ベクトル yn を出力する。
両関数とも sigmoid 関数を呼び出しているが、式 (2.4) に示した標準シグモイド関数の定義通りに実装する。
mlp.scala
def sigmoid(x: Double) = 1.0 / (1.0 + Math.exp(-x))
式 (2.5) に従って、入力ベクトル x1 を受け取り、出力ベクトル yN を計算する predict メソッドを定義する。
mlp.scala
def predict(input: Seq[Double]) = output(hidden(input))
n
以上で、3 層パーセプトロンは分類器として利用可能になった。問題は、どのような手順で重み付け (wij
)を
学習するかである。基本的な考え方は単層パーセプトロンと同じで、最終段の出力誤差が最小になる重みを
勾配法により探索する。ここに、出力誤差を徐々に最前段に伝搬させる誤差逆伝搬法を導入する。最終段の
出力 youtput と正解 yanswer との二乗誤差 E を式 (2.9) で定義する。
E = ||youtput − yanswer ||2 .
(2.9)
第 n 番目の層で、i 次元目から j 次元目への経路 (i → j) の重み wnij を最急降下法により最適化する。
wn′ij = wnij − η
∂E
∂wnij
.
(2.10)
合成関数の微分の公式を用いて式 (2.10) を展開すると、
∂E
∂wnij
=
∂E ∂nij
n
∂nij
n
∂wnij
=
∂E
∂nij
n
xin .
(2.11)
ij i
ただし、xin は第 n 層の入力 xn の i 次元目で、nij
n は経路 (i → j) の寄与 wn xn である。
∂E
∂nij
n
=
∂E ∂ynj
∂ynj ∂nij
n
=
∂E
∂
∂ynj ∂nij
n
fs (nij
n ).
(2.12)
シグモイド関数 fs の定義式 (2.7) を思い出しつつ第 n 層の出力 yn の j 次元目を ynj と定義すると、
∂
e−nn
ij
f (nij ) =
ij s n
∂nn
ij
(1 + e−nn )2
= ynj (1 − ynj ).
(2.13)
ここから、第 n + 1 層から第 n 層への誤差逆伝搬の漸化式を導出することができる。
∂E
∂ynj
=
K
jk
∑
∂E ∂nn+1
k=1
j
∂njk
n+1 ∂yn
=
K
∑
∂E k
jk
k
yn+1 (1 − yn+1
)wn+1
.
k
∂y
n+1
k=1
8
(2.14)
無線部
Scala で実装する機械学習入門
開発班
漸化式 (2.14) の境界条件として、最終段での二乗誤差の勾配を用いる。
∂E
j
∂youtput
j
j
= 2(youtput
− yanswer
).
(2.15)
最終段の誤差の勾配が出力側から入力側へと漸化式 (2.14) に従って伝搬していく様子が朧げにも想像できる
だろうか。数式を並べた解説は以上で終わりにして、さっそく実装してみよう。
mlp.scala
for(i <- 0 until I; h <- 0 until H) w1((i, h)) = 2 * Math.random - 1
for(h <- 0 until H; o <- 0 until O) w2((h, o)) = 2 * Math.random - 1
重みを乱数で初期化する作業を忘れずに入れておこう。全ての重みが 0 だと誤差が逆伝搬できない。
mlp.scala
val eta = 0.0
for(step <- 1 to 10000) samples.foreach{case Data(input , ans) => {
val hid = hidden(input)
val out = output(hid)
for(h <- 0 until H) {
val e2 = for(o <- 0 until O) yield ans(o) - out(o)
val g2 = for(o <- 0 until O) yield out(o) * (1 - out(o)) * e2(o)
for(o <- 0 until O) w2((h, o)) += eta * g2(o) * hid(h)
val e1 = for(o <- 0 until O) yield g2(o) * w2((h, o))
val g1 = for(i <- 0 until I) yield hid(h) * (1 - hid(h)) * e1.sum
for(i <- 0 until I) w1((i, h)) += eta * g1(i) * input(i)
}
}}
以上を先掲の MLP クラスのコンストラクタ内に追記する。さっそく排他的論理和を学習させてみよう。
2.0
1.5
1.5
0.8
0
2.0
0.8
y
0
0.5
0.4
0.20
1.0
0.60
0.5
y
0
0.40
0.6
1.0
0.4
0
0.8
60
0.
0.6
0
0.0
0
0
0.0
−0.5
−0.5
0.0
0.5
x
80
0.
−1.0
−1.0
−0.5
1.0
1.5
2.0
(a) 10,000 steps
−1.0
−1.0
0.
20
0.60
0.
20
0.4
−0.5
0
0.0
0.4
0.5
x
0
1.0
0.2
0
1.5
2.0
(b) 1,000,000 steps
Fig. 2.5: exclusive OR emulated by a 3-layer perceptron.
Fig. 2.5 は、排他的論理和を値域が [0, 1] の実関数として学習した際の出力値の分布である。活性化関数 f に
シグモイド関数を利用しているので、値域は [0, 1] の範囲に収まる。真偽値が欲しい場合は 0.5 を閾値として
分離すれば良い。今回は誤差逆伝搬を 10,000 回繰り返す実装としたが、繰り返し数を増やすと徐々に分布が
変化する。Fig. 2.5(b) は 1,000,000 回繰り返した結果だが、左下に島が出現している。最終的にこの状態で
収束するため、式 (2.9) の誤差関数 E が十分に小さくなった時点で学習を打ち切ると良い。
9
第 3 章 サポートベクターマシン
3.1
ハードマージンの場合
サポートベクターマシンは、ふたつの集団を分離する識別直線を学習する機械である。Fig. 3.1 を見よ。
Fig. 3.1: decision surface with hard margin.
サポートベクターマシンが学習するのは単なる境界線ではない。正集団と負集団の双方からの最短距離 d が
最大になる直線を学習するのである。仮に、既知の訓練データよりも識別直線に近い場所に未知のデータが
出現した場合でも、反対側の集団に誤分類される危険性は低くなる、という期待ができる。既にお気付きの
読者もいるかもしれないが、これは制約付き最大値問題である。点 (xi , yi ) の制約条件は式 (3.1) で与える。
|axi + byi + c| ≥ 1.
(3.1)
ただし、a, b, c は識別直線の係数で、任意の実数値である。集団と識別直線との最短距離は式 (3.2) で求まる。
min d(xi , yi ) = min
|axi + byi + c|
1
√
=√
.
a2 + b2
a2 + b2
(3.2)
あとはラグランジュの未定乗数法で a2 + b2 を最小化するだけである。正解を ti = ±1 として式 (3.3) を得る。
L(a, b, c, λ) =
N
∑
1 2
(a + b2 ) −
λi {ti (axi + byi + c) − 1}.
2
i=1
(3.3)
式 (3.3) を a, b, c で偏微分することにより、L が最小になり最短距離 d が最大になる条件が求まる。
a=
N
∑
λi ti xi , b =
i=1
N
∑
i=1
λi ti yi , 0 =
N
∑
λi ti .
N
1∑
L̃(λ) = min L(a, b, c, λ) =
λi 1 −
λj ti tj (xi xj + yi yj ) .
a,b,c
2 j=1
i=1
N
∑
10
(3.4)
i=1
(3.5)
無線部
Scala で実装する機械学習入門
開発班
あとは双対問題 L̃(λ) を最大化するだけである。最急上昇法による学習が数多く紹介されているが、今回は
逐次最小問題最適化法を採用する。これまで未定乗数法を利用して議論してきたが、厳密には、式 (3.1) の
制約式は不等式であるため、KKT(Karush–Kuhn–Tacker) 条件を正しく適用する必要がある。KKT 条件を
適用することにより、式 (3.5) の最適化問題の解が存在する必要条件は、式 (3.6) で与えられる。
λi {ti (axi + byi + c) − 1} = 0, (λi ≥ 0).
(3.6)
仮に λi > 0 である場合、axi + byi + c = ti が成立する必要がある。このことは即ち、(xi , yi ) が Fig. 3.1 の
点線上に存在することを意味する。逆に、式 (3.4) から明らかなように、λi = 0 の点 (xi , yi ) は、係数 a, b の
決定に何ら寄与しない。サポートベクターマシンの学習では、Fig. 3.1 の識別直線の最近傍のデータのみが
考慮される。このようなデータをサポートベクトルと呼び、サポートベクターマシンの卓越した汎化能力の
所以である。さて、逐次最小問題最適化法では、式 (3.6) の条件を破る点 (xi1 , yi1 ) が存在する限り、任意の
点 (xi2 , yi2 ) を選び、式 (3.4) を根拠として、式 (3.7) を満たす局所的な最適化を行う。
new
old
old
λnew
i1 ti1 + λi2 ti2 = λi1 ti1 + λi2 ti2 .
(3.7)
ここで、1 回の最適化による λi1 , λi2 の増分 δi1 , δi2 を式 (3.8)(3.9) で定義する。
δi1 =λnew
− λold
i1
i1 ,
δi2 =λnew
i2
−
λold
i2 .
(3.8)
(3.9)
式 (3.5) の双対問題 L̃(λi1 , λi2 ) のうち、δi1 , δi2 が関わる部分を抜き出すと、式 (3.10) を得る。
dL̃(δi1 , δi2 ) = + δi1 + δi2 − δi1 δi2 ti1 ti2 (xi1 xi2 + yi1 yi2 )
− δi1 ti1
− δi2 ti2
N
∑
1 2 2 2
2
ti1 (xi1 + yi1
)
λn tn (xi1 xn + yi1 yn ) − δi1
2
n=1
(3.10)
N
∑
1 2 2 2
2
λn tn (xi2 xn + yi2 yn ) − δi2
ti2 (xi2 + yi2
).
2
n=1
式 (3.10) について、特に ti1 = ti2 の場合、式 (3.7) の条件より δi1 = −δi2 が導かれるので、式 (3.11) を得る。
2
(xi1 xi2 + yi1 yi2 )
dL̃(δi1 ) = + δi1
− δi1 ti1
+ δi1 ti1
N
∑
1 2 2
2
λn tn (xi1 xn + yi1 yn ) − δi1
(xi1 + yi1
)
2
n=1
(3.11)
N
∑
1 2 2
2
λn tn (xi2 xn + yi2 yn ) − δi1
(xi2 + yi2
).
2
n=1
L̃(λ) の最大化は、式 (3.10) の極大点を探すことと同義である。式 (3.10) を δi1 で微分し、式 (3.12) を得る。
[N
]
∑
ti1
δi1 = −
λn tn {(xi1 − xi2 )xn + (yi1 − yi2 )yn } − ti1 + ti2 .
(3.12)
(xi1 − xi2 )2 + (yi1 − yi2 )2 n=1
式 (3.12) が成立するとき式 (3.10) が極大値を取る。ti1 = −ti2 の場合も式 (3.12) を得る。ti1 − ti2 なる項を
敢えて追加することで、ti1 = ti2 の場合でも、ti1 ̸= ti2 の場合でも、共通して式 (3.12) を利用できるように
設計されている。訓練データ (xi , yi ) が式 (3.6) の KKT 条件を満たすかどうかを調べるには、係数 a, b, c が
必要である。係数 a, b は式 (3.4) により求まるが、バイアス項 c だけ求められない。ここで、ti (axi + byi ) が
最小になる (xi , yi ) は、Fig. 3.1 の識別直線の最近傍であることに着目しよう。バイアス項 c の値を近似的に
求めるには、そのようなサポートベクトルの候補を探し出して、式 (3.13) を計算すれば良い。
{
}
1
c=−
min (axi + byi ) + max (axj + byj ) .
2 i|ti =+1
j|tj =−1
11
(3.13)
無線部
Scala で実装する機械学習入門
開発班
以上を踏まえ、式 (3.12) を計算すれば λi1 と λi2 の値が更新できるのであるが、ここで油断してはならない。
更新後のラグランジュ乗数もまた、式 (3.6) の条件を満たす必要がある。具体的には、非負性が求められる。
従って、ti1 = ti2 の場合は、式 (3.14)(3.15) を満たす必要がある。
λnew
=λold
i1
i1 + δi1 ≥ 0,
(3.14)
λnew
=λold
i2
i2 − δi1 ≥ 0.
(3.15)
同様に、ti1 ̸= ti2 の場合は、式 (3.16)(3.17) を満たす必要がある。
λnew
=λold
i1
i1 + δi1 ≥ 0,
(3.16)
λnew
=λold
i2
i2 + δi1 ≥ 0.
(3.17)
条件を満たさない場合は δi1 を切り抜く。具体的には、λi1 が負になる場合は δi1 の値を λi1 = 0 となる点で
クリップすれば良い。λi2 が負になる場合も同様である。以上で、サポートベクターマシンの学習に必要な
数式は全て出揃った。さっそく実装に入ろう。手始めに、1 件の訓練データを表す Data クラスを定義する。
svm.scala
case class Data(x: Double , y: Double , t: Int) {
var l = 0.0
}
t は正解ラベル ti で、l はラグランジュ乗数 λi である。続いて SVM クラスを定義する。
svm.scala
class SVM(samples: Seq[Data]) {
var (a, b, c) = (0.0, 0.0, 0.0)
def predict(x: Double , y: Double) = a * x + b * y + c > 0
}
式 (3.13) に従って係数 c の値を求める newc メソッドを SVM クラスに定義する。
svm.scala
def newc = {
val pos = samples.filter(_.t == +1).map(d => a * d.x + b * d.y).min
val neg = samples.filter(_.t == -1).map(d => a * d.x + b * d.y).max
(pos + neg) / 2
}
式 (3.4) に従って係数 a と b を求める newa メソッドと newb メソッドも定義する。
svm.scala
def newa = (for (d <- samples) yield d.l * d.t * d.x).sum
def newb = (for (d <- samples) yield d.l * d.t * d.y).sum
訓練データが式 (3.1) と式 (3.6) を満たすことを確認する kkt メソッドを定義する。
svm.scala
def kkt(d: Data) = d.l match {
case 0 => d.t * (a * d.x + b * d.y + c) >= 1
case _ => d.t * (a * d.x + b * d.y + c) == 1
}
12
無線部
Scala で実装する機械学習入門
開発班
式 (3.12) に基づいてラグランジュ乗数 λi1 と λi2 を更新する update メソッドを定義する。
svm.scala
def update(d1: Data, d2: Data)
val (dx, dy) = (d1.x - d2.x,
val den = dx * dx + dy * dy
val num = (for(d <- samples)
val di1 = clip(d1, d2, -d1.t
d1.l += di1
d2.l -= di1 * d1.t * d2.t
}
{
d1.y - d2.y)
yield d.l * d.t * (dx * d.x + dy * d.y)).sum
* (num - d1.t + d2.t) / den)
式 (3.14)(3.15)(3.16)(3.17) に基づき δi1 , δi2 をクリップするための clip メソッドも定義しておく。
svm.scala
def clip(d1: Data, d2: Data, di1: Double): Double = {
if (d1.t == d2.t) {
if (di1 < -d1.l) return -d1.l
if (di1 > +d2.l) return +d2.l
return di1
} else return Seq(-d1.l, -d2.l, di1).max
}
ここまで来れば、式 (3.6) を満たさない訓練データを見つけて update メソッドに渡すだけである。
svm.scala
while (samples.exists(!kkt(_))) {
val i1 = samples.indexWhere(!kkt(_))
var i2 = 0
do i2 = (math.random * samples.size).toInt while(i1 == i2)
update(samples(i1), samples(i2))
a = newa
b = newb
c = newc
}
ついにサポートベクターマシンの完成である。動作検証を兼ねて、Fig. 3.2 に示すデータを訓練させてみた。
2.0
1.5
1.0
y
0.5
0.0
0.0
−0.5
0
−1.0
−1.5
−2.0
−2.0 −1.5 −1.0 −0.5
0.0
x
0.5
1.0
1.5
2.0
Fig. 3.2: decision surface learned by a hard margin SVM.
正と負の両集団から等距離で、円の中心を結んだ線分と垂直に交差する識別直線が描けていれば成功である。
13
無線部
3.2
Scala で実装する機械学習入門
開発班
ソフトマージンの場合
第 3.1 節に述べたサポートベクターマシンは、正負の集団を完璧に分離する識別直線を習得するまで学習を
継続する。識別直線は全くの誤分類なく、正しく正負を分離する。逆に言えば、サポートベクターマシンは
分離困難な訓練データの出現を適切に想定していない。現実のデータは大域的には単純な分布を示すように
見えても、局所的には入り組んだ分布を示すことがある。Fig. 3.3 に示すように、赤い点で示した正集団に
1 点だけ負の青いデータ点が混ざっている状況は、その例である。
Fig. 3.3: decision surface with soft margin.
Fig. 3.3 の識別直線は、青い点を正しく分離できていない。この点は負の訓練データであるにもかかわらず
正の訓練データとして誤分類されてしまう。サポートベクターマシンは、この誤分類を強引に解消しようと
試みる。しかし、識別直線を跨いで正負が入り組んだ分布になっている場合、必ずしも分離に成功するとは
言えない。縦しんば完璧な分離に成功したとしても、過学習に陥る危険がある。過学習とは、訓練データの
分類精度を限界まで高めようと努力する余り、汎化性能を損ねる状況である。言い換えれば、木を見て森を
見ずの状況であり、未知のデータを適切に分類する能力は犠牲になる。過学習はサポートベクターマシンに
限らず機械学習全般に共通する問題であるが、対策として目先の誤分類に囚われず寛容の精神で許すことが
求められる。サポートベクターマシンの場合はソフトマージンの導入が有効である。これは多少の誤分類を
許容する識別直線である。ただし誤分類された訓練データ xi に対しては個別にペナルティξi を課すものと
する。この設定のもとで、式 (3.1) の制約条件は、式 (3.18) に書き換えられる。
ti (axi + byi + c) ≥ 1 − ξi
(3.18)
ペナルティξi は各訓練データに割り当てられるスラック変数で、式 (3.19) で定義する。
0
if ti (axi + byi + c) > 1,
ξi =
|t − (ax + by + c)| ≥ 1 if t (ax + by + c) ≤ 1.
i
i
i
i
i
(3.19)
i
最小化すべき目的関数は、定数 C > 0 を導入することにより、式 (3.20) で与えられる。
C
N
∑
1
ξi + (a2 + b2 ).
2
i=1
(3.20)
∑
任意の誤分類点 (xi , yi ) について ξi ≥ 1 が成立することから、 i ξi は誤分類された訓練データの個数よりも
∑
大きな値を取る。逆に言えば、 i ξi は誤分類された個数の上限を与える。定数 C は訓練データの誤分類を
許容する程度を間接的に決定する定数として作用することが期待できる。少なくとも、C → ∞ の極限では
14
無線部
Scala で実装する機械学習入門
開発班
誤分類は全く許容されない。従って、第 3.1 節で扱ったハードマージンは、ソフトマージンの特殊な場合と
見なすことができる。式 (3.3) のラグランジュ関数は、式 (3.21) に書き換える。
L(a, b, c, ξ, λ, µ) =
N
N
N
∑
∑
∑
1 2
(a + b2 ) + C
ξi −
λi {ti (axi + byi + c) − 1} −
µ i ξi .
2
i=1
i=1
i=1
(3.21)
式 (3.3) の未定乗数は λi だけであったが、今回は λi と µi が未定乗数である。今回も前回と同様に a, b, c で
偏微分することにより、L が最小になりマージン d が最大化する条件が求まる。
a=
N
∑
λi ti xi , b =
i=1
N
∑
λi ti yi , 0 =
i=1
N
∑
λi ti , λi = C − µi ,
(3.22)
i=1
N
1∑
L̃(λ) = min L(a, b, c, λ) =
λi 1 −
λj ti tj (xi xj + yi yj ) .
a,b,c
2 j=1
i=1
N
∑
(3.23)
式 (3.23) は式 (3.5) と同等であるが、KKT 条件から導かれる必要条件は式 (3.6) と異なる点に注意を要する。
λi {ti (axi + byi + c) + ξi − 1} = 0, (λi ≥ 0, µi ≥ 0, µi ξi = 0).
(3.24)
また、式 (3.22) より次式の制約条件が導出される。
0 ≤ λi ≤ C.
(3.25)
以上の議論に基づき、逐次最小問題最適化法の適用を試みる。まず λi = 0 となる点 (xi , yi ) は式 (3.22) より
µi = C を満たす。C > 0 なので、式 (3.24) を満たすには ξi = 0 が必要である。これは式 (3.19) の定義から
|axi + byi + c| > 1 と同等である。次に、λi = C となる点 (xi , yi ) は式 (3.22) より µi = 0 を満たす。ここで
ξi > 0 となる可能性があり、式 (3.24) により ti (axi + byi + c) ≤ 1 が成立する。最後に、0 < λi < C となる
点 (xi , yi ) は µi > 0 ゆえ ξi = 0 を満たす。この点 (xi , yi ) は正しく分類され、しかも ti (axi + byi + c) = 1 の
境界線上に存在する。以上の考察より、第 3.1 節で実装した kkt メソッドを以下のように書き換える。
svm.scala
def kkt(d: Data) = d.l
case 0 => d.t * (a *
case C => d.t * (a *
case _ => d.t * (a *
}
match
d.x +
d.x +
d.x +
{
b * d.y + c) >= 1
b * d.y + c) <= 1
b * d.y + c) == 1
式 (3.23) は制約条件を除き式 (3.5) と同等なので、式 (3.12) に示した最適化の議論をそのまま流用してよい。
ただし、式 (3.25) の制約を守るために切り抜きが必要になる。ti1 = ti2 の場合、式 (3.26)(3.27) を遵守する。
0 ≤ λnew
=λold
i1
i1 + δi1 ≤ C
(3.26)
0 ≤ λnew
=λold
i2
i2 − δi1 ≤ C
(3.27)
同様に ti1 ̸= ti2 の場合は、式 (3.28)(3.29) を遵守する。
0 ≤ λnew
=λold
i1
i1 + δi1 ≤ C
(3.28)
0≤
(3.29)
λnew
i2
=λold
i2
+ δi1 ≤ C
さっそく、第 3.1 節で定義した clip メソッドを書き換えよう。まず、SVM クラスの冒頭に定数 C を宣言する。
svm.scala
class SVM(samples: Seq[Data]) {
val C = 0.01
}
15
無線部
Scala で実装する機械学習入門
開発班
式 (3.26)(3.27)(3.28)(3.29) に基づき、clip メソッドを以下のように定義する。
svm.scala
def clip(d1: Data, d2: Data, di1: Double): Double = {
if (d1.t == d2.t) {
val L = Math.max(-d1.l, d2.l - C)
val H = Math.min(+d2.l, C - d1.l)
if (di1 < L) return L
if (di1 > H) return H
return di1
} else {
val L = Math.max(-d1.l, -d2.l)
val H = Math.min(-d1.l, -d2.l) + C
if (di1 < L) return L
if (di1 > H) return H
return di1
}
}
以上で、ソフトマージンの実装は完了である。動作検証を兼ねて、Fig. 3.4 に示すデータを訓練させてみた。
2.0
1.5
1.0
y
0.5
0
0.0
.0
0
−0.5
−1.0
−1.5
−2.0
−2.0
−1.5
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
x
Fig. 3.4: decision surface learned by a soft margin SVM.
定数 C の値によって、学習の様子がどのように変化するかを観察してみると面白い。大きな値に設定すると
学習が終わらなくなり、無限ループに陥る。小さな値に設定すれば誤分類データは遺棄され、学習は早期に
終結する。Fig. 3.4 の訓練データは単純なので、境界部分を跨いで複雑に両集団が入り組んだ訓練データを
用意すると良い。興味深いことに、Fig. 3.4 は Fig. 3.2 に比べ、識別直線が負の集団に接近している。実は
サポートベクターマシンには学習の結果が外れ値に左右される難点がある。Fig. 3.4 の例ではふたつの円の
中心にある誤分類点が外れ値であるが、特に識別直線から遠い右上の後分類点に識別直線が引き寄せられて
いる。結果として、訓練データに対する分類能力はこの正負の外れ値を除き正確に保たれるが、汎化能力は
低下する恐れがある。
16
第 4 章 カーネルトリック
第 3 章で解説したサポートベクターマシンは、多層パーセプトロンとは異なり、線型分離不可能な問題には
適用できない。しかし、現にサポートベクターマシンは線型分離不可能な問題にも適用され、最も普及した
機械学習モデルになっている。実は、線型分離不可能な問題を線型分離可能な問題へと変えてしまう便利な
方法があるのだ。本章ではその理論的背景とコーディング例を紹介する。
4.1
高次元への写像
線型分離不可能な問題は、特徴空間を高次元空間に写像すれば、線型分離可能になる可能性がある。次元が
増大するに伴いデータの分布は粗になり、ついには線型分離可能になる。Fig. 4.1 は、2 次元平面上の分布を
例に、この原理を視覚化したものである。同心円に並んだ集団は、直線では分離できないが、座標 (x, y) を
ガウス関数 fg により座標 (x, y, fg (x, y)) に投影すれば、z 軸に垂直な平面で分離できる。
2.0
1.5
1.0
0.5
0.0
z
−0.5
−1.0
−1.5
−2.0
−4
−3
−2
−1
x
0
1
2
3
4
−4
−3
−2
−1
0
1
2
3
4
y
Fig. 4.1: conversion of linearly inseparable problem to separable.
この原理に基づき、第 3.2 節のサポートベクターマシンを改造する。手始めに、高次元空間 Ω への写像 Φ を
定義する。Fig. 4.1 の例では、Φ : (x, y) 7−→ (x, y, fg (x, y)) であるが、Φ の定義については後述する。
Φ : x 7→ Φx (Φx ∈ Ω).
17
(4.1)
無線部
Scala で実装する機械学習入門
開発班
第 3 章で導出した式は、訓練データ (xi , yi ) を Φx に置き換えれば、本章の議論に流用できる。抽象化のため
今後は (xi , yi ) をベクトル x で、係数 a, b をベクトル w で表記する。第 3 章の式のうち、学習時に計算する
必要がある式を列挙すると、式 (3.24)(3.12)(3.13) は式 (4.2)(4.3)(4.4) で再定義できる。⟨·, ·⟩ は内積である。
λi {ti (w · Φx + c) + ξi − 1} = 0, (λi ≥ 0, µi ≥ 0, µi ξi = 0),
∑
n λn tn ⟨Φxi1 − Φxi2 , Φxn ⟩ − ti1 + ti2
δi1 = −ti1
,
⟨Φxi1 , Φxi1 ⟩ − 2⟨Φxi1 , Φxi2 ⟩ + ⟨Φxi2 , Φxi2 ⟩
{
}
1
c=−
min w · Φxi + max w · Φxj .
2 i|ti =+1
j|tj =−1
(4.2)
(4.3)
(4.4)
式 (3.22) により、係数 w が全訓練データに対する総和で計算できることを思い出そう。
w=
N
∑
λi ti Φxi .
(4.5)
i=1
従って、式 (4.2) と式 (4.4) は、内積 ⟨Φxi , Φxj ⟩ を用いて、式 (4.6)(4.7) のように再定義できる。
N
∑
λj tj ⟨Φxi , Φxj ⟩ + c + ξi − 1 = 0,
λi ti
(4.6)
j=1
N
N
∑
∑
1
c=−
. min
λj tj ⟨Φxi , Φxj ⟩ + max
λi ti ⟨Φxi , Φxj ⟩ .
2 i|ti =+1 j=1
j|tj =−1
i=1
(4.7)
第 4.3 節に述べるが、内積を用いることでカーネルトリックと呼ばれる技法の恩恵を享受できる。そもそも
高次元空間への写像 Φ はサポートベクターマシンの学習に大きな計算量を要求する。Φx の次元に比例する
空間計算量を要求する上、時間計算量もそれに比例するためである。カーネルトリックは、高次元空間 Ω で
行うべき計算を、写像前の空間での計算で代用する。ただし、代用できるのは内積 ⟨Φxi , Φxj ⟩ に限られる。
4.2
ヒルベルト空間
カーネル法の議論をする際には、数学的な議論が必要になり、馴染みのない大量の数学用語に困惑するかも
しれない。学生は、どんな目的意識でその理論を持ち出すのか理解できないので、講義を真面に聞くことを
諦めるだろう。本節では、厳密ではないながらも議論の流れを追える程度に整理しつつカーネル法の理解を
目指す。まず、距離が定義された空間である距離空間を導入する。ユークリッド空間から距離の概念のみを
抽出した概念と言える。空間内の任意の 2 点 x, y の距離を表す距離関数 d は式 (4.8) の性質を満たす。
d(x, y) ≥ 0,
x = y ⇔ d(x, y) = 0,
(4.8)
d(x, y) = d(y, x),
d(x, y) ≤ d(x, z) + d(z, y).
式 (4.8) の性質は、それぞれ正定値性、非退化性、対称性、劣加法性と呼ばれる。さて、距離空間ではまだ
抽象的すぎるので、制約を加えて完備距離空間を定義する。距離空間 M が完備であるとは、M 内の任意の
コーシー列が M に属する点に収斂する様子である。コーシー列は式 (4.9) の制約を満たす列である。
lim ||xn − xm || = 0.
n,m→∞
18
(4.9)
無線部
Scala で実装する機械学習入門
開発班
完備距離空間とは別に内積空間も定義する。内積空間ないし計量ベクトル空間とは、文字通り内積を備えた
空間である。ただし、高校数学で習う標準内積よりも抽象化された内積である。まず、対称半双線型形式を
定義する。対称半双線型形式ないしエルミート形式とは、式 (4.10) の 2 性質を満たす写像 ⟨·, ·⟩ である。
⟨x, ay + z⟩ = a⟨x, y⟩ + ⟨x, z⟩,
(4.10)
⟨x, y⟩
= ⟨y, x⟩.
前者は線型性で、後者はエルミート対称性である。式 (4.10) に式 (4.11) の 2 性質を加えたものを内積と呼ぶ。
⟨x, x⟩
≥ 0,
(4.11)
⟨x, x⟩ = 0 ⇒ x = 0.
前者は正定値性であり、後者は非退化性である。高校数学で習う標準内積は、式 (4.10)(4.11) の性質を満たす。
x1
y1
∑
n
x2 y2
x·y = . · . =
xi yi .
.. .. i=1
xn
(4.12)
yn
また、関数 f と g を元とする空間 Ω における内積 ⟨f, g⟩ は測度を µ で表記すると、式 (4.13) で定義される。
∫
⟨f, g⟩ =
f (x)g(x)dµ(x).
(4.13)
Ω
さて、完備距離空間と内積空間を定義したことで、この双方を満たすクラスであるヒルベルト空間の定義が
可能になる。ヒルベルト空間 H は、内積空間であり、かつ完備な距離空間である。実は、ヒルベルト空間は
何ら特殊な空間ではなく、我々が生活する空間を抽象化した概念であるに過ぎない。第 4 章の議論において
重要なのは、内積と距離が定義された空間であるということである。
4.3
再生核カーネル
式 4.14) のように、関数 f を入力に取り関数 Tf を出力するような変換 T を積分変換と呼ぶ。
∫
Tf (s) =
k(s, t)f (t)dt, (s ∈ Ωs ).
(4.14)
Ωt
変換前の関数 f は t を引数とする関数であり、変換後の関数 Tf は s を引数とする関数である。積分変換は
領域 Ωt での求解が難しい問題を別の領域 Ωs に投影して容易に求解し、Ωs から Ωt に逆変換することで解を
得る技法で多用される。フーリエ変換 F やラプラス変換 L がその代表例である。領域 Ωt から領域 Ωs への
仲介を行う関数 k(s, t) をカーネルと呼ぶ。フーリエ変換の場合は e−jωt が、ラプラス変換の場合は e−st が
該当する。中でも再生性を有するカーネルを特に再生核と呼ぶ。再生性とは、式 (4.15) のように積分変換の
前後で同じ関数 f が出現するような性質である。
∫
f (s) =
k(s, t)f (t)dt.
(4.15)
Ωt
再生核を備えたヒルベルト空間を特に再生核ヒルベルト空間と呼ぶ。再生核を使った写像 Φ を考える。
Φ : x 7−→ k(x, xi ).
19
(4.16)
無線部
Scala で実装する機械学習入門
開発班
ただし xi は訓練データの 1 件とする。写像 Φ を用いて、訓練データ xj を高次元ベクトル Φxj に写像する。
ここで、ベクトル Φxj の線型結合 f (xi ) により構成するベクトル空間は、再生核ヒルベルト空間 Hk である。
f (xi ) =
N
∑
aj Φxj =
j=1
N
∑
aj k(xj , xi ).
(4.17)
j=1
再生核ヒルベルト空間 Hk 内でベクトル Φxi とベクトル Φxj との内積 ⟨Φxi , Φxj ⟩ を考えてみよう。
⟨Φxi , Φxj ⟩ = ⟨k(xi , xl ), k(xl , xj )⟩.
ここで関数を元に持つ空間での内積の定義を思い出そう。式 (4.13) より、式 (4.19) が導かれる。
∫
⟨k(xi , xl ), k(xl , xj )⟩ =
k(xi , xl )k(xl , xj )dxl .
(4.18)
(4.19)
xl
カーネル k は再生核であるため、式 (4.15) の再生性に従い、式 (4.20) が導かれる。
∫
k(xi , xl )k(xl , xj )dxl = k(xi , xj ).
(4.20)
xl
ここから、写像 Φ による変換後のベクトル間の内積 ⟨Φxi , Φxj ⟩ は、Φ を経由せずに求まることが判明する。
⟨Φxi , Φxj ⟩ = k(xi , xj ).
(4.21)
この結論は、高次元空間での計算が不要になるという以上の恩恵を我々にもたらす。写像 Φ の定義が困難な
問題でも、カーネルさえ定義できれば線型分離できるのだ。さて、以上でカーネルトリックの理論的背景は
網羅した。この辺りで具体的に既知のカーネルを見ておこう。まず、線型カーネルを紹介する。
kl (xi , xj ) = xi · xj =
D
∑
xid xjd .
(4.22)
d=1
大仰な名前が付いているが、要するに標準内積である。既にお気付きだろうが、線型カーネルとは第 3 章で
解説したサポートベクターマシンそのものである。続いて、多項式カーネルを紹介する。
( D
)p
∑
xid xjd + b .
kp (xi , xj ) = a
(4.23)
d=1
多項式カーネルの利点は、写像後の高次元空間に特徴量の積が出現することである。これにより、特徴量の
積を基底に含めることができる。例えば、2 次元ベクトル x, x′ に対し (x · x′ )2 の基底は、式 (4.24) になる。
2
x1
0
0
(4.24)
0 , x1 x2 , 0 .
0
0
x22
続いて、線型分類器では最も汎用的に用いられるラジアル基底関数もしくはガウシアンカーネルを紹介する。
(
)
|xi − xj |2
krbf (xi , xj ) = exp −
.
(4.25)
2σ 2
ガウシアンカーネルの特徴は、無限次元の特徴空間への写像を表すことにある。試しに級数展開してみよう。
xi · xi
xi · xj
xj · xj
|xi − xj |2
=+
−
+
,
2σ 2
2σ 2
σ2
2σ 2
D
D
∑
∏
x·y
xd yd
xd yd
exp 2 = exp
=
exp 2 ,
2
σ
σ
σ
d=1
exp
(4.27)
d=1
∞
∑
1 ( xd )n 1 ( yd )n
xd yd
√
√
=
.
2
σ
n! σ
n! σ
n=0
20
(4.26)
(4.28)
無線部
Scala で実装する機械学習入門
式 (4.28) は式 (4.29) の要素を持つ無限次元のベクトルの内積を表現することになる。
1 ( x d )n
1 ( yd )n
√
, √
, (n ≥ 0).
n! σ
n! σ
開発班
(4.29)
シグモイドカーネルは、厳密にはヒルベルト空間を構成しないものの、広く使われるカーネルである。
( D
)
∑
ks (xi , xj ) = tanh a
xid xjd + b .
(4.30)
d=1
ここで、サポートベクターマシンの動作を思い出そう。第 3 章にも述べた通り、式 (4.31) がその動作である。
sign(w · Φx + c).
(4.31)
式 (4.30) のシグモイドカーネルを式 (4.31) に適用することにより、第 2.2 節に述べた 3 層パーセプトロンを
模倣できる。入力層では未分類データ x を受け取り、隠れ層では内積 x · xi に式 (2.7) のシグモイド関数を
適用し、出力層では重み λi tj の線型和に定数項 c を足して符号を確認する。ぜひ確かめてみて欲しい。
4.4
分類器への適用
第 4.3 節までの議論を踏まえ、サポートベクターマシンにカーネル法を適用し、線型分離不可能な問題にも
利用可能なように改良する。副次的に、第 3 章の実装が 2 次元にしか対応していなかったのを、任意次元に
対応するようにも改良する。従って、Data クラスから作り直す必要がある。
svm.scala
case class Data(p: Seq[Double], t: Int) {
var l = 0.0
}
SVM クラスを任意のカーネルに対応させるべく、コンストラクタの引数にカーネル k を追加する。第 3 章の
実装では係数 a, b を変数として宣言していたが、今回の実装では使用しない。SVM クラスのフィールドには
未定乗数 L とソフトマージン係数 C、定数項 c のみを宣言する。
svm.scala
class
val
val
var
}
SVM(samples: Seq[Data], k: (Data, Data) => Double) {
L = new Array[Double](samples.size)
C = 1e100
c = 0.0
この SVM クラスの内部に以降のメソッドを追記していく。まず、w · Φx を計算する wx メソッドを定義する。
svm.scala
def wx(x: Data) = (for (s <- samples) yield s.l * s.t * k(x, s)).sum
wx メソッドを利用するように kkt メソッドを定義する。
svm.scala
def kkt(d: Data) = d.l match {
case 0 => d.t * (wx(d) + c) >= 1
case C => d.t * (wx(d) + c) <= 1
case _ => d.t * (wx(d) + c) == 1
}
21
無線部
Scala で実装する機械学習入門
開発班
逐次最小問題最適化法を実行する update メソッドも、第 3.2 節の clip メソッドを流用しつつ、再定義する。
svm.scala
def update(d1: Data, d2: Data) {
val den = k(d1, d1) - 2 * k(d1, d2) + k(d2, d2)
val num = wx(d1 - d2) - d1.t + d2.t
val di1 = clip(d1, d2, -d1.t * num / den)
d1.l += di1
d2.l -= di1 * d1.t * d2.t
}
式 (4.4) を用いて定数項 c を計算する newc メソッドを定義する。
svm.scala
def newc = {
val pos = samples.filter(_.t == +1).map(wx(_)).min
val neg = samples.filter(_.t == -1).map(wx(_)).max
(pos + neg) / 2
}
最後に、式 (3) に基づき、クラスを判定する predict メソッドを定義する。
svm.scala
def predict(p: Data) = wx(p) + c > 0
ここまで完成すれば、後は SVM クラスのコンストラクタで学習の処理を記述するだけである。
svm.scala
while (samples.exists(!kkt(_))) {
val i1 = samples.indexWhere(!kkt(_))
var i2 = 0
do i2 = (math.random * samples.size).toInt while(i1 == i2)
update(samples(i1), samples(i2))
c = newc
}
2.0
1.5
1.0
y
0.5
0.0
−0.5
−1.0
−1.5
−2.0
−2.0 −1.5 −1.0 −0.5
0.0
x
0.5
1.0
1.5
2.0
Fig. 4.2: decision surface learned by a SVM with a Gaussian kernel.
試しに、式 (4.25) のガウシアンカーネルを用いて、線型分離不可能な集団を分離した様子を Fig. 4.2 に示す。
22
第 5 章 単純ベイズ分類器
これまで本稿では、特徴空間に分布するデータを何らかの関数により分離する方式のアルゴリズムについて
議論してきた。第 5 章では視点を変え、条件付き確率を用いて結果から原因を推論する単純ベイズ分類器を
実装する。単純ベイズ分類器は、初歩的な話題判定アルゴリズムであり、しばしばスパムメールフィルタに
応用される。まず、文字列 D を単語 w の集合とし、文字列 D がクラス C に分類される確率を求める。
P (D|C) = P (w1 , . . . , wn |C) =
N
∏
P (wn |C).
(5.1)
n=1
本来ならば、単語間の共起関係も考慮し、式 (5.2) のような複雑かつ膨大で際限のない確率計算を必要とする。
P (w1 , . . . , wn |C) = P (w1 |C)P (w2 , . . . , wn |C, w1 ),
P (w2 , . . . , wn |C, w1 ) = P (w2 |C, w1 )P (w3 , . . . , wn |C, w1 , w2 ),
(5.2)
P (w3 , . . . , wn |C, w1 , w2 ) = P (w3 |C, w1 , w2 )P (w4 , . . . , wn |C, w1 , w2 , w3 ), . . .
単純ベイズ分類器では、単語間の独立性を仮定して計算を大胆に単純化している。文字列 D をクラス C が
生成する確率 P (D|C) を定義したので、次に、文字列 D の原因として最有力のクラス Cj を推測する方法を
考える。直観的には、そのようなクラスは P (C|D) を最大にする筈なので、式 (5.3) で与えられる。
j = arg max P (Ck |D) = arg max
k
k
P (D|Ck )P (Ck )
.
P (D)
(5.3)
P (D) が Ck に依存しないことに着目すると、式 (5.3) の代わりに式 (5.4) で代用することもできる。
j = arg max P (D|Ck )P (Ck ).
(5.4)
k
式 (5.1) の仮定より、単純ベイズ分類器の分類規則として式 (5.5) が導かれる。
j = arg max P (Ck )
k
N
∏
P (wn |Ck ).
(5.5)
n=1
観測不可能な事象 A が原因で観測データ B が得られたとき、条件付き確率 P (B|A) は B の原因が事象 A で
あるという仮説の尤もらしさを表現する。こうした確率 P (B|A) を尤度と呼ぶ。事象 A の確率分布 P (A) の
推定値は、観測データ B を得る前と得た後で変化する。B を得る前の推定値 P (A) を事前確率と呼び、B を
得た後の推定値 P (A|B) を事後確率と呼ぶ。単語 wn と文字列 D が、クラス Ck を原因として生成されたと
考えれば、P (wn |Ck ) と P (D|Ck ) は尤度に、P (Ck |D) は事後確率に相当する。従って、単純ベイズ分類器は
事後確率最大化法の応用と見なせる。そろそろ Data クラスの実装に入ろう。可読性のため、訓練データの
ラベルを表現する Label クラスを別に宣言する。
bayes.scala
case class Label(id: String)
case class Data(words: Seq[String], label: Label)
分類器の本体は NaiveBayes クラスとして定義する。確率は mutable.Map や mutable.Set により保持する。
23
無線部
Scala で実装する機械学習入門
開発班
bayes.scala
import scala.collection.mutable.{Map, Set}
class NaiveBayes(samples: Seq[Data]) {
val labels = Set[Label]()
val vocab = Set[String]()
val prior = Map[Label , Int]()
val denom = Map[Label , Int]()
val numer = Map[(String , Label), Int]()
}
vocab は訓練データ全体に含まれる異なる単語の集合である。prior は事前分布 P (C) を数えるカウンタで
ある。denom はクラス毎の単語数を単語の重複を許して数える。numer は単語とクラスの結合事象について
出現回数を数える。条件付き確率 P (w|C) は以下の pwc メソッドで計算する。
bayes.scala
def pwc(w: String , c: Label) = numer.get((w, c)) match {
case Some(num) => (num + 1.0) / (denom(c) + vocab.size)
case None
=>
1.0 / (denom(c) + vocab.size)
}
未知の単語 w が出現した場合、P (w|C) = 0 であり、式 (5.5) の計算ができない。その場合でも、他の既知の
単語からクラス C を推定できる可能性が残る。そこで、式 (5.6) に示すように、出現頻度に下駄を履かせて
P (w|C) > 0 となるように工夫している。これをラプラススムージングと呼ぶ。
P (w|C) =
numer(w ∩ C) + 1
numer(w ∩ C) + 1
∑
=
.
denom(C)
+ |vocab|
numer(w ∩ C) + 1
(5.6)
w∈vocab
式 (5.1) に従い P (D|C)P (C) を計算する pdc メソッドを実装する。数値計算の都合により、対数を利用する。
bayes.scala
def pdc(words: Seq[String], label: Label): Double = {
val pc = Math.log(prior(label).toDouble / samples.size)
return pc + words.map(w => Math.log(pwc(w, label))).sum
}
NaiveBayes クラスのコンストラクタに、学習の処理を追記する。単語やクラスの頻度を数えるだけである。
bayes.scala
samples.foreach{case Data(words , label) => {
labels += label
words.foreach(w => numer((w, label)) = 0)
}}
labels.foreach(prior(_) = 0)
labels.foreach(denom(_) = 0)
samples.foreach{case Data(words , label) => {
vocab ++= words
prior(label) += 1
words.foreach(w => denom(label) += 1)
words.foreach(w => numer((w, label)) += 1)
}}
最後に、式 (5.5) に基づき P (D|C)P (C) が最大になるクラス C を特定する predict メソッドを実装する。
24
無線部
Scala で実装する機械学習入門
開発班
bayes.scala
def predict(words: Seq[String]) = labels.maxBy(pdc(words , _))
以上で、単純ベイズ分類器が完成した。さっそく遊んでみよう。訓練するにも分類するにも文字列データが
必要だが、某百科事典の XML データが利用できる。文字列を単語に分解するには形態素解析ライブラリが
便利である。日本語向け形態素解析ライブラリは ChaSen や MeCab が著名だが、今回は Lucene GoSen1 を
利用する。100%Java なライブラリなので Scala での使い勝手が良い。形態素解析には辞書データが必要だが
IPA 辞書を内蔵する jar ファイルを選べば良い。jar ファイルにクラスパスを通して下記のコードを実行する。
bayes.scala
import net.java.sen.SenFactory
import net.java.sen.dictionary.Token
val tagger = SenFactory.getStringTagger(null)
val tokens = new java.util.ArrayList[Token]()
val source = Source.fromFile("something.txt")
tagger.analyze(source.mkString , tokens)
実行すると tokens に単語のリストが格納される。形態素の情報を取り出すには下記のように記述すれば良い。
bayes.scala
import collection.JavaConversions._
for (morph <- tokens.map(_.getMorpheme)) {
val ps = morph.getPartOfSpeech
val surface = morph.getSurface
}
百科事典の記事の XML データから本文を取り出すには scala.xml.XML クラスが便利である。
bayes.scala
val xml = XML.loadString(source.mkString)
val txt = (xml \\ "text")(0).text
記事を形態素解析にかけると、無意味な品詞が多く含まれていることに気付く。文章の内容を表現する上で
副詞や助詞、助動詞はほとんど役に立たず、むしろノイズになる。副詞や助詞、助動詞など文法的な機能を
果たすものの語彙的な意味を持たない語を機能語と呼ぶ。反対に、名詞や形容詞、動詞など語彙的な意味を
持つ語を内容語と呼ぶ。情報検索やテキストマイニングでは、機能語を除去して探索空間を小さくするのは
常套手段である。今回は、固有名詞だけを内容語とし、Table 5.1 (a) の記事を学習させて、Table 5.1 (b) の
所在地を推測させてみた。霧島山や鳥海山など複数県に跨がる山がどの県に分類されるかで遊ぶと面白い。
Table 5.1: classification example.
(a) sample data.
(b) test result.
山形県
328 語
赤城山
143 語
群馬県
新潟県
424 語
天城山
89 語
静岡県
群馬県
679 語
榛名山
64 語
群馬県
静岡県
592 語
霧島山
127 語
宮崎県
宮崎県
525 語
妙高山
57 語
新潟県
1 http://github.com/lucene-gosen/lucene-gosen
25
第 6 章 最尤推定法
6.1
混合分布モデル
データマイニングには、単純パーセプトロンやサポートベクターマシンのような線型分類器を用いる手法の
他に、統計的手法も有用である。データが正規分布やポアソン分布などの確率分布に従うと仮定し、平均や
分散といった母数を学習する。学習した確率分布を利用して、例えば、過去の確率分布と現在の確率分布を
比較することで、異常な状態を検出できる。第 6 章で扱う混合分布モデルなら、データが観測された原因を
推論することもできる。混合分布モデルでは、確率分布の線型和により分布が与えられる。Fig. 6.1 に示す
混合正規分布モデルがその代表例である。式 (6.1) の通り、K 個の正規分布の線型和が確率分布を与える。
Fig. 6.1: Gaussian mixture model.
P (x) =
K
∑
wk N (x|µk , Sk ),
k=1
K
∑
wk = 1.
(6.1)
k=1
ただし、µ は平均で、S は分散共分散行列である。混合分布モデルを構成する個別の分布は、モデル内部の
何らかの状態の表出であると見ることができる。即ち、観測者からは見えないものの、内部に K 個の状態が
潜在し、重み wk で状態 k が選択され、N (x|µk , Sk ) に従うデータを生成する。このように、観測者からは
観測できない変数が存在し、確率分布を選択するようなモデルを潜在変数モデルと呼ぶ。
6.2
最尤推定の理論
混合分布モデルを用いるには、潜在状態 k の母数 θk の学習が必須である。式 (6.1) の混合正規分布モデルの
例では θk = (wk , µk , Sk ) を学習する。直観的には、式 (6.2) に示す通り、尤度 L を最大にする θ を学習する。
L(θ) =
N
∏
P (xn |θ).
(6.2)
n=1
第 5 章でも紹介したが、尤度 P (B|A) は、事象 A が原因で観測データ B が得られたという仮説を裏付ける。
観測データの独立同分布性を仮定し、xn の列で表現する場合、尤度は式 (6.2) に示す同時確率で定義される。
26
無線部
Scala で実装する機械学習入門
開発班
式 (6.1) に示した混合正規分布モデルの例で、具体的な定義を示すと、尤度 L は式 (6.3) のように定義する。
L(θ) =
N ∑
K
∏
wk N (xn |µk , Sk ).
(6.3)
n=1 k=1
尤も、数値計算の都合により、式 (6.4) に示すように対数で定義することの方が多く、対数尤度と呼称する。
L(θ) = log
N ∑
K
∏
N
∑
wk N (xn |µk , Sk ) =
n=1 k=1
log
n=1
K
∑
wk N (xn |µk , Sk ).
(6.4)
k=1
尤度 L(θ) を最大にする θ を探索することで確率分布の母数を求める手法を最尤推定と呼ぶ。微分積分学の
知識を用いて解ける筈だが、残念ながら、式 (6.4) を最大化する θ = {wk , µk , Sk } を解析的に求めることは
困難である。試しに、式 (6.4) の尤度関数 L(θ) を平均 µk で偏微分すると、式 (6.5) を得る。
N
K
N
∑
∑
∂L
∂ ∑
=
log
wk N (xn |µk , Sk ) =
γnk Sk−1 (xn − µk ).
∂µk
∂µk n=1
n=1
(6.5)
k=1
寄与率 γnk は xn が潜在状態 k により得られたとする仮説の事後確率 P (k|xn ) であり、式 (6.6) で与えられる。
P (k, xn )
wk N (xn |µk , Sk )
=∑
.
γnk = P (k|xn ) = ∑
wk N (xn |µk , Sk )
P (k, xn )
k
(6.6)
k
式 (6.5) の偏微分を 0 にする µk を、推定平均値 µˆk として計算する。
µˆk =
N
N
∑
1 ∑
γnk xn , Nk =
γnk .
Nk n=1
n=1
(6.7)
同様にして、分散共分散行列の推定値 Sˆk も計算できる。
N
1 ∑
ˆ
Sk =
γnk (xn − µk ) t(xn − µk ).
Nk n=1
(6.8)
重みの推定値 wˆk のみ、式 (6.1) の制約条件を満たす必要があるので、ラグランジュの未定乗数法で求める。
wˆk =
Nk
.
N
(6.9)
寄与率 γnk が求まれば、式 (6.7)(6.8)(6.9) も計算可能だが、そもそも γnk の計算に µk , Sk , wk が必要である。
そこで補助関数法を導入し、目的関数 L の下限たる補助関数 Q の最大化により、間接的に L を最大化する。
L(θ) = max Q(γ, θ).
γ
(6.10)
式 (6.11)(6.12) のように γ と θ を交互に修正しつつ反復的に L を最大化する。ただし、t は試行回数である。
γ t+1 = arg max Q(γ, θt ),
(6.11)
γ
θt+1 = arg max Q(γ t , θ).
(6.12)
θ
ここで、L(θt ) = Q(γ t , θt ) と定義すると、式 (6.11)(6.12) を交互に反復することにより、L は単調増加する。
L(θt+1 ) = Q(γ t+1 , θt+1 ) ≥ Q(γ t+1 , θt ) ≥ Q(γ t , θt ) = L(θt ).
(6.13)
以上の補助関数法を混合正規分布モデルに適用する。まず、イェンゼンの不等式を利用して Q を導出する。
∑
f (x) を定義域と値域が実領域の凸関数、pn を n pn = 1 となる正の実数列とすると、式 (6.14) が得られる。
(
)
∑
∑
pn f (xn ) ≥ f
pn xn .
(6.14)
n
n
27
無線部
Scala で実装する機械学習入門
開発班
log が凹関数であることに注意して、正負を入れ替えつつ式 (6.14) を式 (6.4) に適用すると、式 (6.15) を得る。
L(θ) ≥
N ∑
K
∑
γnk log
n=1 k=1
K
∑
wk N (xn |µk , Sk )
γnk = 1.
= Q(γ, θ),
γnk
(6.15)
k
後は Q を最大化する γ と θ̂ を求めるだけであるが、寄与率 γ はラグランジュの未定乗数法により式 (6.6) を、
推定母数 θ̂ は Q(θ) の偏微分により式 (6.7)(6.8)(6.9) を得る。補助関数を導入せずとも結局は同じ式となり、
補助関数法の意義がないように思われるかもしれないが、式 (6.6) と式 (6.7)(6.8)(6.9) を交互に反復すれば、
尤度 L(θ) が単調増加する。尤度 L(θ) の収束解を保証する上では重要な議論である。
6.3
最尤推定の実装
第 6.2 節に述べたように、最尤推定に補助関数法を組み合わせて収束解を得る手法を期待値最大化法と呼ぶ。
式 (6.6) で γnk を計算する E-射影と、式 (6.7)(6.8)(6.9) で µˆk , Sˆk , wˆk を計算する M-射影を交互に反復する。
なお、今回は分散共分散行列 Sk が対角行列だと仮定して、ベクトルで代用する。◦ はアダマール積である。
{N
}
∑
1
◦2
(6.16)
Sˆk =
γnk xn − µ◦2
k .
Nk
n=1
では、実装作業に入ろう。まず EM クラスを定義する。訓練データは 1 次元とし、潜在状態の数を K とする。
em.scala
class
val
val
val
val
}
EM(samples: Seq[Double], K: Int) {
P = Array.ofDim[Double](K, samples.size)
W = Array.ofDim[Double](K)
M = Array.ofDim[Double](K)
S = Array.ofDim[Double](K)
P は P (k, xn |θk )、W は wk 、M は µk 、S は Sk を格納する。続く normal メソッドは、N (x|µ, Sk ) を求める。
em.scala
def normal(x: Double , m: Double , s: Double) = {
val pow = 0.5 * (x - m) * (x - m) / s
denom / Math.sqrt(s) * Math.exp(- pow)
}
ただし、定数 denom は、下記のように宣言する。
em.scala
val denom = 1 / Math.sqrt(2 * Math.PI)
後はコンストラクタで期待値最大化の本体を実装するだけだが、その前に初期化を忘れずに記述しておこう。
em.scala
for (k <- 0 until K) {
W(k) = Math.random
M(k) = Math.random
S(k) = Math.random
}
val wsum = W.sum
for (k <- 0 until K) W(k) /= wsum
28
無線部
Scala で実装する機械学習入門
開発班
初期化を忘れずに記述したら、期待値最大化による学習ルーチンの本体を記述して EM クラスは完成である。
em.scala
for (step <- 1 to 100) {
for ((x, i) <- samples.zipWithIndex) {
for (k <- 0 until K) P(k)(i) = W(k) * normal(x, M(k), S(k))
val sum = (for (k <- 0 until K) yield P(k)(i)).sum
for (k <- 0 until K) P(k)(i) = P(k)(i) / sum
}
for (k <- 0 until K) {
val nk = P(k).sum
M(k) = (P(k) zip samples).map{case(p, x) => p * x * 1}.sum / nk
S(k) = (P(k) zip samples).map{case(p, x) => p * x * x}.sum / nk
S(k) = S(k) - M(k) * M(k)
}
}
以上で、最尤推定が完成した。動作確認を兼ねて、Table 6.1(a) の母数に従う訓練データを学習させてみた。
Table 6.1: estimation example.
(a) sample data.
(b) test result.
wk
µk
Sk
wk
µk
Sk
0.57454
0.89739
0.69749
0.56005
0.81515
0.57402
0.42546
0.19543
0.83826
0.43995
0.22432
0.82382
1,000 件の訓練データに対して学習ステップを 100 回繰り返した結果が Table 6.1(b) である。収束の速度は
初期化時の乱数生成に依存するため、回数を決め打ちするのではなく、尤度を式 (6.4) に従い計算し、尤度が
収束した時点で学習を打ち切ると良い。今回は主に紙面の都合により訓練データを 1 次元に限定して実装を
単純化したが、Table 6.1(b) は推定の良し悪しを視覚的に判断しにくい。そこで、多次元に対応するように
改良して、Fig. 6.2 に示す分布を学習させた。等高線は 1,000 個の訓練データに対し学習を 2 回繰り返した
結果である。学習結果をアニメーション表示して徐々にモデルが改善されていく様子を観察すると面白い。
3
0.08
2
15
0.
y
1
0.1
0.0
0
2
0.1
0
5
−1
0.
18
0.05
−2
−3
−3
−2
5
0.1
0.08
−1
0.03
0.120.10
0
x
1
2
3
Fig. 6.2: expectation maximization on Gaussian mixture model.
29
© Copyright 2024 ExpyDoc
