統計的機械学習の数理100問（19/20） - 「大人の教養・知識・気付き」を伸ばすブログ

　いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として

統計的機械学習の数理100問 with R (with R) (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

統計的機械学習の数理100問 with Python (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

を用いることにする。

前回

power-of-awareness.com

前回
8.　サポートベクトルマシン

8.　サポートベクトルマシン

8.2　最適化の理論

　主問題、双対問題の最適値をそれぞれ

$\begin{aligned} f^{*}&=\displaystyle{\inf_{\boldsymbol{\beta}}f_0(\boldsymbol{\beta})},\\ g^{*}&=\displaystyle{\sup_{\boldsymbol{\beta}}f_0(\boldsymbol{\beta})} \end{aligned}$

とおくと、 $\boldsymbol{\alpha}\geq\boldsymbol{0},f_1(\boldsymbol{\beta})\leq0,\cdots,f_m(\boldsymbol{\beta})\leq0$ に対して

$\begin{aligned} \displaystyle{\sum_{i=1}^{m}\alpha_if_i(\boldsymbol{\beta})}\leq0 \end{aligned}$

であるから、

$\begin{aligned} f^{*}\geq g^{*}=\displaystyle{\sup_{\boldsymbol{\alpha}\geq\boldsymbol{0}}g(\boldsymbol{\alpha})} \end{aligned}$

が成立する。以下、この不等式の等式が成立していると仮定する*1。

　次に $f_0,f_1,\cdots,f_m:\mathbb{R}^m\rightarrow\mathbb{R}$ が凸関数で $\boldsymbol{\beta}=\boldsymbol{\beta}^{*}$ で微分可能だとする。以下の条件を $\mathrm{KKM}$ 条件という。

$\mathrm{KKM}$ 条件　 $f_1(\boldsymbol{\beta})\leq0,\cdots,f_m(\boldsymbol{\beta})\leq0$ の下で $\boldsymbol{\beta}=\boldsymbol{\beta}^{*}\in\mathbb{R}^p$ が $f_0(\boldsymbol{\beta})$ を最小にすることと、

$\begin{aligned} f_1(\boldsymbol{\beta}^{*}),\cdots,f_m(\boldsymbol{\beta}^{*})\leq0 \end{aligned}$

であって

$\begin{aligned} \alpha_1 f_1(\boldsymbol{\beta}^{*})=\cdots=\alpha_m f_m(\boldsymbol{\beta}^{*})=0,\\ \nabla f_0(\boldsymbol{\beta}^{*})+\displaystyle{\sum_{i=1}^{m}\alpha_i \nabla f_i(\boldsymbol{\beta}^{*})}=0 \end{aligned}$

を満足するような $\alpha_1,\cdots\alpha_m\geq0$ が存在することは同値である。

( $\Longleftarrow$ )　一般に $f:\mathbb{R}^p\rightarrow\mathbb{R}$ が凸であって、 $x=x_0\in\mathbb{R}$ において微分可能であるとき、各 $\boldsymbol{x}\in\mathbb{R}^p$ について

$\begin{aligned} f(\boldsymbol{x})\geq f(\boldsymbol{x_0})+{}^{t}\nabla f(\boldsymbol{x}_0)(f(\boldsymbol{x})-f(\boldsymbol{x})_0) \end{aligned}$

が成立する。このことから、 $\mathrm{KKT}$ 条件の下で $\boldsymbol{\beta}\in\mathbb{R}^p$ について

$\begin{aligned} f_0(\boldsymbol{\beta}^{*})&\leq f_0(\boldsymbol{\beta})-\nabla{}^{t}f_0(\boldsymbol{\beta}^{*})(\boldsymbol{\beta}-\boldsymbol{\beta}^{*})\\ &\leq f_0(\boldsymbol{\beta})+\displaystyle{\sum_{i=1}^{m}\alpha_i{}^{t}\nabla f_i(\boldsymbol{\beta}^{*}) }(\boldsymbol{\beta}-\boldsymbol{\beta}^{*})\\ &\leq f_0(\boldsymbol{\beta})+\displaystyle{\sum_{i=1}^{m}\alpha_i(f_i(\boldsymbol{\beta})-f_i(\boldsymbol{\beta}^{*})) }\\ &=f_0(\boldsymbol{\beta})+\displaystyle{\sum_{i=1}^{m}\alpha_i f_i(\boldsymbol{\beta}) }\\ &\leq f_0(\boldsymbol{\beta}) \end{aligned}$

が成立し、 $\boldsymbol{\beta}^{*}$ が最適解であることが分かる。

( $\Longrightarrow$ )　 $\boldsymbol{\beta}^{*}$ が主問題の最適解であれば、 $f_1(\boldsymbol{\beta})\leq0,\cdots,f_m(\boldsymbol{\beta})\leq0$ の解 $\boldsymbol{\beta}$ が存在することを意味するため、

$\begin{aligned} f_1(\boldsymbol{\beta}^{*}),\cdots,f_m(\boldsymbol{\beta}^{*})\leq0 \end{aligned}$

でなければならない。仮定 $f^{*}=g^{*}$ から

$\begin{aligned} {}^{\exists}\boldsymbol{\alpha}^{*}\geq\boldsymbol{0}\left(f_0(\boldsymbol{\beta}^{*})=g(\boldsymbol{\alpha}^{*})=\displaystyle{\inf_{\boldsymbol{\beta}}\left\{f_0(\boldsymbol{\beta})+\sum_{i=1}^{m}\alpha_i^{*}f_i(\boldsymbol{\beta})\right\}}\leq f_0(\boldsymbol{\beta}^{*})+\displaystyle{\sum_{i=1}^{m}\alpha_i^{*}f_i(\boldsymbol{\beta}^{*})}\leq f_0(\boldsymbol{\beta}^{*})\right) \end{aligned}$

が成り立つ。すなわち2つの不等式はすべて等式となり、

$\begin{aligned} \alpha_1 f_1(\boldsymbol{\beta}^{*})=\cdots=\alpha_m f_m(\boldsymbol{\beta}^{*})=0 \end{aligned}$

が成立する。また $\boldsymbol{\beta}^{*}$ は $f_0(\boldsymbol{\beta})+\displaystyle{\sum_{i=1}^{m} \alpha_i^{*}f_i(\boldsymbol{\beta})}$ の最小値を取ることで得られたものだから、

$\begin{aligned} \nabla f_0(\boldsymbol{\beta}^{*})+\displaystyle{\sum_{i=1}^{m}\alpha_i \nabla f_i(\boldsymbol{\beta}^{*})}=0 \end{aligned}$

も成り立つ。　 $\blacksquare$ )

8.3　サポートベクトルマシンの解

　 $\mathrm{KKT}$ 条件から以下の7つの式が得られる。

$\begin{aligned} \begin{cases} y_i(\beta_0+x_i\beta)-(1-\varepsilon_i)\geq0\\ \alpha_i\{y_i(\beta_0+x_i\beta)-(1-\varepsilon_i)\}=0\\ \varepsilon_i\geq0\\ \boldsymbol{\beta}=\displaystyle{\sum_{i=1}^{n}\alpha_i y_i\boldsymbol{x}_i }\in\mathbb{R}^p\\ \displaystyle{\sum_{i=1}^{n}\alpha_i y_i}=0\\ C-\alpha_i-\mu_i=0 \end{cases} \end{aligned}$

　 $L_P$ の双対問題は

$\begin{aligned} \displaystyle{\sum_{i=1}^{n}\alpha_i}+\displaystyle{\frac{1}{2}\|\boldsymbol{\beta}\|^2}-\displaystyle{\sum_{i=1}^{n}x_i{}^{t}\boldsymbol{\beta}\alpha_i\boldsymbol{y}_i} \end{aligned}$

と書ける。更に $\boldsymbol{\beta}$ で微分して $0$ とおいた式を用いると

$\begin{aligned} \displaystyle{\sum_{i=1}^{n}\alpha_i}+\displaystyle{\frac{1}{2}\|\boldsymbol{\beta}\|^2}-\displaystyle{\frac{1}{2}}{}^{t}\left(\displaystyle{\sum_{i=1}^{n}\alpha_iy_i{}^{t}\boldsymbol{x}_i }\right)\left(\displaystyle{\sum_{j=1}^{n}\alpha_j y_j{}^{t}\boldsymbol{x}_j }\right)-\displaystyle{\sum_{i=1}^{n}\boldsymbol{x}_i}\left(\displaystyle{\sum_{j=1}^{n}\alpha_jy_j{}^{t}\boldsymbol{x}_j}\right)\alpha_iy_i \end{aligned}$

で得られる。以上から $\mathrm{Lagrange}$ 係数だった $\alpha_i,\mu_i\geq0,i=1,\cdots,n$ を入力とする関数

$\begin{aligned} L_D=\displaystyle{\sum_{i=1}^{n}\alpha_i}-\displaystyle{\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_jy_iy_j\boldsymbol{x}_i{}^{t}\boldsymbol{x}_j } \end{aligned}$

を構成できる。 $\boldsymbol{\alpha}$ は

$\begin{aligned} \displaystyle{\sum_{i=1}^{n}\alpha_i y_i}=0,\\ 0\leq\alpha_i\leq C \end{aligned}$

を満たすような範囲で動く。 $\mu_i$ は $L_D$ の中に亜h含まれていないものの、 $\mu_i=C-\alpha_i\geq0$ が $\alpha_i\geq0$ と共に上記の $0\leq\alpha_i\leq C$ として反映されている。

　ここで双対問題を解くに当たり、 $0\leq\alpha_i\leq C$ が以下の3つの場合に分類できることに注意する。

$\alpha_i$ に対する条件

$\begin{aligned} \begin{cases} \alpha_i=0&\Longleftarrow y_i(\beta_0+{}^t\boldsymbol{x}_i\boldsymbol{\beta})\gt1\\ 0\lt\alpha_i\lt C&\Longrightarrow y_i(\beta_0+{}^t\boldsymbol{x}_i\boldsymbol{\beta})=1\\ \alpha_i =C&\Longleftarrow y_i(\beta_0+{}^t\boldsymbol{x}_i\boldsymbol{\beta})\lt1 \end{cases} \end{aligned}$

( $\because$ 　 $\alpha_i=0$ のとき

$\begin{aligned} C-\alpha_i-\mu_i&=0,\\ \mu_i\varepsilon_i&=0,\\ y_i(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})-(1-\varepsilon_i)&\geq0 \end{aligned}$

を順番に用いることで

$\begin{aligned} \alpha_i=0\Longrightarrow \mu_i=C\gt0\Longrightarrow \varepsilon=0\Longrightarrow y_i(\beta_0+{}^t\boldsymbol{x}_i\boldsymbol{\beta})\geq1 \end{aligned}$

が得られる。 $0\lt\alpha_i\lt C$ のとき、

$\begin{aligned} \mu_i\varepsilon_i&=0,\\ C-\alpha_i-\mu_i&=0 \end{aligned}$

を用いることで $\varepsilon_i=0$ を得る。更に $\alpha_i\{y_i(\beta_0+x_i\beta)-(1-\varepsilon_i)\}=0$ より

$\begin{aligned} 0\lt\alpha_i\lt C\Longrightarrow yi(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})-(1-\varepsilon)=0\Longrightarrow y_i(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})=1 \end{aligned}$

