「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。目下、データ分析・語学に力点を置いています。今月(2022年10月)からは多忙につき、日々の投稿数を減らします。

MENU

統計的機械学習の数理100問(19/20)

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

を用いることにする。

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}^pf_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_iL_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:サポートベクトルマシンはこの仮定を満たしている。

プライバシーポリシー お問い合わせ