はじめに

シリーズ第3弾で今回は制約付き非線形計画問題を解く主双対内点法を実装した*1。前回の準ニュートン法の記事はこちら。

yamagensakam.hatenablog.com

制約付き非線形計画問題を式で表すと、以下のような最適化問題になる。

$\displaystyle \begin{align} \min_{{\bf x} \in {\mathbb R}^n} \ & f({\bf x}) \\ {\rm s.t. } \ & g_i({\bf x}) \leq 0 \ (i = 1, \cdots, l) \\ & h_i({\bf x}) = 0 \ (i = 1, \cdots, m) \\ \end{align} \tag{1}$

なお今回取り上げる主双対内点法は、すぐ下に記すKKT条件が最適性の必要条件となることが保証されなければ成立しない。これが保証されるために制約条件が満たすべき条件が存在する。そのような制約条件に関する条件を制約想定と言い、1次独立条件やSlater条件などの種類がある。制約想定に関してはこちらが参考になる*2。

qiita.com

手法

内点法は最適性の1次の必要条件であるKKT条件を満たす点を実行可能領域の内部を通りながら探す手法である。まずKKT条件について簡単に触れた後、主双対内点法の考え方、アルゴリズムを述べる。

ラグランジュ関数とKKT条件

式（1）の最適化問題に対して、以下のような関数を考える。

$\displaystyle L({\bf x}, {\bf u}, {\bf v}) = f({\bf x}) + \sum_{i = 1}^{l} u_i g_i({\bf x}) + \sum_{i = 1}^{m} v_i h_i({\bf x}) \tag{2}$

これはラグランジュ関数と呼ばれる関数で、式（1）の最適化問題が上述の制約想定を満たす場合、この問題における1次の最適性の必要条件はラグランジュ関数を使って、

$\displaystyle \begin{align} & \nabla_{{\bf x}} L({\bf x}^{*}, {\bf u}^{*}, {\bf v}^{*}) = \nabla f({\bf x}^{*}) + \sum_{i = 1}^{l} u_i^{*} \nabla g_i({\bf x}^{*}) + \sum_{i = 1}^{m} v_i^{*} \nabla h_i({\bf x}^{*}) = {\bf 0} \\ & g_i({\bf x}^{*}) \leq 0 , \ \ \ u_i^{*} \geq 0, \ \ \ u_i^{*} g_i({\bf x}^{*}) = 0 \ \ \ (i = 1, \cdots, l) \\ & h_i({\bf x}^{*}) = 0 \ \ \ (i = 1, \cdots, m) \\ \end{align} \tag{3}$

と表すことができる。これはKKT条件と呼ばれ、主双対内点法ではこのKKT条件を満たす $({\bf x}^{*}, {\bf u}^{*}, {\bf v}^{*})$ （に十分に近い点）を探索することになる。

スラック変数とバリア関数の導入

ここで式（1）にスラック変数とバリア関数というものを導入した以下のような最適化問題を考える。

$\displaystyle \begin{align} \min_{{\bf x}, {\bf s}} \ & f({\bf x}) - \rho \sum_{i = 1}^{l} \log s_i\\ {\rm s.t. } \ & g_i({\bf x}) + s_i = 0, \ \ \ s_i \geq 0 \ \ \ (i = 1, \cdots, l) \\ & h_i({\bf x}) = 0 \ \ \ (i = 1, \cdots, m) \\ \end{align} \tag{4}$

$s_i$ がスラック変数でこの変数の導入により不等式制約は等式制約に置き換わる。また、 $- \sum_{i = 1}^{l} \log s_i$ がバリア関数で実行可能領域の内部では有効な値を取り、 $\rho > 0$ であれば制約の境界である $s_i = 0$ に近づくにつれ無限に発散していく。従って、バリア関数は不等式制約を満たさないことに対するペナルティとして働く。 $\rho \rightarrow 0$ で式（4）は元の最適化問題式（1）と等価になる。そのため、主双対内点法では式（4）の停留点方向に移動しながら徐々に $\rho$ を $0$ に近づけていくという操作を繰り返すことで式（1）の停留点を見つける。

主双対内点法のアルゴリズム

主双対内点法では現時点が $({\bf x}^{(k)}, {\bf s}^{(k)}, {\bf u}^{(k)}, {\bf v}^{(k)})$ だとして、 $({\bf x}^{(k + 1)}, {\bf s}^{(k + 1)}, {\bf u}^{(k + 1)}, {\bf v}^{(k + 1)}) = ({\bf x}^{(k)} + \alpha {\bf \Delta x}, {\bf s}^{(k)} + \alpha {\bf \Delta s}, {\bf u}^{(k)} + \alpha {\bf \Delta u}, {\bf v}^{(k)} + \alpha {\bf \Delta v})$ と更新する。以下、更新方向、ステップ幅、 $\rho$ をどのように更新するか述べる。