である。
　 $\alpha_i=C$ のとき $\varepsilon_i\geq0$ および $\alpha_i\{y_i(\beta_0+x_i\beta)-(1-\varepsilon_i)\}=0$ より

$\begin{aligned} \alpha_i=C\Longrightarrow yi(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})-(1-\varepsilon)=0\Longrightarrow y_i(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})\leq1 \end{aligned}$

が成り立つ。更に $\alpha_i\{y_i(\beta_0+x_i\beta)-(1-\varepsilon_i)\}=0$ より

$\begin{aligned} y_i(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})\lt1\Longrightarrow \alpha_i=0 \end{aligned}$

である。他方で、

$\begin{aligned} y_i(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})-(1-\varepsilon_i)&\geq0,\\ \mu_i\varepsilon_i&=0,\\ C-\alpha_i-\mu_i&=0 \end{aligned}$

を用いることで

$\begin{aligned} y_i(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})\lt1\Longrightarrow \varepsilon_i\gt0\Longrightarrow \mu_i=0\Longrightarrow \alpha_i=C \end{aligned}$

が成立する。　 $\blacksquare$ )

　このような性質を利用して、 $\boldsymbol{\beta}$ から $\beta_0$ を求める方法を検討する。
　まず $y_i(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})\neq1$ であるような $i$ が少なくとも1つ存在することを示す。上記の命題から、 $a_1=\cdots=a_n=0$ かつ $y_i(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})\gt1,$ $i=1,\cdots,n$ を意味する。この場合、

