非線形計画ソルバーMINOSの実装とSVMへの適用

Posted at 2024-01-21

はじめに

機械学習の手法の一つであるSVM(サポートベクトルマシン)を実装するためには、二次計画問題を解く必要があります(詳しくは後述します）。SVMを、特に大規模データセットに適用するために、効率的に二次計画問題を解く手法がいくつか提案されていますが、その一つである分解法(Osuna et al.,1996)について調べようと思い、論文を読んでいたら、分解法の手法の中で大規模問題を分割した、部分問題を解くための手法として、非線形計画ソルバーMINOS(Murtagh et al., 1976)が紹介されていたので、まずはそこからと思い、MINOSの論文を調べpythonを用いて自分で実装してみました。この記事を最後まで読むのは、相当なもの好きしかいないかもしれませんが、少しでも興味のある方はどうぞお付き合いいただければ幸いです。（本記事を読むに当たり、KKT条件やシンプレックス法、準ニュートン法など最適化の最低限の知識はあった方がいいかと思います。)

[参考文献]

LARGE-SCALE LINEARLY CONSTRAINED OPTIMIZATION
PRML 7章

非線形計画問題

MINOSが扱う非線形計画問題は以下のように定式化されます。

\displaylines{
\text{minimize} \quad F(\textbf{x})=f(\textbf{x}) \\
\text{subject to} \quad A\textbf{x}=b, l\leqq \textbf{x} \leqq u \\
\text{where} \quad A \quad \text{is} \quad m \times n,\quad m \leqq n
}

$f(\textbf{x})$は二次関数を含む非線形関数で、制約等式は$m$個、変数の数は$n$個です。全ての変数には下限と上限があります（矩形制約）。(下限だけ、上限だけでも可）下限、上限は変数ごとに異なっていてもよいです。また、等式制約の個数は変数の個数よりも小さいものとします。

MINOSの基本的アプローチ

MINOSでは変数$\textbf{x}$を$m$個の$x_{B}$, $s$個の$x_{S}$,$n-m-s$個の$x_{N}$に分割します。$m$は等式制約の個数、$s$はこちらから与えるパラメータです。$x_{B}$をbasic変数、$x_{S}$をsuper-basic変数、$x_{N}$をnon-basic変数と呼ぶことにします。non-basic変数$x_{N}$は上限値、または下限値の値をとり、super-basic変数$x_{S}$はアルゴリズムに従って、等式制約で表される実行可能領域内を動き、残りのbasic変数$x_{B}$は$x_{S}, x_{N}$の値が定まった状態で、等式制約から一意に定まります。
変数の分割に応じて、行列$A$も$A= \lbrack B,S,N \rbrack$と分割します。

KKT条件

MINOSのアルゴリズムを導出するために、制約付き最適化問題の最適性のための必要条件である、KKT条件を求めます。まず、今考えている問題のラグランジュ関数は以下のようになります。

L(\textbf{x}, \mu, \lambda) = f(\textbf{x})-\mu^T(A\textbf{x}-b) + \lambda_{1}^T(\textbf{l}-\textbf{x}) + \lambda_{2}^T(\textbf{x}-u)

よってKKT条件は以下のようになります。

\begin{align}
& \nabla f(\textbf{x}) - A^T\mu - \lambda_{1} + \lambda_{2} = 0 \tag{1} \\
& \lambda_{1}^T(\textbf{l}-\textbf{x}) = 0 \tag{2} \\
& \lambda_{2}^T(\textbf{x}-\textbf{u}) = 0 \tag{3} \\
& \lambda_{1} \geqq 0, \lambda_{2} \geqq 0 \tag{4} \\
& A\textbf{x} = b \tag{5} \\
& \textbf{l} \leqq \textbf{x} \leqq \textbf{u} \tag{6}\\
\end{align}

式$(1)$はラグランジュ関数を最適化変数$x$で偏微分したものです。式$(2), (3)$は相補性条件、式$(4)$は不等式制約に対するラグランジュ乗数は正であるという条件、式$(5), (6)$は元の問題の制約条件です。これらを合わせたものがKKT条件となります。KKT条件は最適性のための必要条件なので、問題の最適解を求めるために、KKT条件を満たす$x^* $,$\mu^* $, $\lambda_{1}^* $,$\lambda_{2}^* $を求めることを目標とします。

MINOSの基本等式の導出

MINOSのアルゴリズムの導出の前に、KKT条件の式$(1)$を少し変形しておきます。$\lambda_{1} - \lambda_{2} = \lambda'$とおきます。$x_{B}$と$x_{S}$は $\textbf{l} < \boldsymbol{x}<\textbf{u}$ を満たすため、相補性条件$(2), (3)$ より,　$\lambda_{1}, \lambda_{2}$の上から$m+s$個の要素は$0$となります。また、$x_{Ni}=l$の場合は$\lambda_{1}$の$m+s+i$個目が、$x_{Ni}=u$の場合は$\lambda_{2}$の$m+s+i$個目が$0$以上の値をもちます。相補性条件も考慮すると、$\lambda_{1}$と$\lambda_{2}$は同じインデックスで、$0$よりも大きい値を同時に持つことはないことが分かります。
以上のことから

\lambda' = \begin{pmatrix} 0 \\0\\:\\:\\0\\\lambda\end{pmatrix}

$\lambda \in \mathbb{R}^{N}$と表されます。

ここで、$x_{Ni}=l$を満たす$i$の集合を$L$, $x_{Ni}=u$を満たす$i$の集合を$U$とすると$i \in L$のとき$\lambda_i > 0$, $i \in U$のとき$\lambda_i < 0$が成り立ちます。

MINOSのアルゴリズムの導出

MINOSではまず、等式制約$A\textbf{x}=b$を満たす初期解$x_0$から始め、$x \rightarrow x+\alpha p_x$と更新していくことで、上記のKKT条件を満たす$x$を求めます。$\alpha$は更新幅で、$p_x$は更新方向です。MINOSのアルゴリズムを導くために、一旦$\alpha$を無視し、$x+p_x$がKKT条件の$(1)$式と$(5)$式を満たすと仮定します。まず$(5)$式を満たすとすると、

\displaylines{
A(x+p_x) = b \\
Ap_x = 0 \\
\begin{pmatrix} B & S & N \end{pmatrix} \begin{pmatrix} p_B \\ p_S \\ p_N \end{pmatrix} = 0
}

となります。ここで$x$は$Ax=b$を満たすことを用いました。
また、$x_N$は変数の制約の境界値に固定されているため、$p_N=0$となることを合わせると以下の式が得られます。

Property 1. \quad
\begin{pmatrix} B & S & N \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} p_B \\ p_S \\ p_N \end{pmatrix} = 0

