# 知能機械情報学レポート
<div style="text-align: right;">03-180265 鎌田啓路</div>

## カーネルトリックを使ったSVM

今回はSVMを、ソフトマージンとガウスカーネルによるカーネルトリックを使って実装しました

理論的な説明としては、授業ノートにもあったとおり、ソフトマージンを入れた場合のSVMの目的関数は以下のようになります

\begin{align}
\min_{w, b, \xi} & \left(\frac{1}{2} ||w||^2 + C\sum_{n=1}^N \xi_n \right) \\
s.t.&~~ t_n y(x_n) \ge 1-\xi_n \\
&\xi_n \ge 0
\end{align}

$\xi_n$が小さいほど目的関数も小さくなるので、これはまとめて次のようにかけます
\begin{align}
\min_{w, b} \left( \frac{1}{2} ||w||^2 + C\sum_{n=1}^N \max(0, 1-t_n y(x_n) ) \right)
\end{align}

ここで、識別関数$y(x)$を、非線形関数$\phi$を用いて
\begin{equation}
y(x_n) = w^T\phi(x_n)+b
\end{equation}
のように変更します。
ここで、
\begin{equation}
w = \sum_{m=1}^N \theta_m \phi(x_m)
\end{equation}
の形に限定することで、以下の最適化問題に帰着できます。

\begin{align}
\min_\theta &\left( \frac{1}{2} \theta^T K \theta + C\sum_{n=1}^N \max(0, 1-t_n y_\theta(x_n) )\right) \\
&\phi(x_i)^T\phi(x_j) = K(x_i, x_j) \\
&y_\theta(x_n) = \sum_{m=1}^N \theta_m K(x_n, x_m) + \theta_0
\end{align}
(以下では便宜上$\theta_0$は無視することとします)

今回は、カーネル関数として、ガウスカーネル
\begin{equation}
K(x, x') = \exp\left(-\frac{||x-x'||^2}{2h^2}\right)
\end{equation}
を用いることとします。

この目的関数の最適化は劣勾配法で行うことができます

\begin{align}
\theta \gets \theta - \epsilon\left(K\theta +C\sum_{n=1}^N \partial_\theta \max(0, 1-t_n y_\theta(x_n)) \right)
\end{align}

この劣微分の第$m$成分を実際に計算すると以下のようになります。