$\begin{aligned} L_P=\displaystyle{\frac{1}{2}\left\|\boldsymbol{\beta}\right\|_2}+C\displaystyle{\sum_{i=1}^{n}\varepsilon_i}-\displaystyle{\sum_{i=1}^{n}\alpha_i\left\{y_i(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})-(1-\varepsilon_i)\right\} }-\displaystyle{\sum_{i=1}^{n}\mu_i\varepsilon_i} \end{aligned}$

の最初の2項は小さく出来るから、 $\alpha_1,\cdots,\alpha_n$ が最適であることに矛盾する。
　次に $y_i(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})\neq1$ の下で $\alpha_i=C$ は $\varepsilon_i\gt0$ を、 $\alpha_i=0$ は $\varepsilon=0$ を意味する。そこで $\varepsilon_{*}=\displaystyle{\min_{i}\varepsilon_i}$ とおき、各 $\varepsilon_i$ を $\varepsilon_i-\varepsilon_{*}$ に、 $\beta_0$ を $\beta_0+y_i\varepsilon_{*}$ に置き換えても $y_i\pm1$ となるから、

$\begin{aligned} y_i(\beta_0+x_i\beta)-(1-\varepsilon_i)\geq0 \end{aligned}$

の等号が成立し、 $\varepsilon_i\gt0$ が $\varepsilon_i\geq0$ に置き換わるのみであるから、7つの $\mathrm{KKT}$ 条件は成立する。しかし $L_P$ は第3項が $\varepsilon_{*}\displaystyle{\sum_{i=1}^{n}\alpha_i}\gt0$ だけ減少するから、上記の命題から、少なくとも1つの $i$ について $y_i(\beta_0+{}^{t}\boldsymbol{x}_i\boldsymbol{\beta})=1$ を得る。