$Property1$は$x+p_x$が実行可能領域に存在することを保証します。
次に$x+p_x$が式$(1)$を満たすとし、テイラー展開を用いると、以下のように式変形できます。

\displaylines{
\nabla f(x+p_x) - A^T\mu - \lambda' = 0 \\
\sim \nabla (f(x)+\nabla f(x)^T p_x )  - A^T\mu - \lambda' = 0 \\
\nabla f(x) + \nabla^2 f(x) p_x = A^T\mu + \lambda'
}

$\nabla f(x) = g, \nabla^2f(x) = G$とおくと以下の式が得られます。

Property2. \quad \begin{pmatrix} g_B \\ g_S \\ g_N \end{pmatrix} +G\begin{pmatrix} p_B \\ p_S \\ p_N \end{pmatrix} = \begin{pmatrix} B^T & 0 \\ S^T & 0 \\ N^T & I \end{pmatrix} \begin{pmatrix} \mu \\ \lambda \end{pmatrix}

さて、$Property1$からは$Bp_B = -Sp_S$ が得られます。また、最適解の周辺では$\lVert p_s \lVert \rightarrow 0$となることを用いると、$Property2$の一行目、三行目からはそれぞれ、$B^T \pi = g_B$, $\lambda = g_N - N^T \pi $が、また二行目からは

\displaylines{
Z^TGZp_s = -h \\
(h=Z^Tg = g_s -S^T\pi) \\
(Z=\begin{pmatrix} -B^{-1}S \\ I \\ 0 \end{pmatrix} )
}

が得られます。
$m \times m $の正方行列$B$は$Bx=y, B^Tx=y$の形の方程式を解いたり、$B^{-1}$を直接求めたりする必要があるため、$LU$分解しておきます。また、$G$を直接計算するには手間がかかるため、$s \times s$の上三角行列$R$を用いて、$R^TR \simeq Z^TGZ$と近似します。

MINOSのアルゴリズム