\begin{align}
\partial_{\theta_m} \max(0, 1-t_n y_\theta(x_n)) = \left\{
\begin{array}{l}
-t_n K(x_n, x_m) &(1-t_n y_\theta(x_n) > 0)\\
0 & (1-t_n y_\theta(x_n) \le 0)\\
\end{array}
\right.
\end{align}

## 実験

今回は、このSVMを用いて、以下の式で生成された正例データ$x_p$と負例データ$x_m$を識別するタスクを行いました。
（二値分類のタスクで線形分離が難しいデータセットを見つけるのが難しかったので人工データを用いることにしました）

\begin{align}
x_p &= ~(a_i \cos a_i+ \epsilon_{p,i,1}, ~ a_i \sin a_i+ \epsilon_{p,i,2})\\
x_m &= ~((a_i+\pi) \cos a_i+ \epsilon_{m,i,1}, ~ (a_i + \pi) \sin a_i+ \epsilon_{m,i,2})\\
a_i &= \frac{10(i-1)}{n-1}\\
\epsilon &\sim N(0,1)
\end{align}

最適化はアルミホ条件を満たすまでバックトラッシュを行っています。

実装したコードは以下のようになります。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

def make_data(n):
    """
    データの生成
    input:
        n: データ数
    output:
        X[perm]: 学習データ (n, 2) (x,y軸方向の座標)
        T[perm]: 教師ラベル (n,) (1 or -1)
    """
    a = np.array([10*i/(n//2-1) for i in range(n//2)])
    x_p_1 = a*np.cos(a) + np.random.rand(n//2)
    x_p_2 = a*np.sin(a) + np.random.rand(n//2)
    x_p = np.vstack([x_p_1, x_p_2])
    t_p = np.ones_like(x_p_1)
    x_m_1 = (a+np.pi)*np.cos(a) + np.random.rand(n//2)
    x_m_2 = (a+np.pi)*np.sin(a) + np.random.rand(n//2)
    x_m = np.vstack([x_m_1, x_m_2])
    t_m = -np.ones_like(x_m_1)
    X = np.vstack([x_p.T, x_m.T])
    T = np.hstack([t_p.T, t_m.T])
    perm = np.random.permutation(n)
    return X[perm], T[perm]

def calc_K(X1, X2, h=1):
    """
    ガウスカーネルの計算
    input: 
        X1: (1, N, 2)
        X2: (M, 1, 2)
    output:
        K: (N, M) (X1, X2の要素の全ての組み合わせをガウス関数に通した行列)
    """
    K = np.exp(-np.sum(((X1 - X2))**2, axis=2)/(2*h*h))
    return K

def y_theta(X, theta, h=1):
    """
    識別関数y_\thetaの計算
    input: 
        X: 入力データ (n, 2)
        theta: カーネル関数の重み (n,)
    output: y_\thetaの値 (n,)
    """
    K = calc_K(X[None,:,:], X[:,None,:], h)
    return np.sum(theta * K, axis=1)

def y_theta_test(test_X, X, theta, h=1):
    """
    テスト時のy_\thetaの計算
    test_Xはメッシュを入力
    """
    K = calc_K(X[None,:,:], test_X[:,None,:], h)
    return np.sum(theta * K, axis=1)

def f(X, T, theta, C=0.1, h=1):
    """
    目的関数fの計算
    input: 
        X: 入力データ (n, 2)
        T: 教師ラベル (n,)
        theta: カーネル関数の重み (n,)
    output:
        目的関数の値(スカラー)
    """
    K = calc_K(X[None,:,:], X[:,None,:], h)

    return C * np.sum(np.fmax(1-y_theta(X, theta)*T, 0)) +\
            theta@K@theta / 2.0

def nabla_f(X, T, theta, C=0.1, h=1):
    """
    目的関数fの劣微分の計算
    input: 
        X: 入力データ (n, 2)
        T: 教師ラベル (n,)
        theta: カーネル関数の重み (n,)
    output:
        目的関数の各thetaにおける劣微分値 (n,)
    """
    K = calc_K(X[None,:,:], X[:,None,:], h)
    partial = 1-y_theta(X, theta)*T
    return C*np.sum(- T * (partial > 0)*K, axis=1) + K@theta

def run(N=100, alpha=0.5, beta=0.9, eps=0.5, plot=True):
    N = 100
    alpha = 0.5
    beta = 0.9
    eps = 0.5

    X, T = make_data(N)
    theta = np.random.rand(N)

    for _ in range(80):
        # アルミホ条件を満たすまでバックトラッシュ
        while f(X, T, theta - eps*nabla_f(X, T, theta)) \
                - f(X, T, theta) > \
                - alpha * eps * np.sum(nabla_f(X,T,theta)**2):
            eps = eps * beta
        theta = theta - eps*nabla_f(X, T, theta)
    if plot:
        plot_result(X, T, theta)
    return

def plot_result(X, T, theta):
    plt.scatter(X[:,0][T>0], X[:,1][T>0], c='red')
    plt.scatter(X[:,0][T<0], X[:,1][T<0], c='blue')
    plt.grid(True)

    pred_T = y_theta(X, theta)
    x1_mesh = np.arange(-15, 15, 0.1)
    mesh_X = np.array(np.meshgrid(x1_mesh, x1_mesh)).transpose(1,2,0)
    pred_T = y_theta_test(mesh_X.reshape(-1,2), X, theta)
    pred_T = pred_T.reshape((len(x1_mesh), len(x1_mesh)))
    pred_T = np.where(pred_T > 0, -1, 1)
    plt.contourf(x1_mesh, x1_mesh, pred_T, alpha=0.4, \
                 cmap=ListedColormap(('red', 'blue')))
    plt.xlim(-15, 15)
    plt.savefig('result_svm.png')
    
if __name__ == '__main__':
    run()

## 実験結果
結果を表示すると以下のようになり、線形分離が難しいデータにおいても決定境界をうまく引くことができていることがわかります。
今回は、人工データを用いて実験を行ったので、汎化性能の評価は行っていません。

以下のコードでデータ数Nに対する計算時間の変化を計測しました。

import time

def calc_time(range_N=[50, 100, 200, 400]):
    result = []
    for N in range_N:
        print(N)
        times = []
        for _ in range(10):
            start_time = time.time()
            run(N=N, plot=False)
            end_time = time.time()
            times.append(end_time - start_time)
            print(end_time - start_time)
        result.append([N, sum(times) / len(times)])
    return np.array(result)
    
def plot_calc_time():
    calc_times = calc_time()
    plt.figure()
    plt.scatter(calc_times[:, 0], calc_times[:, 1], label='measured time')
    x = np.linspace(0,400,400)
    res = np.polyfit(calc_times[:, 0], calc_times[:, 1], 2)
    pred = np.poly1d(res)(x)
    plt.plot(x, pred, label='approx curve')
    plt.xlabel('N')
    plt.ylabel('calculation time (sec)')
    plt.legend()
    plt.savefig('calctime.png')

結果は以下の図のようになり、計算量は$O(N^2)$とわかりました。各行列演算の計算量が$O(N^2)$だと思うので（numpyの中身の詳しい実装はわからない...）、理論値と一致すると思います。