　具体的に $0\lt\alpha_i\lt C$ となるような $i$ を集めて、 $\beta_0=y_i-{}^{t}\boldsymbol{x}_i\boldsymbol{\beta}$ の算術平均を取ればよい。
　双対問題

$\begin{aligned} L_D=\displaystyle{\sum_{i=1}^{n}\alpha_i}-\displaystyle{\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_jy_iy_j\boldsymbol{x}_i{}^{t}\boldsymbol{x}_j },\\ \displaystyle{\sum_{i=1}^{n}\alpha_i y_i}=0,\\ 0\leq\alpha_i\leq C \end{aligned}$

を二次計画法の解ソルバーで解く。これには

$\begin{aligned} L_D&=-\displaystyle{\frac{1}{2}{}^{t}\boldsymbol{\alpha}D_{\mathrm{mat}}}\boldsymbol{\alpha}+{}^{t}\boldsymbol{d}_{\mathrm{vec}}\boldsymbol{\alpha},\\ A_{\mathrm{mat}}\boldsymbol{\alpha}&\geq\boldsymbol{b}_{\mathrm{vec}} \end{aligned}$

となるような $D_{\mathrm{mat}}\in\mathbb{R}^{n\times n},A_{\mathrm{mat}}\in\mathbb{R}^{m\times n}, \boldsymbol{d}_{\mathrm{vec}}\in\mathbb{R}^n,\boldsymbol{b}_{\mathrm{vec}}\in\mathbb{R}^m,m\geq1$ を指定する必要がある。

8.4　カーネルを用いた拡張

　サポートベクトルマシンを双対問題で解く必要があるのは、 $L_D$ が内積 $\left<\cdot,\cdot\right>$ を用いて

$\begin{aligned} L_D=\displaystyle{\sum_{i=1}^{n}\alpha_i}-\displaystyle{\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_j y_i y_j}\left< x_i,x_j\right> \end{aligned}$

と書ける。 $V$ を内積が定義されたベクトル空間だとして、 $\phi:\mathbb{R}^p\rightarrow V$ を用いると、内積 $\left< x_i,x_y\right>$ を $k(x_i,x_j)=\left<\phi(x_i),\phi(x_j)\right>$ に置き換えることができる。この場合、 $(\phi(x_1),y_1),\cdots,(\phi(x_n),y_n)$ から非線形の分類規則が得られる。
　 $(i,j)$ 成分が $\phi(x_i),\phi(x_j)\in V$ の内積であるような行列 $K$ によってカーネルを構成する。実際、任意の $\boldsymbol{z}\in\mathbb{R}^n$ に対して、二次形式