最適解においては$\lVert p_s \lVert \rightarrow 0$となるため、$Property2$の二行目から導かれる式から$\lVert h \lVert \rightarrow 0$となります。$\lVert p_s \rVert \simeq 0$即ち、$\lVert h \lVert \simeq 0$が満たされている状態ではKKT条件の$(1), (5)$式が満たされています。$(2),(3)$は$\lambda'$の定義から満たされ、また$(6)$を満たすように$\alpha$を求め$x$を更新していきます。よってあとは$(4)$が満たされればKKT条件を満足し、その時の$x$が最適解ということになります。$(4)$を満たすかどうかは、以下のアルゴリズムの$step2$で調べることになります。

アルゴリズム

$prepare$

$s$を定める
$A\textbf{x}=b$を満たす$\textbf{x}_{0}$を求める($x_N=$$l$ or $u$,$l < x_B, x_S < u$)
$calculate \quad g$
$LU \quad decomposition$ $of$ $B$
$set$ $R=I$
$solve$ $B^T\pi=g_B$
$compute$ $h=g_S-S^T\pi$

$step1$

$\lVert h \lVert \simeq 0$でなければ$step3$へ

$step2$

$compute$ $\lambda = g_N-N^T\pi$
$\lambda_{q1}=min(\lambda_i),(i\in L)$, $\lambda_{q2}=max(\lambda_i),(i\in U)$とおいた時、$\lambda_{q1}<-\epsilon$ $or$ $\lambda_{q2}>\epsilon$ を満たさなければ終了、$x$が最適解
$\lambda_{q1} < -\epsilon$ $or$ $\lambda_{q2} > \epsilon$　なら
(1)$|\lambda_{q}| = max(|\lambda_{q1}|, |\lambda_{q2}|)$
(2)$N_q$を$S$の列に追加
(3)$\lambda_{q}$を$h$に追加
(4)$R$に新たな列を追加
(5)$s$を１増やす

$step3$

solve $R^TRp_s = -h$
solve $LUp_B = -Sp_s$
$p_x^T = (p_B^T \quad p_S^T \quad 0)$

$step4$

$x+\alpha p_x$が実行可能な領域に収まる範囲内で最大の$\alpha_{max}$を見つける
もし$\alpha_{max}=0$なら$step7$へ

\newcommand{\min}{\mathop{\rm min}\limits}

$step5$(直線探索)

$f(x+\alpha p_x) = \min_{0<\theta \leqq \alpha_{max}} f(x+\theta p_x)$
$x \rightarrow x+ \alpha p_x$
$calculate$ $g$

$step6$

$solve$ $B^T\pi = g_B$
$compute$ $\overline{h} = g_S - S^T\pi$
$update$ $R$
$set$ $h=\overline{h}$
もし $\alpha < \alpha_{max}$なら $step1$ へ

$step7$
( $\alpha=\alpha_{max}$ $or$ $\alpha_{max}=0$ )

$detect$ $p$ ($x_p$ $is$ $on$ $the$ $bounds$, $0<p\leqq m+s$)

$(a)$ $if$ $0<p\leqq m$

\begin{pmatrix} B \\ x_B^T \end{pmatrix} \quad \begin{pmatrix} S \\ x_S^T \end{pmatrix}

$interchange$ $pth$ $and$ $qth$ $column$ ($q$の探し方については後述します)
$modify$ $L$, $U$, $R$ $and$ $\pi$
$compute$ $h=g_S - S^T\pi$
$go$ $to$ $(b)$

$if$ $m<p\leqq m+s$, $\quad$ $q=p-m$
$(b)$ $qth$ $variable$ を$nonbasic$にする

$delete$ $qth$ $column$

\begin{pmatrix} S \\ x_S^T \end{pmatrix} and \begin{pmatrix} R \\ h^T \end{pmatrix}

$restore$ $R$ $to$ $trianglar$ $form$
$decrease$ $s$ $by$ $1$ $and$ $go$ $to$ $step1$

step2 Rに加える列について

$step2$で条件を満たした場合$nonbasic$変数から$superbasic$変数へ移行させるため、$R$に列を加える必要があります。ここではどんな列を加えればいいのかについて書きます。
まず$R^TR$が近似する行列は、$superbasic$変数を一つ増やした場合次のように変化します。
まず$Z \rightarrow \overline{Z}$となり、$\overline{Z}$は

\displaylines{
\overline{Z} = \begin{pmatrix} Z & z \end{pmatrix} \\
where \quad z= \begin{pmatrix} -B^{-1}a_q \\ e_s \\ 0 \end{pmatrix}
}