更新方向の計算

式（4）の停留点 ${\bar{\bf x}}, {\bar {\bf s}}$ ではKKT条件より、

$\displaystyle \begin{align} & \nabla f({\bar{\bf x}}) + \sum_{i = 1}^{l} \bar{u}_i \nabla g_i({\bar{\bf x}}) + \sum_{i = 1}^{m} {\bar v}_i \nabla h_i({\bar{\bf x}}) = {\bf 0} \\ & {\bar u}_i {\bar s}_i= \rho, \ \ \ g_i({\bar{\bf x}}) + {\bar s}_i = 0 , \ \ \ {\bar u}_i > 0, \ \ \ {\bar s}_i > 0 \ \ \ (i = 1, \cdots, l) \\ & h_i({\bar{\bf x}}) = 0 \ \ \ (i = 1, \cdots, m) \\ \end{align} \tag{5}$

が満たされる。そのため更新方向は式（5）を近似的に満たすように決める。つまり、

$\displaystyle \begin{align} & \nabla f({\bf x}^{(k)} + {\bf \Delta x}) + \sum_{i = 1}^{l} {u}_i^{(k)} \nabla g_i({\bf x}^{(k)} + {\bf \Delta x}) + \sum_{i = 1}^{m} v_i^{(k)} \nabla h_i({\bf x}^{(k)} + {\bf \Delta x}) = {\bf 0} \\ & (u_i^{(k)} + {\Delta u}_i) (s_i^{(k)} + {\Delta s}_i) = \rho_k, \ \ \ g_i({\bf x}^{(k)} + {\bf \Delta x}) + (s_i^{(k)} + {\Delta s}_i) = 0 \\ & (u_i^{(k)} + {\Delta u}_i) > 0, \ \ \ (s_i^{(k)} + {\Delta s}_i) > 0 \ \ \ (i = 1, \cdots, l) \\ & h_i({\bf x}^{(k)} + {\bf \Delta x}) = 0 \ \ \ (i = 1, \cdots, m) \\ \end{align} \tag{6}$

を満たすように( ${\bf \Delta x}, {\bf \Delta s}, {\bf \Delta u}, {\bf \Delta v})$ を決めたいのだが、式（6）の方程式を解くのは困難なので1次近似して解く。1次近似なので、

$\displaystyle \begin{align} \nabla f({\bf x}^{(k)} + {\bf \Delta x}) & \simeq \nabla f({\bf x}^{(k)} ) +\nabla^2 f({\bf x}^{(k)}) {\bf \Delta x} \\ \nabla g_i({\bf x}^{(k)} + {\bf \Delta x}) & \simeq \nabla g_i({\bf x}^{(k)} ) +\nabla^2 g_i({\bf x}^{(k)}) {\bf \Delta x} \\ \nabla h_i({\bf x}^{(k)} + {\bf \Delta x}) & \simeq \nabla h_i({\bf x}^{(k)} ) +\nabla^2 h_i({\bf x}^{(k)}) {\bf \Delta x} \\ g_i({\bf x}^{(k)} + {\bf \Delta x}) & \simeq g_i({\bf x}^{(k)} ) +(\nabla g_i({\bf x}^{(k)}))^T {\bf \Delta x} \\ h_i({\bf x}^{(k)} + {\bf \Delta x}) & \simeq h_i({\bf x}^{(k)} ) +(\nabla h_i({\bf x}^{(k)}))^T {\bf \Delta x} \end{align} \tag{7}$

を式（6）の等式に代入し、2次項 $\Delta u_i \Delta s_i$ を無視すると、以下の1次連立方程式が得られる。

$\displaystyle \left( \begin{array}{rrrr} \nabla_{{\bf x}{\bf x}}^2 L & {\bf O} & {\bf J}_g^T & {\bf J}_h^T \\ {\bf O} & {\bf D}_u & {\bf D}_s & {\bf O} \\ {\bf J}_g & {\bf I} & {\bf O} & {\bf O} \\ {\bf J}_h & {\bf O} & {\bf O} & {\bf O} \end{array} \right) \left( \begin{array}{r} {\bf \Delta x} \\ {\bf \Delta s} \\ {\bf \Delta u} \\ {\bf \Delta v} \\ \end{array} \right) = \left( \begin{array}{r} -{\nabla}_{\bf x} L \\ {\bf \rho_k {\bf 1} - {\bf s} \odot {\bf u}} \\ -({\bf g} + {\bf s}) \\ -{\bf h} \\ \end{array} \right) \tag{8}$