$\begin{aligned} {}^{t}\boldsymbol{z}K\boldsymbol{z}&=\displaystyle{\sum_{i=1}^{n}\sum_{j=1}^{n}z_i\left<\phi(x_i),\phi(x_j) \right> z_j}\\ &=\left< \displaystyle{\sum_{i=1}^{n}z_i\phi(x_i)},\displaystyle{\sum_{j=1}^{n}z_j\phi(x_j)}\right>\\ &=\left\|\displaystyle{\sum_{i=1}^{n}z_i\phi(x_i)}\right\|^2\geq0 \end{aligned}$

が非負定値で対称であるから、狭義のカーネルである。

　ある $\phi:\mathbb{R}^p\rightarrow V$ を用いて、 $\boldsymbol{x}_i\in\mathbb{R}^p,i=1,\cdots,n$ をすべて $\phi(\boldsymbol{x}_i)\in V$ に置き換えるものとする。したがって $\boldsymbol{\beta}\in\mathbb{R}^p$ は $\boldsymbol{\beta}=\displaystyle{\sum_{i=1}^{n}\alpha_i y_i\phi(\boldsymbol{x}_i)}\in V$ であり、 $L_D$ の定義の中の内積 $\left<\boldsymbol{x}_i,\boldsymbol{x}_j\right>$ は $\phi(\boldsymbol{x}_i)$ と $\phi(\boldsymbol{x}_j)$ の間の内積、すなわちカーネル $K(\boldsymbol{x}_i,\boldsymbol{x}_j)$ となる。

$\begin{aligned} \end{aligned}$

8.5　Rによるシミュレーション

8.5.1　サポートベクトルマシンの例

library("quadprog")  #二次計画法を解くパッケージ

svm_1 <- function(X, y, C){
  eps <- 10^(-3)
  n <- nrow(X)
  meq <- 1
  Dmat <- matrix(nrow = n, ncol = n)
  for(i in 1:n){
    for(j in 1:n){
      Dmat[i,j] <- sum(X[i,] * X[j,]) * y[i] * y[j]
    }
  }
  
  Dmat <- Dmat + eps * diag(n) # 対角成分に微小な値を加えることで正則にする
  dvec <- rep(1:n)
  Amat <- matrix(nrow = (2*n +1), ncol = n)
  Amat[1,] <- y
  Amat[2:(n+1), 1:n] <- -diag(n)
  Amat[(n+2):(2*n+1),1:n] <- diag(n)
  Amat <- t(Amat)
  
  bvec <- c(0, rep(-C,n), rep(0,n))
  
  alpha <- solve.QP(Dmat, dvec, Amat, bvec = bvec, meq=1)$solution
  beta <- drop((alpha * y) %*% X)
  index <- (1:n)[eps < alpha & alpha < C - eps]
  beta_0 <- mean(y[index] - X[index, ] %*% beta)
  return(list(beta = beta, beta_0 = beta_0))
}

####
library("ggplot2")

set.seed(1)

a <- rnorm(1)
b <- rnorm(1)
n <- 300

X <- matrix(rnorm(n*2), nrow = n, ncol = 2)
y <- sign(a * X[,1] + b * X[,2] + 0.1*rnorm(n))

df_test <- data.frame(X,y)
colnames(df_test) <- c("x1","x2","y")

qq <- svm_1(X,y,10)

# グラフ
g <- ggplot(df_test,aes(x = x1, y = x2, fill = y, colour = y))
g <- g + geom_point() + geom_abline(intercept = -qq$beta_0/qq$beta[1],slope = -qq$beta[1]/qq$beta[2])
g <- g + labs(fill = "",color="",x="x1", y = "x2")
g <- g + theme_classic()
g <- g + theme(plot.title = element_text(hjust = 0.5),legend.position = "none",
               legend.title=element_text(size = 12),
               legend.text=element_text(size = 12))