と表されます。
よって近似する行列は

\displaylines{
\overline{Z}^TG\overline{Z} = \begin{pmatrix} Z^T \\ z^T \end{pmatrix}G\begin{pmatrix} Z & z \end{pmatrix} \\
=\begin{pmatrix} Z^TGZ & Z^TGz \\ z^TGZ & z^TGz \end{pmatrix} \tag{7}
}

となります。
また

\overline{R} = \begin{pmatrix} R & r \\ 0 & \rho \end{pmatrix}

とおくと

\displaylines{
\overline{R}^T\overline{R} = \begin{pmatrix} R^T & 0 \\ r^T & \rho \end{pmatrix} \begin{pmatrix} R & r \\ 0 & \rho \end{pmatrix} \\
= \begin{pmatrix} R^TR & R^Tr \\ r^R & r^Tr+\rho^2 \end{pmatrix} \tag{8}
}

となるので式$(7)$と式$(8)$を見比べるのですが、ここでベクトル$v$を導入します。

\displaylines{
v = \frac{g(x+\delta z) - g(x)}{\delta} \\
= \frac{\nabla f(x+\delta z) - \nabla f(x)}{\delta} \\
\simeq Gz
}

よって式$(7)=(8)$が成り立つためには、$v \simeq Gz$であることを利用すると、以下の二つの式が成り立つ必要があることが分かります。
まず

R^Tr = Z^Tv

また

\displaylines{
r^Tr+\rho^2 = z^TGz \\
\Leftrightarrow \rho^2=z^Tv-\lVert r \lVert^2 = \sigma \\
\therefore \rho = \sigma^{\frac{1}{2}}
}

この二つの式から$r$と$\rho$が求まります。
また、$R$の対角成分の$max$を$d_{max}$,
$min$を$d_{min}$とすると、$\lVert v \lVert$が$d_{max}$よりひどく大きいか、$d_{min}$よりひどく小さい時、また$\sigma$が負の時

\displaylines{
R = \begin{pmatrix} R & 0 \\ 0 & \rho \end{pmatrix} \\
\rho = (z^Tv)^{\frac{1}{2}} \quad or \quad 1.0
}

とします。
また、$h$の新たな要素として$\lambda_q$を加えられる理由については、$h$と$\lambda$の定義式を比べれば分かると思います。

step5 直線探索(linear search)について

$step5$では$0<\alpha \leqq \alpha_{max}$の範囲の中で$f(x+\alpha p_x)$が最小となる$\alpha$を求めます。これは$\Phi (\alpha) = f(x+\alpha p_x)$とおくと、一変数の最適化問題となり、直線探索と呼ばれます。直線探索を行うにはアルミホ($Almijo$)条件やウルフ($Wolfe$)条件を満たす$\alpha$を求めれば良いです。今回はウルフ条件を用いました。ウルフ条件とは以下の二つの条件のことを指します。

\displaylines{
1. \quad \Phi(\alpha) \leqq \Phi(0)+C_1\alpha \Phi'(0) \\
2. \quad \Phi'(\alpha) \geqq C_2\Phi'(0)\\
 0 < C_1 < C_2 < 1
}

この二つの式が意味するところは、$\alpha$を横軸に取ったグラフを書いてみると分かりやすいと思います。これらのウルフ条件を満たす$\alpha$が存在するためには、$\Phi'(0) < 0$を満たす必要があります。
また$\Phi'(\alpha)$は