ただし、

$\displaystyle \begin{align} &\nabla_{{\bf x}{\bf x}}^2 L = \nabla^2 f({\bf x}^{(k)}) + \sum_{i = 1}^{l} {u}_i^{(k)} \nabla^2 g_i({\bf x}^{(k)}) + \sum_{i = 1}^{m} {v}_i^{(k)} \nabla^2 h_i({\bf x}^{(k)}), \\ &{\bf J}_g =\left( \nabla g_1({\bf x}^{(k)}), \cdots, \nabla g_l({\bf x}^{(k)}) \right)^T, \\ &{\bf J}_h =\left( \nabla h_1({\bf x}^{(k)}), \cdots, \nabla h_m({\bf x}^{(k)}) \right)^T, \\ &{\bf D}_u = {\rm diag} ({\bf u^{(k)}}), \\ &{\bf D}_s = {\rm diag} ({\bf s^{(k)}}), \\ &{\nabla}_{\bf x} L =\nabla f({\bf x}^{(k)}) + \sum_{i = 1}^{l} {u}_i^{(k)} \nabla g_i({\bf x}^{(k)}) + \sum_{i = 1}^{m} {v}_i^{(k)} \nabla h_i({\bf x}^{(k)}) ,\\ &{\bf g} =\left( g_1({\bf x}^{(k)}), \cdots, g_l({\bf x}^{(k)}) \right)^T, \\ &{\bf h} =\left( h_1({\bf x}^{(k)}), \cdots, h_m({\bf x}^{(k)}) \right)^T ,\\ &{\bf 1} = \left( 1, \cdots, 1 \right)^T \end{align} \tag{9}$

であり、 $\odot$ はアダマール積を表す。この式（8）を解くことで更新方向が得られる。

ステップ幅の決め方

次にステップ幅 $\alpha$ だが、まず ${\bf s}^{(k)} + \alpha {\bf \Delta s} > 0, {\bf u}^{(k)} + \alpha {\bf \Delta u} > 0$ となる必要があるため、これを満たす $0 \lt \alpha \leq 1$ を見つける*3。更にその条件を満たす中でもなるべく適切な $\alpha$ を選びたいので、以下の目的関数と制約のペナルティを組み合わせたメリット関数を定義し、

$\displaystyle \phi_{(\eta, \rho)}({\bf x}, {\bf s}) = f({\bf x}) - \rho \sum_{i = 1}^{l} \log s_i+ \eta \sum_{i = 1}^{l} |g_i({\bf x}) + s_i| + \eta \sum_{i = 1}^{m} |h_i({\bf x})| \tag{10}$

$\phi_{(\eta, \rho_k)}({\bf x}^{(k)} + \alpha {\bf \Delta x} , {\bf s}^{(k)} + \alpha {\bf \Delta s} )$ を最小にする $\alpha$ を直線探索により見つける*4。 $\eta > 0$ はペナルティの重みを表すパラメータで、今回の実装では固定値をあらかじめ与えている。

$\rho$ の更新

最後に $\rho$ の更新だが、これは ${\bf x}^{(k + 1)}, {\bf s}^{(k + 1)}, {\bf u}^{(k + 1)}, {\bf v}^{(k + 1)}$ が式（6）を満たすのなら $s_i^{(k + 1)} u_i^{(k + 1)} = \rho_k$ なので*5、 $\rho_k = \frac{({{\bf s}^{(k + 1)}})^{T} {\bf u}^{(k + 1)}}{l}$ となっているのが理想的である。このことから、適切なパラメータ $0 \lt \delta \lt 1$ をあらかじめ与えて、

$\displaystyle \rho_{k + 1} = \delta \frac{({{\bf s}^{(k + 1)}})^{T} {\bf u}^{(k + 1)}}{l} \tag{11}$

として更新する。

アルゴリズムまとめ

主双対内点法のアルゴリズムは以下の通り。

初期点 ${\bf x}^{(0)}, {\bf s}^{(0)}, {\bf u}^{(0)}, {\bf v}^{(0)}$ とバリア関数のパラメータの初期値 $\rho_0$ を定める。
$\rho_k$ が十分に小さければ終了。
式（8）を解いて探索方向 $( {\bf \Delta x}, {\bf \Delta s}, {\bf \Delta u}, {\bf \Delta v})$ を決める。
ステップ幅 $\alpha$ を決める。
$({\bf x}^{(k + 1)}, {\bf s}^{(k + 1)}, {\bf u}^{(k + 1)}, {\bf v}^{(k + 1)}) = ({\bf x}^{(k)} + \alpha {\bf \Delta x}, {\bf s}^{(k)} + \alpha {\bf \Delta s}, {\bf u}^{(k)} + \alpha {\bf \Delta u}, {\bf v}^{(k)} + \alpha {\bf \Delta v})$ として更新する。
式（11）で $\rho_{k+1}$ を求める。
2に戻る。