g <- g + scale_x_continuous(limits = c(-3,3),labels = seq(-3,3,1),breaks = seq(-3,3,1))
g <- g + scale_y_continuous(limits = c(-3,3),labels = seq(-3,3,1),breaks = seq(-3,3,1))
g <- g + geom_vline(xintercept = 0, linetype = 1)
g <- g + geom_vline(xintercept = seq(-3,3,1), linetype = 3)
g <- g + geom_hline(yintercept = 0, linetype = 1)
g <- g + geom_hline(yintercept = seq(-3,3,1), linetype = 3)

plot(g)

8.5.2　カーネルを用いたサポートベクトルマシンの拡張

svm_2 <- function(X, y, C, K) {
  eps <- 10^(-4)
  n <- nrow(X)
  Dmat <- matrix(nrow = n, ncol = n)
  for (i in 1:n){
    for(j in 1:n){
      Dmat[i, j] <- K(X[i, ], X[j, ]) * y[i] * y[j]
    }
  }
  Dmat <- Dmat + eps * diag(n)
  dvec <- rep(1, n)
  Amat <- matrix(nrow = (2 * n + 1), ncol = n)
  Amat[1, ] <- y
  Amat[2:(n + 1), 1:n] <- -diag(n)
  Amat[(n + 2):(2 * n + 1), 1:n] <- diag(n)
  Amat <- t(Amat)
  bvec <- c(0, -C * rep(1, n), rep(0, n))
  meq <- 1
  
  alpha <- solve.QP(Dmat, dvec, Amat, bvec = bvec, meq = 1)$solution
  index <- (1:n)[eps < alpha & alpha < C - eps]
  beta <- drop((alpha * y) %*% X)
  beta_0 <- mean(y[index] - X[index, ] %*% beta)
  return(list(alpha = alpha, beta_0 = beta_0))
}

K_linear <- function(x, y) {
  return(t(x) %*% y)
}

K_poly <- function(x, y) {
  return((1 + t(x) %*% y) ^ 2)
}

plot_kernel <- function(K, lty) {  # 引数ltyで，線の種類を指定する
  qq <- svm_2(X, y, 1, K)
  alpha <- qq$alpha
  beta_0 <- qq$beta_0
  
  f <- function(u, v) {
    x <- c(u, v)
    S <- beta_0
    for (i in 1:n)
      S <- S + alpha[i] * y[i] * K(X[i, ], x)
    return(S)
  }  # fは(x, y)でのz方向の高さを求める。これから輪郭を求めることができる。
  
  u <- seq(-2, 2, 0.1)
  v <- seq(-2, 2, 0.1)
  w <- array(dim = c(41, 41))
  for (i in 1:41)
    for (j in 1:41)
      w[i, j] <- f(u[i], v[j])
  contour(u, v, w, level = 0, add = TRUE, lty = lty)
}

# 実行
a <- 3
b <- -1
n <- 200
X <- matrix(rnorm(n * 2), ncol = 2, nrow = n)
y <- sign(a * X[, 1] + b * X[, 2] ^ 2 + 0.3 * rnorm(n))

plot(-3:3, -3:3, xlab = "X[, 1]", ylab = "X[, 2]", type = "n")
for (i in 1:n) {
  if (y[i] == 1) {
    points(X[i, 1], X[i, 2], col = "red")
  } else {
    points(X[i, 1], X[i, 2], col = "blue")
  }
}
plot_kernel(K_linear, 1)
plot_kernel(K_poly, 2)

8.6　Pythonによるシミュレーション

8.6.1　サポートベクトルマシンの例

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib
import scipy
from scipy import stats
from numpy.random import randn #正規乱数

import cvxopt
from cvxopt import matrix

a = randn(1)
b = randn(1)
n = 100
X = randn(n,2)
y = np.sign(a * X[:,0] + b * X[:,1] + 0.1*randn(n))
y = y.reshape(-1,1)

# 
def svm_1(X, y, C):
    eps = 0.0001
    n = X.shape[0]
    P = np.zeros((n,n))
    
    for i in range(n):
        for j in range(n):
            P[i,j] = np.dot(X[i,:], X[j,:]) * y[i] * y[j]
    
    P = matrix(P + np.eye(n) * eps)
    A = matrix(-y.T.astype(np.float))
    b = matrix(np.array([0]).astype(np.float))
    h = matrix(np.array([C]*n+[0]*n).reshape(-1,1).astype(np.float))
    G = matrix(np.concatenate([np.diag(np.ones(n)),np.diag(-np.ones(n))]))       
    q = matrix(np.array([-1]*n).astype(np.float))
    
    # ソルバーの実行
    res = cvxopt.solvers.qp(P,q,A=A, b=b,G=G, h=h)
    alpha = np.array(res['x'])
    beta = ((alpha*y).T@X).reshape(2,1)
    index = (eps < alpha[:, 0]) & (alpha[:, 0] < C - eps)
    beta_0 = np.mean(y[index]-X[index,:]@beta)
    return {'beta':beta, 'beta_0':beta_0}

