「サポート・ベクタ・マシーン ( Support Vector Machine )」は、「Vladimir Vapnik」と「Alexey Chervonenkis」によって 1963 年に発表されたクラス分類手法です。当初は線形分離しかできませんでしたが、「カーネル法 ( Kernel Method )」との組み合わせによる非線形な分類が 1992 年に提案されて以降、非常に注目される機械学習モデルの一つとなりました。ここでは、サポート・ベクタ・マシーンの概要やカーネル法、具体的な実装方法の「逐次最小問題最適化法 ( Sequential Minimal Optimization ; SMO )」について紹介します。
「線形識別モデル」は特徴ベクトルを超平面でクラスに分類する手法でした。クラスが二つだった場合、超平面は線形識別関数
を使って g(x) = 0 で表すことができ、g(x) の正負によってクラスを判別することができるのでした。ここで w は重みベクトル、w0 はバイアス・パラメータになります。
線形識別モデルの一つとして「パーセプトロン」を以前紹介しました。パーセプトロンは、線形分離可能であるときは必ず収束することが保証されています。しかし、得られた超平面が最適であるとは限りません。平面の引き方は無数にありえるので、どのような面が得られるかは重みベクトルとバイアス・パラメータの初期値や学習パターンの入力順番などによって変わります。そこで、データが正しく分離できているという条件下で分割面とデータの間の距離を調べ、その中で最短の距離 ( この距離のことを「マージン ( Margin )」といいます ) が最大になったときを最適な分割面とみなします。このとき、二つのクラスの中の少なくともそれぞれ一つのデータがマージンを持つことになります。
ある学習パターン ( x1, x2, ... xN ) が、超平面 g(x) = 0 によって線形分離可能であると仮定します。また、それぞれの学習パターンに対する識別関数の出力値が正の場合、その教師データは +1 であり、逆に負の場合は -1 になっているとします。超平面 g(x) = 0 と i 番目の入力ベクトル xi の間の距離は g(xi) / ||w|| で与えられるのでした ( *1-1 )。g(xi) は正負のいずれの場合もありますが、対応する教師データ ti を掛けることで必ず正になります。よって、マージンは
で表されます。但し、min は得られた全ての値の最小値を返す関数とします。この値を最大にするような w, w0 を求めればいいのですが、tig(xi) / ||w|| は w, w0 によって変化するので、この問題をそのまま解くことは非常に困難です。そこで、w, w0 を適当な値で定数倍しても tig(xi) / ||w|| は変化しないことを利用して、マージンを持つ点 xi に対して
が成り立つように w, w0 を調整します。このとき、正しく分離されていれば全てのデータに対して
が成り立つことに注意して下さい。マージンは 1 / ||w|| となるのでこれを最大化する問題になり、さらにこれは ||w||2 を最小化する問題と同じなので、結局
を制約条件 ti( wTxi + w0 ) ≥ 1 下で最小化すればよいことになります。但し、係数の 1 / 2 は後の計算を単純化するために導入したものです。最小化する値にバイアス・パラメータ w0 は存在しませんが、制約条件の中にあるため w から間接的に得られることに注意して下さい。この問題は「二次計画法 ( Quadratic Programming ; QP )」と呼ばれる、線形不等式による制約条件下での二次関数の最適化問題です。この問題を解くためには、ラグランジュの未定乗数法を利用する必要があります ( 補足 1 )。制約条件 ti( wTxi + w0 ) ≥ 1 ごとにラグランジュ乗数 ai ≥ 0 を定義すると、最適化したい関数 F( w, w0, a ) は
となります。但し、a = ( a1, a2, ... aN )T とします。制約条件 ti( wTxi + w0 ) ≥ 1 より、二項目は減算とすることで制約条件領域内の最大化方向に ai が向かうことに注意して下さい。このとき、
∇wF = w - Σi{1→N}( aitixi )
∂F / ∂w0 = Σi{1→N}( aiti )
となります。但し、∇wF = ( ∂F / ∂w1, ∂F / ∂w2, ... ∂F / ∂wM ) とします。極値をとる点においては ∇wF = 0, ∂F / ∂w0 = 0 となるので、
w = Σi{1→N}( aitixi )
Σi{1→N}( aiti ) = 0
という条件が得られ、これらを F( w, w0, a ) に代入すると
F( a ) | = | (1/2)Σi{1→N}( Σj{1→N}( aiajtitjxiTxj ) ) - Σi{1→N}( aitiΣj{1→N}( ajtjxiTxj ) + aitiw0 - ai ) |
= | Σi{1→N}( ai ) - (1/2)Σi{1→N}( Σj{1→N}( aiajtitjxiTxj ) ) |
となって、w と w0 を消去することができます。これは元の式の「双対表現 ( Dual Representation )」と呼ばれ、重みベクトルとバイアス・パラメータがラグランジュ乗数に置き換わり、学習パターンの内積で表現された式になります。制約条件は
ai ≥ 0
Σi{1→N}( aiti ) = 0
となって、この問題は再び二次計画法となります。F( a ) を 最大化する ai が得られれば、g(x) は w を ai で置き換えることによって
g(x) | = | wTx + w0 |
= | Σi{1→N}( aitixiTx ) + w0 |
で計算でき、その正負によってクラスの判別が可能になります。
F( a ) の「Karush-Kuhn-Tucker 条件 ( KKT 条件 )」 ( 補足 1 ) は、
ai ≥ 0
tig(xi) ≥ 1
ai[ tig(xi) - 1 ] = 0
なので、全ての学習パターンに対して ai = 0 または tig(xi) = 1 が成り立ちます。tig(xi) = 1 となる学習パターンはマージンを持つデータを意味し、「サポート・ベクトル ( Support Vector )」と呼ばれます。それ以外のデータは ai = 0 なので g(x) の計算に対して何も影響せず、学習が完了した後はサポート・ベクトル以外の学習パターンは予測値の計算に必要ないため、処理時間を短くすることができます。
ai が得られた後は、バイアス・パラメータ w0 を求める必要があります。任意のサポート・ベクトルは tig(xi) = 1 を満たすことから、あるサポート・ベクトル xj を使って
と表せるので、tj を両辺に掛けて ( tj2 = 1 が常に成り立つことに注意して ) 整理すると
という結果が得られます。但し、i ∈ S とはサポート・ベクトルを和の対象にすることを表しています。数値計算上の誤差を考慮して、実際には全てのサポート・ベクトルを使い
と計算することが一般的なようです。ここで NS はサポート・ベクトルの数を表します。
今までは、学習パターンが線形分離可能であることを前提に話を進めてきましたが、実際にはそれぞれのクラスの分布が重なっている場合もあり得るので、それを考慮したモデルのほうが適切である可能性もあります。そのため、誤分類に対してはある程度のペナルティを与えるだけでその存在は許容するような式を検討します。分割面とデータとの距離がマージンより小さくなる場合、tig(xi) < 1 となるので、「スラック変数 ( Slack Variable)」ξi を導入してそのようなデータに対しては
になるようにします。なお、スラック変数とは、不等式制約を等式制約に変換するための変数を意味します。もし、データが分割面からはみ出さず正しく分類されていたのであれば、0 ≤ tig(xi) < 1 なので 0 < ξi ≤ 1 が成り立ちます。しかし、誤分類していた場合、tig(xi) < 0 なので ξi > 1 となります。分割面とデータとの距離がマージン以上であれば ξi = 0 としておけば影響しないので、結局、制約条件は
ξi ≥ 0
tig(xi) + ξi ≥ 1
とすることができます。このようにして得られたマージンは「ソフト・マージン ( Soft Margin )」と呼ばれ、対して先に述べたマージンは「ハード・マージン ( Hard Margin )」になります。全データに対して ξi = 0 が成り立てば、両者は等しくなるのは明らかです。
ξi はできるだけ小さい方がいいので、最小化する問題を以下のようにします。
但し、C > 0 とします。C はペナルティに対する影響度を制御するパラメータで、大きくなるほどマージンより小さい距離を持つデータを許容しなくなり、C → ∞ の極限においてはハード・マージンと等しくなります。このとき最適化したい関数は、
となります。但し、ai ≥ 0, μi ≥ 0 はラグランジュ乗数で、ξ = ( ξ1, ξ2, ... ξN )T, a = ( a1, a2, ... aN )T, μ = ( μ1, μ2, ... μN )T をそれぞれ表します。KKT 条件 ( 補足 1 ) は、
ai ≥ 0
tig(xi) + ξi - 1 ≥ 0
ai[ tig(xi) + ξi - 1 ] = 0
μi ≥ 0
ξi ≥ 0
μiξi = 0
で、w, w0, ξi による偏微分から極値をとる点は
∇wL = w - Σi{1→N}( aitixi ) = 0
∂L / ∂w0 = Σi{1→N}( aiti ) = 0
∂L / ∂ξi = C - ai - μi = 0
なので、F( w, w0, ξ, a, μ ) の双対表現は
F( a ) | = | (1/2)Σi{1→N}( Σj{1→N}( aiajtitjxiTxj ) ) + CΣi{1→N}( ξi ) - Σi{1→N}( Σj{1→N}( aiajtitjxiTxj ) + aiξi - ai ) - Σi{1→N}( μiξi ) |
= | Σi{1→N}( ai ) - (1/2)Σi{1→N}( Σj{1→N}( aiajtitjxiTxj ) ) + Σi{1→N}( ( C - ai - μi )ξi ) | |
= | Σi{1→N}( ai ) - (1/2)Σi{1→N}( Σj{1→N}( aiajtitjxiTxj ) ) |
となって、ハード・マージンの場合と同じになります。但し、ai = C - μi ≤ C なので、ai の制約条件は
0 ≤ ai ≤ C
Σi{1→N}( aiti ) = 0
となります。これを「矩形制約 ( Box Constraint )」といいます。ハード・マージンの場合と同様に、ai = 0 の場合は g(x) の計算に何の影響も及ぼしません。ai < C のとき、μi = C - ai > 0 なので ξi = 0 が成り立ち、分割面との距離がマージンに等しいことになります。また、ai = C ならば μi = 0 であり、先の議論から ξi ≤ 1 のときは正しく分類されているのに対し、ξi > 1 ならそのデータは誤分類されていることになります。
*1-1) 詳細は、「パターン認識 (4) パーセプトロン ( Perceptron )」の「1) 最近傍決定則 ( NN 法 ) ( Nearest Neighbor Rule ; NN Rule )」の最後の箇所を参照
サポート・ベクタ・マシーンやパーセプトロンは線形分離が原則でしたが、学習パターンを変換することで非線形な境界にも対応することは原理的に可能です。例えば対数やべき乗へ変換することで線形分離な形にすることができるかもしれません。以下に示す一次元データは線形分離することができませんが、10 を中心とした距離で判別ができることに気付けば f(x) = ( x - 10 )2 と変換することで線形分離可能となります。下図の下側は二次元空間 f(x) = ( x, ( x - 10 )2 ) にデータを写像した結果を示しています。
一般的には、学習パターン x = ( x1, x2, ... xK ) を非線形関数 f(x) = ( f1(x), f2(x), ... fM(x) ) に変換した上で処理を行うことになりますが、問題はどのように変換すれば線形分離可能になるのか事前にはわからないという点です。データが高次元になればグラフで確認して最適な変換方法を見つけるようなこともできないので、とりあえず考えられる変換を全て要素として列挙して試すような形になり、その数が膨大になれば計算時間も飛躍的に増大することになります。
サポート・ベクタ・マシーンの双対表現は
となるのでした。x を f(x) に置き換えると
となります。ここで
とさらに置き換えて
とします。この関数 K( xi, xj ) を「カーネル関数 ( Kernel Function )」といいます。f(x) のサイズが非常に大きいとしても、その内積が他の形で表されるのなら処理時間は短縮できることが期待できます。実際、多項式カーネル
は展開すると、多項定理より [ M! / Πk{1→K}( rk! )rc! ]crcΠk{1→K}( xikrkxjkrk ) の和で表されます。但し、和の各項は Σk{1→K}( rk ) + rc = M を満たす全ての組み合わせから成ります。したがって、[ M! / Πk{1→K}( rk! )rc! ]crc = A として A1/2Πk{1→K}( xikrk ), A1/2Πk{1→K}( xjkrk ) を要素とするベクトルの内積に置き換えることができます。また、内積がどのような形になるのか知らなくてもカーネル関数さえ分かっていれば計算することが可能で、さらには無限次元の内積についても処理できるようになります。これを「カーネル・トリック ( Kernel Trick )」または「カーネル置換 ( Kernel Substitution )」といいます。他によく利用されるカーネル関数として、以下に示すガウス・カーネルがあります。
ガウス・カーネルは無限次元の内積となるカーネル関数の例になります。
カーネル関数のサンプル・プログラムを以下に示します。
/* カーネル法基底クラス */ struct Kernel_Base { // カーネル関数の初期化 virtual void init() = 0; // カーネル関数による内積計算 // x1, x2 はベクトルの各要素を表す virtual void operator()( double x1, double x2 ) = 0; // カーネル関数による内積計算後の後始末 virtual void cleanUp() = 0; // カーネル関数による計算結果を返す virtual double val() const = 0; }; /* 線形カーネル */ class Kernel_Linear : public Kernel_Base { double val_; // 内積計算結果を保持する変数 public: // カーネル関数の初期化 virtual void init() { val_ = 0; } // カーネル関数による内積計算 virtual void operator()( double x1, double x2 ) { val_ += x1 * x2; } // カーネル関数による内積計算後の後始末 virtual void cleanUp() {} // カーネル関数による計算結果を返す virtual double val() const { return( val_ ); } }; /* 多項式カーネル */ class Kernel_Polynomial : public Kernel_Base { double p_; // 多項式の指数 double c_; // 多項式の定数 double val_; // 内積計算結果を保持する変数 public: // 多項式の指数を指定して構築 // // p 多項式の指数 // c 多項式の定数 Kernel_Polynomial( double p, double c ) : p_( p ), c_( c ) { assert( p > 0 ); } // カーネル関数の初期化 virtual void init() { val_ = c_; } // カーネル関数による内積計算 virtual void operator()( double x1, double x2 ) { val_ += x1 * x2; } // カーネル関数による内積計算後の後始末 virtual void cleanUp() { val_ = std::pow( val_, p_ ); } // カーネル関数による計算結果を返す virtual double val() const { return( val_ ); } }; /* ガウシアンカーネル */ class Kernel_Gaussian : public Kernel_Base { double s_; // 標準偏差 double val_; // 内積計算結果を保持する変数 public: // ガウス関数のパラメータ(標準偏差) s を指定して構築 Kernel_Gaussian( double s ) : s_( s ) { assert( s > 0 ); } // カーネル関数の初期化 virtual void init() { val_ = 0; } // カーネル関数による内積計算 virtual void operator()( double x1, double x2 ) { val_ += std::pow( x1 - x2, 2 ); } // カーネル関数による内積計算後の後始末 virtual void cleanUp() { val_ = std::exp( -val_ / ( 2 * s_ * s_ ) ); } // カーネル関数による計算結果を返す virtual double val() const { return( val_ ); } }; /* カーネル関数の計算 kernel カーネル関数への参照 s1, e1 配列1の範囲 s2 配列2の開始位置 */ template< class In1, class In2 > typename In1::value_type Kernel( Kernel_Base& kernel, In1 s1, In1 e1, In2 s2 ) { kernel.init(); for ( ; s1 != e1 ; ++s1, ++s2 ) kernel( *s1, *s2 ); kernel.cleanUp(); return( kernel.val() ); }
カーネル関数を計算するためのクラスは 3 つの処理部分 init, operator(), cleanUp から成ります。init は値を初期化する部分で、operator() が反復処理を実行する部分、最後の cleanUp で後始末を行います。例えば多項式カーネルなら init で定数 c を保持し、operator() で xik と xjk の積を加算していきます。ここまでで Σk{1→K}( xikxjk ) + c の値が得られるので、最後の cleanUp で M 乗して計算を完了します。このクラスを使って一連の処理を行う関数が Kernel です。
サポート・ベクタ・マシーンは二次計画法を解くことによって解を得ることができます。一般的な二次計画法のアルゴリズムの場合、未知数と同程度のサイズの行列演算が必要になりますが、学習パターンが多くなれば処理時間やメモリ使用量などが増大するため実用的ではなくなります。そこで、サポート・ベクタ・マシーン専用に考案されたアルゴリズムがいくつかあります。今回は、その中から最もよく利用されている手法の一つである「逐次最小問題最適化法 ( Sequential Minimal Optimization ; SMO )」を紹介します。SMO は、一回のステップで最適化するラグランジュ乗数を最小にすることで行列演算などを回避して処理を単純化し、高速化を実現しています。サポート・ベクタ・マシーンには制約条件として Σi{1→N}( aiti ) = 0 があるため、最適化するラグランジュ乗数の最小値は 2 となります。このとき、あるラグランジュ乗数の値が増減すれば、制約条件に従ってもう一方のラグランジュ乗数が増減することになります。
サポート・ベクタ・マシーンの双対表現に対し、二つのラグランジュ乗数 ai, aj 以外の値は全て定数であると仮定します。つまり、最適化時に二つのラグランジュ乗数のみを変数として扱うことになります。簡略化のため Kij = K( xi, xj ) で表すとして、N = 3, i = 1, j = 2 のときに
F( a1, a2 ) | = | a1 + a2 - (1/2)Σk{1→3}( aka1tkt1Kk1 + aka2tkt2Kk2 + aka3tkt3Kk3 ) |
= | a1 + a2 + a3 - (1/2)( | |
a1a1t1t1K11 + a1a2t1t2K12 + a1a3t1t3K13 + | ||
a2a1t2t1K21 + a2a2t2t2K22 + a2a3t2t3K23 + | ||
a3a1t3t1K31 + a3a2t3t2K32 + a3a3t3t3K33 ) | ||
= | a1 + a2 + a3 - (1/2)a12K11 - (1/2)a22K22 - (1/2)a32K33 - a1a2t1t2K12 - a2a3t2t3K23 - a3a1t3t1K31 |
になることから類推すると ( このとき、ti2 = 1、Kij = Kji となることに注意してください )、一般式は
F( ai, aj ) | = | ai + aj - (1/2)ai2Kii - (1/2)aj2Kjj - aiajtitjKij - aitiΣk{1→N;k≠i,j}( aktkKik ) - ajtjΣk{1→N;k≠i,j}( aktkKjk ) + Const. |
= | ai + aj - (1/2)ai2Kii - (1/2)aj2Kjj - aiajtitjKij - aitivi - ajtjvj + Const. |
となります。但し Const. は ai, aj を含まない定数項です。また、途中から
vi = Σk{1→N;k≠i,j}( aktkKik )
vj = Σk{1→N;k≠i,j}( aktkKjk )
としています。制約条件から
となりますが、tj2 = 1 となることを利用して両辺に tj を掛け、s = titj と表せば
とすることができます。ここで γ は定数です。この式を使って aj を消去することができて
と表すことができるので、これを ai で微分して
dF / dai | = | 1 - s - aiKii + ( γ - sai )sKjj - ( γ - sai )sKij + aiKij - tivi + stjvj |
= | ai( 2Kij - Kii - Kjj ) + γs( Kjj - Kij ) - tivi + titj2vj + 1 - s | |
= | ai( 2Kij - Kii - Kjj ) + γs( Kjj - Kij ) + ti( vj - vi ) + 1 - s |
になります。ここで、
であることを利用すると
vi = g(xi) - w0 - aitiKii - ajtjKij
vj = g(xj) - w0 - aitiKij - ajtjKjj
となります。但し、ここに出現した g(x), w0 や ai, aj は前の反復処理で得られた値を意味します。これを ti( vj - vi ) に代入すると
ti( vj - vi ) | = | ti[ g(xj) - g(xi) - aitiKij - ajtjKjj + aitiKii + ajtjKij ] |
= | ti[ g(xj) - g(xi) ] - aiKij - ajsKjj + aiKii + ajsKij |
また、γ も前の反復処理で得られた値だとして γs( Kjj - Kij ) を計算すると
γs( Kjj - Kij ) | = | ( sai + aj )( sKjj - sKij ) |
= | aiKjj - aiKij + ajsKjj - ajsKij |
なので
と求めることができます。ここで、ai には前の値であることを示す "old" を付加しています。最後に、
なので、極値をとる点、すなわち dF / dai = 0 のときは
となります。
とおくと、η は d2F / dai2 に等しくなり、
とすれば、Ei は予測値と実測値との差異を表します。これらを使って
という漸化式が得られます。
漸化式から得られた ai は矩形範囲 0 ≤ ai ≤ C を満たす必要があります。もう一つの制約条件 sai + aj = γ から、aj = 0 のとき ai = sγ、aj = C のとき ai = sγ - sC なので、矩形範囲も考慮すると、ai の取りうる範囲 [ L, H ] は ti = tj のとき
ti ≠ tj のとき
となります。γ = saiold + ajold で計算すれば γ の値を得ることができるので、これを使って極値をとる ai が境界内に存在するか確認することができます。もし、ai ≥ H ならば ai = H とし、ai ≤ L ならば ai = L とします。最後に、γ が定数であることから saiold + ajold = sai + aj であることを利用すれば
から aj が得られます。
極値をとる点において F は最大値をとるはずなので、η は負にならなければなりません。しかし、学習パターンが一致していた場合など、そうならないときも発生します。その場合、最大値をとるのは境界であるはずなので、境界 L, H での F の値を直接計算して両者を比較し、より大きくなった方を採用します。
w0 は、
を利用して計算します。vi は反復処理前後において値が変わらないので
が成り立ちます。ここで "old" が付加された値は反復処理前の値であることを表します。この式を整理すると
となります。更新された ai が 0 < ai < C の範囲にあった場合、ξi = 0 が成り立つので、tig(xi) + ξi = 1 より g(xi) = ti となり、上式は
と変形することができます。同様に、aj を使う場合は
を使えば計算ができます。もし、ai, aj のどちらも矩形の境界上にあった場合は、両者から計算した w0 の値の範囲は KKT 条件を満たすので、二値の平均値を採用します ( 補足 2 )。
SMO では二つのラグランジュ乗数のみを一度に最適化するのでした。そこで、どのラグランジュ乗数を選択するかが問題になります。まず、一つ目のラグランジュ乗数を選択するために、KKT 条件を利用します。ラグランジュ乗数は KKT 条件を満たさなければならないので、そうでないラグランジュ乗数は最適化の候補とします。具体的には、ai = 0 のとき tig(xi) ≥ 1、ai = C のとき tig(xi) = 1 - ξi ≤ 1 ( μi = 0 より ξi ≥ 0 が成り立つので ) となるので、
ai > 0 かつ tig(xi) > 1
ai < C かつ tig(xi) < 1
のいずれかとなるようなラグランジュ乗数 ai を探します。処理開始時、候補は全ての学習パターンの中から選択されますが、一度全データに対してチェックが完了したら、次の処理では高速化のため 0 でも C でもない ( つまり境界上にない ) ラグランジュ乗数だけを対象にします。そして、境界上にないラグランジュ乗数が一つも更新されない状態になったら、再び全データを対象に候補を選び出します。もし、全データに対して一つも更新されるラグランジュ乗数が存在しなくなったら処理が完了となります。
二つ目のラグランジュ乗数は、境界上にないものの中から ai の更新値が最大になる値を選びます。更新値は前述の漸化式から ti( Ej - Ei ) / η となりますが、η にはカーネル関数が含まれているので高速化のため | Ej - Ei | で評価します。もし、更新がされなかった場合は、次に境界上にないものをランダムな位置から順番に選択し、それでも更新されなければ境界上にあるラグランジュ乗数をランダムな位置から順番に選択して処理を行います。よって、更新がされない限り最終的には全てのデータが評価されることになります。
高速化のため、Ei はキャッシュ化されます。更新された ai が境界上になければそれはサポート・ベクトルであることを意味するので Ei = 0 となります ( w0 の計算で、更新された g(xi) が ti に等しいことを前提にしていることから成り立つことに注意してください )。ai, aj 以外のラグランジュ乗数 am ( m ≠ i, j ) は、
vm | = | Σk{1→N;k≠i,j}( aktkKmk ) |
= | g(xk) - w0 - aitiKmi - ajtjKmj |
より vm が反復処理前後で変化しないことを利用して
となるので、この式から
という結果が得られ、この値を Em に加算することで更新ができます。なお、Em が初期化されるのは am が境界上にない場合のみなので、この更新は am が境界上にない場合のみについて行います。当然、キャッシュを利用するのはラグランジュ乗数が境界上にない場合だけに限られます。
SMO のサンプル・プログラムを以下に示します。
/* SMO : 逐次最小問題最適化法 ( SMO ) による SVM の解の計算用クラス */ class SMO { const Matrix< double >& x_; // 訓練データ std::vector< double > t_; // 訓練データの解 Kernel_Base& kernel_; // カーネル関数 double c_; // ソフト・マージン用パラメータ(大きいほど誤分類へのペナルティが増す) double eps_; // ラグランジュ乗数評価時の余裕値 double tolerance_; // KKT条件評価時の余裕値 double w0_; // バイアス・パラメータ std::vector< double > a_; // ラグランジュ乗数 std::vector< std::vector< double >::size_type > sv_index_; // サポートベクトル(0でないa[])のインデックス std::map< size_t, double > err_cache_; // エラー値のキャッシュ // ラグランジュ乗数の範囲チェック bool isNotBound( double a ) { return( a > eps_ && a < ( c_ - eps_ ) ); } bool isBound( double a ) { return( ! isNotBound( a ) ); } // a_[i] の更新評価 ( KKT 条件のチェック ) bool checkKKT( bool alldata, size_t i, double* ei ); // a_[i] との更新ペア a_[j] を探し、見つかったら更新処理を行う bool innerLoop( size_t i, double ei ); // 目的関数の値を計算する double getObjFunc( std::vector< double >::size_type i, std::vector< double >::size_type j, double ai_new, double aj_new, double kii, double kij, double kjj ); // a_[i],a_[j]を更新する bool update( size_t i, size_t j, double si, double ej ); // 学習パターンから求められる理論値を計算する double f( size_t i ) const; public: // SMO コンストラクタ : 学習パターン x と教師データ t を入力して訓練を行う SMO( const Matrix< double >& x, const std::vector< bool >& t, Kernel_Base& kernel, double c, double eps, double tolerance, unsigned int maxCnt ); // output : 入力値から理論値を計算する template< class In > double output( In s ) const; }; /* SMOコンストラクタ : 学習パターン x と教師データ t を入力して訓練を行う kernel : カーネル関数 c : ソフト・マージン用パラメータ(大きいほど誤分類へのペナルティが増す) eps : ラグランジュ乗数評価時の余裕値 tolerance : KKT条件評価時の余裕値 maxCnt : 更新処理の最大回数 */ SMO::SMO( const Matrix< double >& x, const vector< bool >& t, Kernel_Base& kernel, double c, double eps, double tolerance, unsigned int maxCnt ) : x_( x ), kernel_( kernel ), c_( c ), eps_( eps ), tolerance_( tolerance ), w0_( 0 ) { ErrLib::CheckLinearModel( x_, t.begin(), t.end() ); size_t xSize = x.rows(); // 学習パターンの数 // 正例と負例を数値変換する for ( vector< bool >::const_iterator cit = t.begin() ; cit != t.end() ; ++cit ) t_.push_back( ( *cit ) ? 1 : -1 ); // ラグランジュ乗数の初期化 a_.assign( xSize, 0.0 ); // 外側のループ処理 bool alldata = true; // 全データ対象に処理するなら true bool changed; // a[i] が更新されたら true double ei; // i 番目の論理値と実測値の誤差 for ( unsigned int loop = 0 ; loop < maxCnt ; ++loop ) { changed = false; for ( size_t i = 0; i < xSize; ++i ) { if ( ! checkKKT( alldata, i, &ei ) ) continue; if ( innerLoop( i, ei ) ) changed = true; } if ( alldata ) { alldata = false; if ( ! changed ) break; } else { if ( ! changed ) alldata = true; } } // 有効なラグランジュ乗数のインデックスを取得 for ( size_t i = 0 ; i < xSize ; ++i ) { if ( a_[i] > eps_ ) { sv_index_.push_back( i ); } } } /* SMO::checkKKT : a[i] の更新評価 ( KKT 条件のチェック ) 同時に Ei を取得する */ bool SMO::checkKKT( bool alldata, size_t i, double* ei ) { if ( ( ! alldata ) && isBound( a_[i] ) ) return( false ); // 教師データの理論値との誤差を求める if ( isNotBound( a_[i] ) ) { *ei = err_cache_[i]; // 0 < a < C ならキャッシュを代入 } else { *ei = f( i ) - t_[i]; } double tFi = *ei * t_[i]; // ti・f(i) - 1 // KKT条件のチェック return( ( a_[i] < ( c_ - eps_ ) && tFi < -tolerance_ ) || ( a_[i] > eps_ && tFi > tolerance_ ) ); } /* SVM::innerLoop : a[i] との更新ペア a[j] を探し、見つかったら更新処理を行う */ bool SMO::innerLoop( size_t i, double ei ) { size_t xSize = x_.rows(); // 0 < aj < C を満たし | Ej - Ei | が最大のものを選択 double maxDelta = 0.0; size_t max_j = xSize; for ( size_t j = 0 ; j < xSize ; ++j ) { if ( isBound( a_[j] ) ) continue; double ej = err_cache_[j]; if ( std::abs( ej - ei ) > maxDelta ) { maxDelta = std::abs( ej - ei ); max_j = j; } } if ( max_j < xSize ) { if ( update( i, max_j, ei, err_cache_[max_j] ) ) { return( true ); } } // 0 < aj < C を満たすものをランダムに選択 size_t offset = (size_t)( ( (double)rand() / (double)RAND_MAX ) * (double)( xSize - 1 ) ); for ( size_t j = 0 ; j < xSize ; ++j ) { size_t pos = ( j + offset ) % xSize; if ( isNotBound( a_[pos] ) ){ if ( update( i, pos, ei, err_cache_[pos] ) ) { return( true ); } } } // aj = 0, C のものをランダムに選択 offset = (size_t)( ( (double)rand() / (double)RAND_MAX ) * (double)( xSize - 1 ) ); for ( size_t j = 0 ; j < xSize ; ++j ) { size_t pos = ( j + offset ) % xSize; if ( isBound( a_[pos] ) ) { if ( update( i, pos, ei, f( pos ) - t_[pos] ) ) { return( true ); } } } return( false ); } /* SMO::getObjFunc : 目的関数の値を計算する 計算結果に a_[i], a_[j] 以外のラグランジュ係数成分から成る定数項は含まないことに注意 */ double SMO::getObjFunc( vector< double >::size_type i, vector< double >::size_type j, double ai_new, double aj_new, double kii, double kij, double kjj ) { double vj = f( j ) - w0_ - t_[j] * a_[j] * kjj - t_[i] * a_[i] * kij; double vi = f( i ) - w0_ - t_[j] * a_[j] * kij - t_[i] * a_[i] * kii; return( aj_new + ai_new - kjj * aj_new * aj_new / 2.0 - kii * ai_new * ai_new / 2.0 - t_[j] * t_[i] * kij * aj_new * ai_new - t_[j] * aj_new * vj - t_[i] * ai_new * vi ); } /* SMO::update : a[i],a[j]を更新する */ bool SMO::update( size_t i, size_t j, double ei, double ej ) { if ( i == j ) return( false ); // 線形制約対角線の端点 double boundL, boundH; if ( t_[i] != t_[j] ) { boundL = std::max( 0.0, a_[i] - a_[j] ); boundH = std::min( c_, a_[i] - a_[j] + c_ ); } else { boundL = std::max( 0.0, a_[i] + a_[j] - c_ ); boundH = std::min( c_, a_[i] + a_[j] ); } if ( boundL == boundH ) return( false ); // 目的関数の二階導関数の値 ( 符号は逆転 ) double kii = Kernel( kernel_, x_[i], x_[i].end(), x_[i] ); double kjj = Kernel( kernel_, x_[j], x_[j].end(), x_[j] ); double kij = Kernel( kernel_, x_[i], x_[i].end(), x_[j] ); double eta = 2.0 * kij - kii - kjj; // a[i] の更新 double ai_new; if ( eta < 0.0 ) { ai_new = a_[i] - ( t_[i] * ( ej - ei ) / eta ); if ( ai_new > boundH ) { ai_new = boundH; } else if ( ai_new < boundL ) { ai_new = boundL; } } else { // a[i] = boundL のときの目的関数の値 double Lobj = getObjFunc( i, j, boundL, a_[j] + t_[i] * t_[j] * ( a_[i] - boundL ), kii, kij, kjj ); // a[i] = boundH のときの目的関数の値 double Hobj = getObjFunc( i, j, boundH, a_[j] + t_[i] * t_[j] * ( a_[i] - boundH ), kii, kij, kjj ); if ( Lobj > Hobj + eps_ ) { ai_new = boundL; } else if ( Lobj + eps_ < Hobj ) { ai_new = boundH; } else { return( false ); } } if ( std::abs( ai_new - a_[i] ) < eps_ * ( ai_new + a_[i] + eps_ ) ) { return( false ); } // a[j] の更新 double aj_new = a_[j] + t_[i] * t_[j] * ( a_[i] - ai_new ); // w0 更新 double w0_old = w0_; if ( isNotBound( ai_new ) ) { w0_ -= ei + ( ai_new - a_[i] ) * t_[i] * kii + ( aj_new - a_[j] ) * t_[j] * kij; } else if ( isNotBound( aj_new ) ) { w0_ -= ej + ( ai_new - a_[i] ) * t_[i] * kij + ( aj_new - a_[j] ) * t_[j] * kjj; } else { w0_ -= ( ei + ( ai_new - a_[i] ) * t_[i] * kii + ( aj_new - a_[j] ) * t_[j] * kij + ej + ( ai_new - a_[i] ) * t_[i] * kij + ( aj_new - a_[j] ) * t_[j] * kjj ) / 2.0; } // err_cache の更新 ( i, j 以外 ) for ( std::map< size_t, double >::iterator it = err_cache_.begin() ; it != err_cache_.end() ; ++it ) { size_t m = it->first; if ( m == i || m == j ) { continue; } else if ( isNotBound( a_[m] ) ) { it->second += t_[j] * ( aj_new - a_[j] ) * Kernel( kernel_, x_[j], x_[j].end(), x_[m] ) + t_[i] * ( ai_new - a_[i] ) * Kernel( kernel_, x_[i], x_[i].end(), x_[m] ) + w0_ - w0_old; } } a_[i] = ai_new; a_[j] = aj_new; // err_cache の更新 ( i, j ) if ( isNotBound( a_[i] ) ) err_cache_[i] = 0.0; if ( isNotBound( a_[j] ) ) err_cache_[j] = 0.0; return( true ); } /* SMO::f : 学習パターンから求められる理論値を計算する */ double SMO::f( size_t i ) const { double sum = 0; for ( size_t j = 0 ; j < x_.rows() ; ++j ) { if ( a_[j] < eps_ ) continue; sum += a_[j] * t_[j] * Kernel( kernel_, x_[j], x_[j].end(), x_[i] ); } sum += w0_; return( sum ); } /* SMO::output : 入力値から求められる理論値を計算する s : 入力値の開始位置 */ template< class In > double SMO::output( In s ) const { double sum = 0; for( std::vector< std::vector< double >::size_type >::const_iterator it = sv_index_.begin() ; it != sv_index_.end() ; ++it ) { sum += a_[*it] * t_[*it] * Kernel( kernel_, x_[*it], x_[*it].end(), s ); } sum += w0_; return( sum ); }
サンプル・プログラムは前述したアルゴリズムをそのまま実装した形になります。行列を表すクラス Matrix は実装していませんが、x[i] で i 番目の行の開始を持つ反復子を返し、x[i].end() で行の終了の次の位置を得ることができます。また、メンバ関数 rows でその行列の行数を返します。コンストラクタが訓練を行う箇所で、ここでの処理によってラグランジュ乗数が決定します。処理の最初にある ErrLib::CheckLinearModel は、学習パターン x と教師データ t のデータ数が一致するかどうかをチェックするための関数で、サンプル・プログラムの中では実装はされていません。コンストラクタ内にあるループ処理では、最初に alldata を true にして処理を開始します。alldata = true のときは、全てのラグランジュ乗数が評価対象になります。その後、alldata = true のときは false に変更した上で更新があったかをチェックし、なければ処理を終了します。また、alldata = false のときは更新がない場合のみ alldata を true にします。これにより、一度全てのラグランジュ乗数を評価した後は、更新がされなくなるまで境界上にないラグランジュ乗数のみが更新対象になります。更新するかどうかの判定は checkKKT で行い、ここで KKT 条件を満たしていないかどうかをチェックします。もし、そうであれば innerLoop へ処理を移し、ペアとなるラグランジュ乗数を探索して update で更新処理を行います。
ラグランジュ乗数が境界上にあるかどうかの判定は厳密には行わず、ある程度の余裕を持ちます。この余裕値がコンストラクタで渡される変数 eps です。同様に、変数 tolerance は tig(xi) - 1 に対する余裕値になります。余裕値は、参考文献ではだいたい 10-2 から 10-3 くらいが適しているとしています。
サンプル・プログラムを使い、有名な「アヤメの品種データ ( iris flower data set )」で学習を行った結果を以下に示します。アヤメの品種データは 3 品種に対して 4 つのデータを持ちますが、今回は 2 クラスの分類になるため、品種は "Iris-setosa" と "Iris-versicolor" の二種類とし、データも "Sepal Length ( がく片の長さ )" と "Sepal Width ( がく片の幅 )" の 2 つとしています。処理は、ペナルティに対する影響度 C を 1000、ラグランジュ乗数と tig(xi) - 1 に対するしきい値への余裕度 eps, tolerance を 0.001 として行いました。まずは線形カーネルを使って学習した結果です。線形カーネルはいわゆる通常の内積です。よって、データは線形分離されることになります。
各学習パターンを表すプロットの中で、×の付いたものはサポート・ベクトルを表しています。データは線形分離可能であるため、境界線によって誤分類なくきれいに分割されていることがわかります。また、二つのクラスの分布に対し、ちょうど中央に境界線がある様子を確認することができます。
次は多項式カーネルを使った結果です。多項式の次数は 2、定数項はゼロとしました。
境界面は非線形となりますが、元々線形分離可能なデータだったためほぼ直線に近い状態となりました。
最後に、ガウスカーネルを使った結果です。標準偏差は 1 としました。
境界面によって、"Iris-versicolor" 側の領域は囲まれたような形になっています。線形分離可能なデータに対してこれでは正しく判定されません。このデータにガウスカーネルを使うのはやり過ぎなように思います。
今回は学習パターンのサイズが 2 だったので図によって確認することができましたが、サイズがもっと大きくなればそれは困難になります。どのカーネルが最も適切に分離をしてくれるのか、その時のパラメータはいくつがいいのか、学習とテストを繰り返して試行錯誤することは避けられません。
パーセプトロンが公開されたのは 1962 年、サポート・ベクタ・マシーンはその翌年に発表されたことになります。当初は線形分離しかできないということで日の目を見なかったわけですが、カーネル法との組み合わせで非線形な分類にも対応可能なことが 1992 年に発表されたことで一躍脚光を浴びるようになりました。これからも、過去にあった手法がまた注目を浴びるようなことはあるのかもしれません。
「ラグランジュの未定乗数法 ( Lagrange Multiplier )」については以前にも紹介したことがありますが、そのときは制約条件が等式でした。ここでは、制約条件が不等式になった場合を考えてみます。
制約条件 g(x) ≤ 0 において f(x) を最大化したいとします。このとき、あり得る場合として二つの可能性が考えられます。まず、求めたい最大値が g(x) < 0 内にある場合で、このときは制約条件は存在しないのと同様なので、∇f(x) = 0 から x を求めることができます。この場合を「無効制約 ( Inactive Constraint )」といいます。もう一つの可能性は求めたい最大値が g(x) = 0 上にある場合で、このときは制約条件が等式であるときの f(x) を解けばよいことになります。この場合は「有効制約 ( Active Constraint )」といいます。制約が有効な場合、ある定数 λ があって
が成り立つのでした。ここで λ = 0 になれば、無効制約と等しくなるので、制約条件が等式の場合と同様に
の極値を計算することでいいことになります。今回は、制約条件が g(x) ≤ 0 なので ∇g(x) は領域外の方向へ向いており、f(x) の最大値が制約条件の領域外ならば ∇f(x) も同じく領域外を向くことになるので λ > 0 となることに注意してください。上式の場合、制約が有効な場合と無効な場合を含めると以下の条件が成り立ちます。
g(x) ≤ 0
λ ≥ 0
λg(x) = 0
これを「Karush-Kuhn-Tucker条件 ( KKT条件 )」といいます。もし求めたい値が最小値だった場合、λ ≥ 0 を満たすためには F( x, λ ) の第二項の前の符号は逆転することに注意が必要です。
次に、制約条件が複数ある場合を考えます。x = ( x1, x2, ... xN )T としておきます。制約条件が全て等式 gi(x) = 0 ( i = 1, 2, ... M ) の場合、それぞれの制約条件は N 次元空間での曲面であり、それらが交差した面が存在するならば、極値を求めたい関数 f(x) が極値 c をとる面はその交差面に接しているはずです。もし接していなければ、交差した位置から交差面に沿って値が小さくなる、あるいは大きくなることができるからです。各制約条件による曲面の法線ベクトル ∇gi = ( ∂gi / ∂x1, ∂gi / ∂x2, ... ∂gi / ∂xN )T は交差面に対して垂直です。このとき、交差面上の f(x) = c の接点における法線ベクトル ∇f = ( ∂f / ∂x1, ∂f / ∂x2, ... ∂f / ∂xN )T は ∇gi の線型結合で表さなければなりません。すなわちある定数 λi が存在して
となります。つまり、∇f は ∇gi から成る M 次元の部分空間に含まれるということを意味します。もしそうならなければ、∇f は部分空間に直交する成分を持つことになります。その成分は交差面に接する方向を持つので、その方向に f の値が増減することを意味し、極値とはならないことになります。
N = 3, M = 2 のとき、g1(x), g2(x) は 3 次元空間上の曲面で、その交差した部分は曲線になります。∇g1, ∇g2 は交線に直交しており、f(x) = c が交線に接しているならば ∇f も交線に直交していることになるので、これら 3 つの法線ベクトルは同一平面上にあることになります。もし ∇f が ∇g1, ∇g2 を含む平面上になければ、その平面に直交する成分を ∇f が持つことになり、その分 ∇f は "傾いて" いることになります。つまり、f(x) = c が交線と交差していることを意味し、その位置では極値をとりません。
制約条件が不等式である場合も含めてさらに一般化すると、gi(x) = 0 ( i = 1, 2, ... M ) と hj(x) ≤ 0 ( j = 1, 2, ... K ) という制約条件下で f(x) を最大化したい場合は、ラグランジュ乗数を λ = ( λ1, λ2, ... λM )T, μ = ( μ1, μ2, ... μK )T として
を最適化すればよいことになります。但し、不等式制約に対しては以下の KKT 条件を満たす必要があります。
hj(x) ≤ 0
μj ≥ 0
μjhj(x) = 0
バイアス・パラメータは、0 < ai < 1 が満たされれば g(xi) = ti となるので、Ei を使って求めることができました。0 < aj < 1 も成り立っていれば、両者から計算した結果は等しくなります。しかし、ai, aj のいずれも矩形制約の境界上にある場合はこの二つの値は異なるものになります。
ai, aj のいずれも矩形制約の境界上にある場合、線形制約 sai + aj = γ は矩形制約の対角線上を通らなければなりません。このとき、s = 1 ならば、ai + aj = C であり、( ai, aj ) = ( C, 0 ) または ( 0, C ) です。また、s = -1 ならば -ai + aj = 0 なので、( ai, aj ) = ( 0, 0 ) または ( C, C ) となります。
ai = 0, ti = 1 の場合、KKT 条件から g(xi) ≥ ti > 0 が成り立ちます。よって、
より、Ei を使って得られた w0 を wi で表すと
という結果が得られます。このとき、aj の取りうる値は 0 か C のいずれかで、aj = 0 のときは s = -1 すなわち tj = -1 なので、KKT 条件から g(xj) ≤ tj < 0 であり、
より、Ej を使って得られた w0 を wj で表すと
となります。また、aj = C のときは s = 1 すなわち tj = 1 なので、KKT 条件から g(xj) = tj - ξj ≤ tj であり、上記と同じ結果が得られます。ai = 0, ti = -1 の場合からスタートしたときは、wi と wj の大小関係は逆転しますが、いずれにしても w0 は wi と wj の間に存在することになります。逆に、wi と wj の間に存在する値は全て KKT 条件を満たすことになります ( この範囲から外れると g(xi) と ti、または g(xj) と tj の大小関係が逆転するので条件を満たさなくなります )。よって、文献では二つの間の平均値を使うことで w0 の近似値としています。
前に戻る | タイトルに戻る |