その他補足

実装する上で困ったこと、まだよく分かってないこと、その他細かい話など。

初期点について

内点法は基本的に実行可能領域の内部を通って解を見つけるアルゴリズムなので初期点も内部に置かないといけないと思ったのだが、実は今回述べたアルゴリズムだと初期点を実行可能領域外の点にしても試した範囲では実行可能領域の停留点に収束した。

今回の方法が実行可能領域外を初期点にしても収束するという理論的保証があるかはよく分かっていないが、雰囲気的には $\rho$ が十分に大きければ式（4）の最適解は実行可能領域内の境界から遠く離れた点であり、式（8）を満たす方向に点を動かすということはその式（4）の最適解に近づけることを意味するので、たとえ初期点が実行可能領域外だったとしても内部に入る方向に動き最終的には内部の停留点に収束するという感じだろう*6。ということで実装では、適当な値を初期値としている。

非凸な最適化問題の場合

ラグランジュ関数のヘッセ行列 $\nabla_{{\bf x}{\bf x}}^2 L$ （式（9）参照）が正定値でない場合、式（8）を満たす ${\bf \Delta x}$ が下降方向とは限らない。一方、ステップ幅を探索するときは、なるべくメリット関数が小さくなる $\alpha$ を探すので、 ${\bf \Delta x}$ が上昇方向であれば限りなく $0$ に近い $\alpha$ が得られなかなか収束しないということが起こりうる。そのため、今回の実装では上述の通り $\alpha$ が小さくなりすぎないよう探索の範囲を絞っている。また、そのほかの対策方法として、ヘッセ行列の代わりに準ニュートン法同様にヘッセ行列を近似する正値対称行列を用いる場合もある模様。

実装

以下が実装した主双対内点法のpythonコード。

import numpy as np
import copy

# 以下の最適化問題を定義するクラス
# min f(x) : objective
# s.t. g_i(x) <= 0 (i = 1..l) : inequality_constraints_list
#      h_i(x) = 0  (i = 1..m) : equality_constraints_list
# objective, inequality_constraints_list[i], equality_constraints_list[i]は、
# それぞれ以下のメンバを持ってる必要あり
#      func(x)   : return f(x)
#      grad(x)   : return ∇f(x)
#      hessian(x): return ∇^2 f(x)
#      dim       : 変数の次元
class problem():
    def __init__(self, objective, 
                inequality_constraints_list,
                equality_constraints_list):

        for constraint in inequality_constraints_list:
            if objective.dim != constraint.dim:
                raise Exception('dimension mismatch')
        
        for constraint in equality_constraints_list:
            if objective.dim != constraint.dim:
                raise Exception('dimension mismatch')
        
        self.dim = objective.dim
        self.objective = objective
        self.inequality_constraints_list = inequality_constraints_list
        self.equality_constraints_list = equality_constraints_list

def merit_function(problem, rho, eta, x, s):
    res = problem.objective.func(x)
    res -= rho * np.sum(np.log(s))
    g_sum = 0
    for i, constraint in enumerate(problem.inequality_constraints_list):
        g_sum += np.abs(constraint.func(x) + s[i])
    res += eta * g_sum
    h_sum = 0
    for i, constraint in enumerate(problem.equality_constraints_list):
        h_sum += np.abs(constraint.func(x))
    res += eta * h_sum 
    return res 

def backtracking_line_search(problem, x, s, u, dx, ds, du, rho, eta, beta):
    phi = lambda x, s:merit_function(problem, rho, eta, x, s)

    alpha1 = 1.0
    alpha2 = 1.0

    # αの初期値はs + α ds > 0 and u + α du > 0を満たす最大値
    # ただし、α <= 1
    if(np.min(ds) < 0):
        alpha1 = np.min(-s[ds < 0]/ds[ds < 0])
    if(np.min(du) < 0):
        alpha2 = np.min(-u[du < 0]/du[du < 0])
    alpha = min(1.0, alpha1 * beta, alpha2 * beta)

    min_merit = np.inf
    res_alpha = alpha
    init_alpha = alpha

    # 探索範囲はαの初期値の1/10
    while(alpha > init_alpha * 0.1):
        new_x = x + alpha * dx
        new_s = s + alpha * ds

        m = phi(new_x, new_s)
        if (m < min_merit):
            min_merit = m
            res_alpha = alpha
        alpha *= beta
    return res_alpha