# テスト
a=randn(1);b=randn(1)
n=100
X=randn(n,2)
y=np.sign(a*X[:,0]+b*X[:,1]+0.1*randn(n))
y=y.reshape(-1,1)
for i in range(n):
    if y[i]==1:
        plt.scatter(X[i,0],X[i,1],c="red")
    else :
        plt.scatter(X[i,0],X[i,1],c="blue")
res=svm_1(X,y,C=10)
def f(x):
    return -res['beta_0']/res['beta'][1]-x*res['beta'][0]/res['beta'][1]
x_seq=np.arange(-3,3,0.5)
plt.plot(x_seq,f(x_seq))
res

8.6.2　カーネルを用いたサポートベクトルマシンの拡張

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib
import scipy
from scipy import stats
from numpy.random import randn #正規乱数

import cvxopt
from cvxopt import matrix

# 
def K_linear(x,y):
    return x.T@y
def K_poly(x,y):
    return (1+x.T@y)**2

def svm_2(X,y,C,K):
    eps = 0.0001
    n = X.shape[0]
    P = np.zeros((n,n))
    for i in range(n):
        for j in range(n):
            P[i,j] = K(X[i,:],X[j,:]) * y[i] * y[j]
    #パッケージにあるmatrix関数を使って指定する必要がある
    P = matrix(P+np.eye(n)*eps)
    A = matrix(-y.T.astype(np.float))
    b = matrix(np.array([0]).astype(np.float))
    h = matrix(np.array([C]*n+[0]*n).reshape(-1,1).astype(np.float))
    G = matrix(np.concatenate([np.diag(np.ones(n)),np.diag(-np.ones(n))]))       
    q = matrix(np.array([-1]*n).astype(np.float))
    res = cvxopt.solvers.qp(P,q, A=A, b=b,G=G, h=h)
    alpha = np.array(res['x'])  #xが本文中のalphaに対応
    beta = ((alpha*y).T@X).reshape(2,1)
    index = (eps < alpha[:, 0]) & (alpha[:, 0] < C - eps)
    beta_0 = np.mean(y[index]-X[index,:]@beta)
    return {'alpha':alpha, 'beta':beta, 'beta_0':beta_0}

# テスト
a = 3
b = -1
n = 200
X = randn(n,2)
y = np.sign(a * X[:,0] + b * X[:,1]**2 + 0.3 * randn(n))
y = y.reshape(-1,1)

def plot_kernel(K,line):
    res = svm_2(X,y,1,K)
    alpha = res['alpha'][:,0]
    beta_0 = res['beta_0']
    def f(u,v):
        S = beta_0
        for i in range(X.shape[0]):
            S = S + alpha[i] * y[i] * K(X[i,:],[u,v])
        return S[0]
    uu = np.arange(-2,2,0.1)
    vv = np.arange(-2,2,0.1)
    ww = []
    for v in vv:
        w = []
        for u in uu:
            w.append(f(u,v))
        ww.append(w)
    plt.contour(uu,vv,ww,levels=0,linestyles=line)
for i in range(n):
    if y[i]==1:
        plt.scatter(X[i,0],X[i,1],c="red")
    else:
        plt.scatter(X[i,0],X[i,1],c="blue")
plot_kernel(K_poly,line="dashed")
plot_kernel(K_linear,line="solid")

*1:サポートベクトルマシンはこの仮定を満たしている。

前回

8. サポートベクトルマシン

8.2 最適化の理論

8.3 サポートベクトルマシンの解

8.4 カーネルを用いた拡張

8.5 Rによるシミュレーション

8.5.1 サポートベクトルマシンの例

8.5.2 カーネルを用いたサポートベクトルマシンの拡張

8.6 Pythonによるシミュレーション

8.6.1 サポートベクトルマシンの例

8.6.2 カーネルを用いたサポートベクトルマシンの拡張