引き続き『カーネル多変量解析』 赤穂昭太郎 (岩波書店)
大文字のXと小文字のxが出てきますが、大文字は入力データを意味し、
以下の議論において所与とする。
X=(X_1,\ X_2,\ ...,\ X_n)
復習
先だっての第1章では、カーネル関数を用いた回帰の例を挙げていました。
線形回帰と同様、データ(X,Y)に対してy=f(x)で近似しますが、
カーネル関数を使った近似では、このf(x)が特殊な形をしている。
回帰係数をw_1~nとして、
f(x)=\sum_{i=1}^n w_i\ k(X_i, x) \tag{0}
特徴量
2章ではまず、特徴量という概念が説明されています。
多項式近似のイメージを使って説明してみましょう。
set.seed(111)
N <- 100
x <- runif(N, -2, 2)
y <- 2*x^2+3*x+rnorm(N, 0, 2)
こんなデータ(X, Y)を考えます。
入力は一次元の変数Xで、観測量がY。
Xの一次式で表せる関数(線形)では、直線的なカンケイしか表せないので、
曲線の式y=f(x)には、次数dの多項式を考える。
f(x)=\sum_{m=1}^d α_m x^m
この多項式は、xの線形ではないですが、
推定すべきαについて線形になっている事に注意しておきます。
サンプルデータでは2次式(d=2)と考えて、下記を用いる。
y=α_1*x+α_2*x^2
そんな事は知っているよ！だから何だよ！（ーー；）
ええ。だから、ですね。
# 必要に応じて、事前に install.packages("rgl")
library("rgl")
x2 <- x^2
plot3d(cbind(x, x2, y), col="Red", size=6)
これをグリグリ回すと分かると思います。
「xとyのカンケイ」は、単純な線や面になっていませんので、面倒。
「xとx²とyのカンケイ」にすると、明らかに平面に乗っている。
線や平面や超平面は、線形和で書けるはずなので簡単に解ける。
1次元の入力Xから2次元の特徴量x, x²を抽出した空間でデータを評価する事で、より高精度な近似が(単純な線形回帰により)得られる。
x_k <- seq(-2,2, 0.1)
y_1 <- predict(lm(y~x), data.frame(x=x_k))
y_2 <- predict(lm(y~x2+x), data.frame(x=x_k, x2=x_k^2))
par(tcl=-0.2, mgp=c(1.5, 0.3, 0))
plot(x,y, pch=16, col="Red")
lines(x_k, y_1)
lines(x_k, y_2)
しつこいですが、d個の特徴量のベクトルをφと表すと、
\begin{eqnarray}
φ(x)&=&\big(φ_1(x)=x,\ φ_2(x)=x^2, ...,\ φ_d(x)=x^d \big) \\
α&=&\big(α_1,\ α_2,\ ...,α_d\big)
\end{eqnarray}
の元、
\begin{eqnarray}
f(x)&=&\sum_{m=1}^d α_m x^m=\sum_{m=1}^d\ α_m\ φ_m(x) \\
&=&α^Tφ(x)
\tag{1}
\end{eqnarray}
という事です(xの線形ではなくφの線形)。
特徴ベクトルからみたカーネル関数
データ集合Xに対し、カーネル関数は、含まれる2点X_i, X_jにおける特徴ベクトル同士の内積として定義できる。
上の例では、2つの特徴ベクトルx, x²を用いているので、
\begin{eqnarray}
k(X_i, X_j)&=&φ(X_i)^Tφ(X_j) \\
&=& \big( φ_1(X_i),\ φ_2(X_i) \big)^T\big( φ_1(X_j),\ φ_2(X_j) \big) \\
&=& X_iX_j+X^2_iX^2_j
\end{eqnarray}
Rで内積マトリックスKを書いてみると、
K_{ij}=k(X_i, X_j)
dat <- data.frame(x, x2)
K <- matrix(0, N, N)
for(i in 1:N){
for(j in 1:N){
K[i, j] <- sum(dat[i,]*dat[j,])
}}
image(1:N, 1:N, K, col=grey(1:100/100))
当然、対称行列ですね。
この行列は、半正定値性を満たすという重要な性質がありますが、ここでは踏み込まない。
さて、任意のxにおける入力Xに対するカーネル関数の積和は、
f(x)=\sum_{i=1}^n k(X_i,x)
y_k <- rep(0, length(x_k))
for(k in 1:length(x_k)){
for(i in 1:N){
y_k[k] <- y_k[k]+x[i]*x_k[k]+x[i]^2*x_k[k]^2
}
}
適切な重みw_iを用いる事で、データ(X,Y)を近似できそうだ。
\begin{eqnarray}
f(x)&=&\sum_{i=1}^n w_i\ k(X_i, x) \tag{2}\\
&=&\sum_iw_i\ φ(X_i)^T\ φ(x)\\
&=&\Big(\sum_iw_i\ φ(X_i)\Big)^T\ φ(x) \tag{3}\\
\end{eqnarray}
(0)と(2)が同じ形をしている事に注意。
実は、これは自明ではなく、証明が必要であるが、ひとまず割愛。
(1)(3)を見比べて、
\begin{eqnarray}
α&=&\sum_{i=1}^n w_i\ φ(X_i) \\
α_m&=&\sum_{i=1}^n w_i\ φ_m(X_i) \\
&=&\sum_{i=1}^n w_i\ X_i^m &(m=1,2,...,d)
\end{eqnarray}
と置く事で、
f(x)=\sum_{m=1}^d α_m x^m=\sum_{i=1}^n w_i\ k(X_i, x) \tag{4}\\
すなわち、多項式近似において次数d個の係数αを推定することは、
カーネル関数を用いてデータ数n個の係数wを推定する事と表せる。
実装
前回の議論は、今回のカーネル関数においてもまったく同様なので、(4)を用いた最小二乗近似を満たすwは、
w=K^{-1}Y
で求まるはず。
例によってsolveの問題があるので、逆行列はMASS::ginvを用いてムーアペンローズの疑似逆行列で代用。
library(MASS)
set.seed(111)
N <- 100
x <- runif(N, -2, 2)
x2 <- x^2
y <- 2*x^2+3*x+rnorm(N, 0, 2)
dat <- data.frame(x, x2)
K <- matrix(0, N, N)
for(i in 1:N){
for(j in 1:N){
K[i, j] <- sum(dat[i,]*dat[j,])
}}
w <- ginv(K) %*% y
x_k <- seq(-2, 2, 0.1)
y_k <- rep(0, length(x_k))
for(k in 1:length(x_k)){
for(i in 1:N){
y_k[k] <- y_k[k]+w[i]*(x[i]*x_k[k]+x[i]^2*x_k[k]^2)
}}
par(tcl=-0.2, mgp=c(1.5, 0.3, 0))
plot(x,y, col="Red", pch=16)
lines(x_k, y_k, lwd=5)
lines(x_k, y_2, lty=2, lwd=3, col="green") # 大昔に作った多項式近似
カーネル近似が黒、多項式が緑破線。
ふ〜。
for文の嵐。
enjoy!!