# problem : 最適化問題
# eps     : 収束判定基準
# eta     : メリット関数パラメータ
# beta    : 直線探索の減衰パラメータ
# t       : ρ更新時のパラメータ
def interior_point(problem, eps = 1e-9, eta =0.1, beta = 0.9, t = 0.5):
    
    l = len(problem.inequality_constraints_list) #不等式制約数
    m = len(problem.equality_constraints_list) #等式制約数
    n = problem.dim #変数の次元

    rho = 1
    # 初期値xはとりあえず制約を満たさなくてもよい
    x = np.ones(n)
    # s の初期値は s > 0
    s = np.ones(l)
    # s_i u_i = ρを満たすようにuを設定
    u = rho / s

    v = np.zeros(m)
    
    while(1):
        if rho < eps:
            break
        H = problem.objective.hessian(x)
        d = problem.objective.grad(x)
        Jg = np.zeros((l, n))
        g = np.zeros(l)
        # 不等式制約の関数がらみの計算
        for i, constraint in enumerate(problem.inequality_constraints_list):
            H += constraint.hessian(x) * u[i]
            Jg[i, :] = constraint.grad(x)
            d += Jg[i, :] * u[i]
            g[i] = constraint.func(x)
        
        Jh = np.zeros((m, n))
        h = np.zeros(m)
        # 等式制約の関数がらみの計算
        for i, constraint in enumerate(problem.equality_constraints_list):
            H += constraint.hessian(x) * v[i]
            Jh[i, :] = constraint.grad(x)
            d += Jh[i, :] * v[i]
            h[i] = constraint.func(x)

        Du = np.diag(u)
        Ds = np.diag(s)

        # make array P       
        # |H   O   Jg.T  Jh.T|
        # |O   Du  Ds    O   |
        # |Jg  I   O     O   |
        # |Jh  O   O     O   | 
        P = np.zeros((n + 2 * l + m, n + 2 * l + m))
        P[:n, :n] = H
        P[:n, n + l: n + 2 * l] = Jg.T
        P[:n, n + 2 * l:] = Jh.T
        P[n:n + l, n:n + l] = Du
        P[n:n + l, n + l:n + 2 * l] = Ds
        P[n + l:n + 2 * l, :n] = Jg
        P[n + l:n + 2 * l, n:n + l] = np.eye(l)
        P[n + 2 * l:, :n] = Jh

        # make vector r
        r = np.concatenate([-d, rho - s * u, -(g + s), -h])

        delta = np.linalg.solve(P, r)

        dx = delta[:n]
        ds = delta[n:n + l]
        du = delta[n + l:n + 2 * l]
        dv = delta[n + 2 * l:]

        alpha = backtracking_line_search(problem, x, s, u,
                                            dx, ds, du, rho, eta, beta)

        # x,s,u,vの更新
        x += alpha * dx
        s += alpha * ds
        u += alpha * du
        v += alpha * dv

        # ρの更新
        rho =  t * u @ s / l

    return x, s, u, v

式（1）の最適化問題はproblemというクラスに記述し、そのインスタンスを関数interior_pointに渡すと停留点を探して出力する。problemのobjectiveに目的関数を表すクラスのインスタンスをセットし、inequality_constraints_listには不等式制約の関数（式（1）の $g_i({\bf x})$ ）を表すクラスのインスタンスのリスト、inequality_constraints_listには等式制約の関数（式（1）の $h_i({\bf x})$ ）を表すクラスのインスタンスのリストをそれぞれセットする。各クラスは $f({\bf x}), g_i({\bf x}), h_i({\bf x})$ の値を返すメソッドfunc、 $\nabla f({\bf x}), \nabla g_i({\bf x}),\nabla h_i({\bf x})$ を返すメソッドgrad、 $\nabla^{2} f({\bf x}), \nabla^{2} g_i({\bf x}),\nabla^{2} h_i({\bf x})$ を返すメソッドhessianを実装してる必要がある。

なお、式（9）の連立方程式はnp.linalg.solveで解いているが、零要素が多いためもう少し効率的に解く方法があるかもしれない。また、この実装だと不等式制約が存在しない問題は解けない（等式制約が存在しない問題は解ける）。

ついでにSVMの最適化問題を内点法で解く

SVMは分類問題を解く教師あり学習の一種であり、その「学習」は最適化問題を解くことに対応する。通常はSMOなどの専用の効率が良い最適化アルゴリズムなどで解くが、今回はこれを主双対内点法で解いてみる。