\displaylines{
\Phi'(\alpha) = \frac{df(x+\alpha p_x)}{d\alpha} \\
= \sum_{i} \frac{\partial f(x')}{\partial x_i'} \frac{dx_i'}{d\alpha} \\
= \nabla f(x')^Tp_x
}

と表され、$g$を用いて計算することができます。
ウルフ条件を満たす$\alpha$はバックトラック法を用いて求めることができます。

バックトラック法

$0<\rho<1$を満たす$\rho$を決める。$\alpha=\alpha_{max}$とおく。
$\alpha$がウルフ条件の二式を満たせば$\alpha$を出力
$\alpha = \rho \alpha $とする。

step6 Rの更新について

$step6$における$R$の更新は$R^TR$の形で行われます。$R$の更新式は以下の通りです。

\displaylines{
\overline{R}^T\overline{R} = R^TR + \frac{1}{\alpha y^Tp_s}yy^T + \frac{1}{h^Tp_s}hh^T \\
(y = \overline{h} - h)
}

step7-(a) について

$step7(a)$では$basic$変数と$superbasic$変数を交換します。$basic$変数$x_{Bp}$と交換する$superbasic$変数のインデックス$q$の求め方と、行列$B$,$S$が変化することに伴う、$\pi$, $R$の更新の方法について書きます。
まず$q$の決め方についてです。

\displaylines{
B^T\pi_p = e_p \\
y=S^T\pi_p \\
y_{max} = max|y_j| \\
d_j = min \lbrace |x_j-l_j|, |x_j-u_j| \rbrace, (0<j<s) \\
d_q = max \lbrace d_j |  |y_j| \geqq 0.1y_{max} \rbrace \\
}

$e_p$は単位ベクトルです。
$\pi$の更新、$R$の更新はそれぞれ以下で表されます。

\displaylines{
\overline{\pi} = \pi + (\overline{h_q}/y_q)\pi_p \\
\overline{R}^T\overline{R} = (I+ve_q^T)R^TR(I+e_qv^T) \\
(v = -\frac{1}{y_q}(y+e_q))
}

step7-(b)について

$superbasic$変数を$nonbasic$変数に移すことにより、$s$が一つ小さくなるため$R$の$q$番目の列を取り除きますが、$R$は上三角行列に保持する必要があります。そのためには単位行列の$(i,i), (i,j), (j,i), (j,j)$番目の要素をそれぞれ$cos(\theta), sin(\theta), sin(\theta), -cos(\theta)$に変更した行列$Q_{i,j}$を用いて以下のように変形をします。

Q_{s-1,1}\cdots Q_{q,q+1} \begin{bmatrix} R \quad with \\ qth \quad column \\ deleted \end{bmatrix} = \begin{bmatrix} \overline{R} \\ 0 \end{bmatrix}

このように変形してもよいのは、$q$番目の列を取り除いた後の行列を$R'$とすると、$R'^TR'= \overline{R}^T\overline{R}$が満たされるからです。
$R$の更新は$R^TR$の形で行われるため、$R^TR$が同じであれば、$R$にどのような変更を加えても良いということです。また、$R$は$R^TR$を分解して得られます。

SVMの二次計画問題

これまで説明してきた、非線形計画問題の最適化手法MINOSをSVMに適用するために、二値分類SVMについて簡単に紹介したいと思います。
SVMでは二つのクラスの分類の際、特徴量空間において、識別境界から最短距離にあるデータ点までの距離(マージン）が最大となるように識別境界を決定します。
マージンの中にもデータ点が入り込むことを許すソフトマージンSVMでは、マージン内に入り込んだそれぞれのデータ点に対して、ペナルティを課し、そのペナルティが最小に、マージンが最大になるように目的関数を設計し最適化します。最適化の際には、主問題をそのまま扱うのではなく、ラグランジュ関数をKKT条件を用いて変形し、最適化変数を変更した双対問題を最適化します。ソフトマージンSVMの最適化問題（双対問題）は以下のように定式化されます。

\displaylines{
minimize \quad L(\textbf{a}) = -\sum_{n} a_n +\frac{1}{2}\sum_{n}\sum_{m} a_n a_m t_n t_m \phi^T(x_n) \phi(x_m)\\
\leftrightarrow W(\Gamma) = -\Gamma^TI+\frac{1}{2}\Gamma^T D\Gamma, \quad (D_{i,j} = t_it_jk(x_i, x_j)) \\ 
s.t. \quad 
\left\{ \begin{array}{l}
\Gamma^T\textbf{t}=0 \\
0 \leqq \Gamma \leqq C
\end{array}
\right.
}

（ここでは分量の都合上これ以上詳しい話は書きませんが、もっと詳しく知りたい場合は、PRMLの７章などをご覧ください。）
$a_i$はサンプルごとに与えられるラグランジュ乗数で、$t_i$はサンプルのクラスラベル（$t_i \in \lbrace 1,-1 \rbrace$）です。また、クラスの識別関数$y$は$\textbf{a}$を用いて計算されます。
上記の定式化から、SVMの学習が二次最適化問題であることが分かったので、上で説明したMINOSを適用することができます。
$LU$分解、$B$の逆行列の計算、$R^TR$の分解、一次方程式の解の計算、インデックスの管理、MINOSのアルゴリズム本体などのコードについては、以下に置いておきます。
https://colab.research.google.com/drive/1LE7uGNfEellFnxjShg4_qel2ks0-eAxw?usp=sharing!