SVMで解く最適化問題

SVMの学習で解くことになる最適化問題を簡単に説明する。SVMの学習では教師データ集合 $\left\{({\bf x}_i, y_i) \right\}_{i=1}^N, ({\bf x}_i \in \mathbb{R}^n, y_i \in \left\{-1, 1 \right\})$ が得られてるとして、それをうまいこと分離できる関数 $f$ とバイアス項 $b$ を求める。具体的にはヒンジ損失関数に正則化項を加えた以下の最小化問題を解く。

$\displaystyle \min_{f, b} \left\{ \frac{1}{2} \parallel f \parallel_{\mathcal H}^2+ C \sum_{i = 1}^{N} \max(1 - y_i (f({\bf x}_i) + b), 0) \right\} \tag{12}$

この問題は関数の $f$ に関する最適化問題であり、加えて第2項（ヒンジ損失関数）は微分できないので、このまま解くのは難しい。そこでリプレゼンター定理などを使ってゴニョゴニョしたりヒンジ損失関数の制約条件への変換を行ったりすると、以下の凸2次計画問題に変形できる。

$\displaystyle \begin{align} \min_{{\bf a}, b, {\bf \xi}} \ & \frac{1}{2}{\bf a}^T {\bf K} {\bf a} + C \sum_{i=1}^{N} \xi_i \\ {\rm s.t.} \ & - y_i \left( \sum_{j = 1}^N a_j k({\bf x}_i,{\bf x}_j)+ b \right) + 1 - \xi_i \leq 0, \ \xi_i \geq 0 \ (i = 1, \cdots, N) \end{align} \tag{13}$

ここで $k({\bf x}, {\bf x}')$ はカーネル関数と呼ばれ適切な正定値対称関数を選ぶ（例えば、 $k({\bf x}, {\bf x}') = \exp(-\gamma \parallel {\bf x} - {\bf x}' \parallel^{2})$ など）。また ${\bf K} \in{\mathbb R}^{N \times N}$ は $i$ 行 $j$ 列目に $k(x_i, x_j)$ を持つグラム行列である。これらの理屈などについては、以下の資料や参考文献として挙げた『カーネル多変量解析』という書籍を参考にするとよい*7。

https://www.ism.ac.jp/~fukumizu/ISM_lecture_2004/Lecture2004_kernel_method.pdf

この最適化問題式（13）の解 ${\bf a}^{*}$ と $b^{*}$ を用いて、新たに与えられた入力 ${\hat{\bf x}}$ は

$\displaystyle {\rm sign} \left( \sum_{i = 1}^N a_i^{*} k({\bf x}_i, {\hat {\bf x}})+ b^{*} \right) \tag{14}$

としてクラスを分類できる。

SVMの実装

以下が実装したSVMとそれを用いて検証用データの学習と分類を行うpythonコード。上述の主双対内点法の関数を利用している。

import numpy as np
import interior_point_method as ipm
import matplotlib.pyplot as plt

# 目的関数を表すクラス
class objective_function:
    def __init__(self, K, C, var_dim, constraint_num):
        self.K = K
        self.C = C
        self.a_dim = var_dim
        self.z_dim = constraint_num # ξはzで表現
        self.dim = var_dim + constraint_num + 1 # aの次元 + ξの次元 + bの次元

    # a_z_bはaとξとbを1つに並べたベクトル（array）
    def func(self, a_z_b):
        return 0.5 * a_z_b[:self.a_dim] @ self.K @ a_z_b[:self.a_dim] + \
                self.C * np.sum(a_z_b[self.a_dim:self.a_dim + self.z_dim])
           
    def grad(self, a_z_b):
        return np.concatenate([self.K @ a_z_b[:self.a_dim],
                               self.C * np.ones(self.z_dim),
                               [0]])

    def hessian(self, a_z_b):
        H = np.zeros((self.dim, self.dim))
        H[:self.a_dim, :self.a_dim] = self.K
        return H

# 教師データによる不等式制約クラスを生成するクラス
class train_data_constraint_creator:
    def __init__(self, K, y):
        self.K = K
        self.y = y

    # 教師データによる不等式制約を表すクラス
    class train_data_constraint:
        def __init__(self, K, y, idx):
            self.idx = idx
            self.K = K
            self.y = y
            self.a_dim = K.shape[0]
            self.z_dim = y.shape[0]
            self.dim = self.a_dim + self.z_dim + 1
        
        def func(self, a_z_b):
            return -a_z_b[self.a_dim + self.idx] + 1 - \
                    self.y[self.idx] * \
                    (a_z_b[:self.a_dim] @ self.K[self.idx, :] + a_z_b[-1])
             
        def grad(self, a_z_b):
            return np.concatenate([-self.K[self.idx, :] * self.y[self.idx],
                                    np.ones(self.z_dim) * -1,
                                    [-self.y[self.idx]]])

        def hessian(self, a_z_b):
            return np.zeros((self.dim, self.dim))

    def create(self, idx):
        return self.train_data_constraint(self.K, self.y, idx)

# 非負制約クラスを生成するクラス
class nonnegative_constraint_creator:
    def __init__(self, K, y):
        self.constraint_num = y.shape[0]
        self.var_dim = K.shape[0]
        self.constraint_num = y.shape[0]
        self.dim = self.var_dim + self.constraint_num + 1

    # 非負制約を表すクラス
    class function_nonnegative_constraint:
        def __init__(self, var_dim, constraint_num, idx):
            self.idx = idx
            self.a_dim = var_dim
            self.z_dim = constraint_num
            self.dim = self.a_dim + self.z_dim + 1
        
        def func(self, a_z_b):
            return -a_z_b[self.a_dim + self.idx]
        
        def grad(self, a_z_b):
            res = np.zeros(self.dim)
            res[self.a_dim + self.idx] = -1
            return res

        def hessian(self, a_z_b):
            return np.zeros((self.dim, self.dim))
             
    def create(self, idx):
        return self.function_nonnegative_constraint(self.var_dim,
                                                self.constraint_num, idx)

class svm:
    def __init__(self, C, kernel_func):
        self.C = C
        self.K = None
        self.X = None
        self.kernel_func = kernel_func

    def fit(self, X, y):
        self.X = X
        self.K = np.zeros((X.shape[0], X.shape[0]))

        # グラム行列の計算
        for i in range(X.shape[0]):
            for j in range(i, X.shape[0]):
                self.K[i, j] = self.kernel_func(X[i, :], X[j, :])
                self.K[j, i] = self.K[i, j]

        # 各種不等式制約インスタンス生成用
        tdcc = train_data_constraint_creator(self.K, y)
        ncc = nonnegative_constraint_creator(self.K, y)
        # 最適化問題のインスタンスを生成
        Pro = ipm.problem(objective_function(self.K, self.C,
                                            self.K.shape[1], y.shape[0]), 
                        [tdcc.create(i) for i in range(y.shape[0])] + 
                        [ncc.create(i) for i in range(y.shape[0])],
                        [])
        
        a_z_b_opt, s_opt, u_opt, v_opt = ipm.interior_point(Pro, eps = 1e-15)
        self.a = a_z_b_opt[:self.K.shape[0]]
        self.b = a_z_b_opt[-1]
        self.u = u_opt
        self.support_vectors = np.where(self.u[:X.shape[0]] > 1e-8)

    def predict(self, X_):
        if X_ is None:
             raise Exception('not learned')
            
        K_ = np.zeros((self.K.shape[0], X_.shape[0]))
        for i in range(self.K.shape[0]):
            for j in range(X_.shape[0]):
                K_[i, j] = self.kernel_func(self.X[i], X_[j])
        return np.sign(self.a @ K_ + self.b)


# テストデータ生成用
def make_data(xmin, xmax, n):
    np.random.seed(1111)
    cnt = 0
    t = 0
    N = n * 9
    X_train = np.zeros((N, 2))
    y_train = np.zeros(N)
    center1 =  np.linspace(-(-xmin[0] + xmax[0])/3.5, 
                            (-xmin[0] + xmax[0])/3.5, 3)
    center2 =  np.linspace(-(-xmin[1] + xmax[1])/3.5, 
                            (-xmin[1] + xmax[1])/3.5, 3)
    for i in range(3):
        for j in range(3):
            y_train[t:t + n] = 1 if cnt % 2 == 0 else -1
            X_train[t:t + n, :] = 1.5 * np.random.randn(n, 2) + \
                                np.array([center1[i], center2[j]])
            t += n
            cnt += 1

    #テスト入力点生成
    T1 = np.linspace(xmin[0], xmax[0], 60)
    T2 = np.linspace(xmin[1], xmax[1], 60)
    X_test = np.zeros((T1.shape[0] * T2.shape[0], 2))
    i = 0
    for t1 in T1:
        for t2 in T2:
            X_test[i, :] = np.array([t1, t2])
            i += 1
    return X_train, y_train, X_test

def main():
    #教師データとテストデータの生成
    xmin = np.array([-9, -9])
    xmax = np.array([9, 9])
    n = 20
    X_train, y_train, X_test = make_data(xmin, xmax, n)

    # カーネル関数（ガウスカーネル）
    gamma = 0.5
    kernel_func = lambda x_i, x_j : np.exp(-gamma * (x_i - x_j) @ (x_i - x_j))

    sv = svm(0.1, kernel_func)
    # SVM学習
    sv.fit(X_train, y_train)
    # テストデータの予測
    y_pred = sv.predict(X_test)
    y_pred[y_pred == -1] = 0

    # 以下結果のプロット
    size = 100
    sc = plt.scatter(X_test[:, 0], X_test[:, 1],
                        size, y_pred, marker = ',', cmap='bwr')

    plt.scatter(X_train[y_train == 1, 0], X_train[y_train == 1, 1] ,
                edgecolors="black", marker="o" , color='r',
                label= "train data (y = 1)")
    plt.scatter(X_train[y_train == -1, 0], X_train[y_train == -1, 1] ,
                edgecolors="black", marker="o" , color='b',
                label= "train data (y = -1)")
    plt.axis([xmin[0], xmax[0], xmin[1], xmax[1]])
    plt.legend(loc='upper left')
    plt.show()

if __name__ == '__main__':
    main()

先述のproblemクラスとinterior_point関数を使ってSVMを実装している。基本的にはこれまで述べてきたことをそのまま書いてるだけなので特筆すべきことはないと思う。これを実行すると以下の結果が得られる。

f:id:YamagenSakam:20210124133035p:plain

ガウスカーネルを使っているので非線形な分離できている。ちなみに、グラム行列 ${\bf K}$ が正則でない場合（例えば、 $\gamma$ を小さくしすぎた場合や全く同じ学習サンプルが2つ以上ある場合などにあり得る）、主双対内点法のコードにおけるPも正則でなくなり計算が破綻する。このあたりの弱点はSMO法などで克服される。

ここに挙げた主双対内点法のコードやSVMのコードに加えて、線形のSVMやその他の最適化問題で確認したコードもあるので、それらは全て以下のgithub リポジトリに上げている。

github.com

最後に

今回は制約付き非線形計画問題を解く主双対内点法を実装した。制約付き非線形計画問題を解く手法としては他にも逐次2次凸計画法や乗数法（拡張ラグランジュ関数法）などがある。ちなみに乗数法については、以前その改良の1つである交互方向乗数法（ADMM）を実装した記事を書いている。

yamagensakam.hatenablog.com

ということで、非線形計画問題については一旦ここまでにして次から整数計画問題に移る予定。

参考文献

いつも通り今回の内容もこちらの書籍を主に参考にして書いた。
しっかり学ぶ数理最適化モデルからアルゴリズムまで (KS情報科学専門書)
- 作者:梅谷俊治
- 発売日: 2020/10/26
- メディア: 単行本
SVMやカーネル法、リプレゼンター定理などについてはこちらの書籍がわかりやすい。
カーネル多変量解析―非線形データ解析の新しい展開 (シリーズ確率と情報の科学)
- 作者:赤穂昭太郎
- 発売日: 2008/11/27
- メディア: 単行本

*1:主双対内点法は数ある内点法の一種だが、単に内点法と書いて主双対内点法のことを指してる文献もある。

*2:制約想定以外にも主双対内点法について色々書かれているので、今書いてる自分の記事は不要じゃないかとも思ったが...まぁいいか。

*3: ${\bf s}^{(k)} > 0, {\bf u}^{(k)} > 0$ なので、条件を満たす $\alpha$ は必ず存在する。

*4:実装では条件を満たす中でもなるべく大きい $\alpha$ を初期値としてバックトラッキングを行っている。ただし、探索範囲の下限を小さくすると、極端に小さい $\alpha$ が得られることがあるので探索範囲の下限は初期値の $1/10$ 程度までにとどめている。

*5:1次近似して式（6）を解いているので一般には満たされない。

*6:実行可能とは限らない初期点を使うことのできる内点法を『インフィージブル内点法』というらしい。今回がその一手法に当たるかは不明。

*7:かなり説明を端折ったのでもう少し補足すると、SVMの学習をソフトマージンのマージン最大化問題として定式化→それは式（12）の正則化付きヒンジ損失最小化と等価→これにリプレゼンター定理を適用して有限次元の最適化問題に帰着→これをソフトマージンの問題に戻すと有限次元の2次計画問題である式（13）が出てくるという流れ。詳しくは貼った資料のP39~P45参照。

甲斐性なしのブログ

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

数理最適化をしっかり学ぶために主双対内点法とそれを使ったSVMを実装