読者です 読者をやめる 読者になる 読者になる

RのコードをCに書き変えてみたら1500倍くらい早くなったけど、早すぎて将来が不安

C R

昨日の続きのところ。Cでナダラヤワトン推定量*1計算させて、Rのoptimを使って計算したら早くなったっぽい気がした。でも、一変数はRだけでまともなコードがなかった(文字列evalとかのやつしか残っていなかった)ので速度比較ができなかった。

これじゃ最適化関数もCで書かないといけないのか、そこまでしなくていいいのかの判断材料がないので、Rで書いて比較することにした。Cのプログラムは昨日のところのまま。

dyn.load("/tmp/hoge.so")
cv <- function(y,x1,h) {
.C("cv",
   as.double(y), 
   as.double(x1),
   as.integer(length(y)),
   as.double(h),
   result = 0)$result 
}

m <- function(y,x,x1,h) {
.C("m",
   as.double(y),
   as.double(x), 
   as.double(x1),
   as.integer(length(y)),
   as.double(h),
   result = 0)$result 
}

epanechnikov <- function(X1,h1){
  X1 <- X1
  return(function(x1){
    u <- (x1 - X1)/h1
    if(1 > abs(u)){
      return(0.75 * (1-u^2))
    }else{
      return(0)
    }
  })
}

y.exp <- function(X1,Y,h1){
  X1 <- X1
  Y <- Y
  return(function(x1){
    u <- (x1 - X1)/h1
    if(1 > abs(u)){
      return(0.75 * (1-u^2) * Y)
    }else{
      return(0)
    }
  })
}

n <- 1000
s <- seq(-30,30,length.out=n)
y <- (function(x){sin(x)})(s) + rnorm(n,0,0.1)
plot(s,y)
(h <- optimize(function(h){cv(y,s,h)},c(0.1,0.4))$minimum)
lines(s,mapply(function(x){m(y,s,x,h)},s))
lines(lowess(s,y),col="red")

mean.h <- function(n){
  X1 <- seq(-3,3,length.out=n)
  Y <- (function(x){sin(x)})(X1) + rnorm(n,0,0.1)
  min <- 0.001
  max <- 1
  t_c <- system.time(h_c <- optimize(function(h){cv(Y,X1,h)},c(min,max))$minimum)[3]
  t_r <- system.time(h_r <- optimize(function(h1){
    sum(mapply(function(X1i,Yi){
      (Yi -
       (Reduce("+",Map(function(f){f(X1i)},mapply(function(X1,Y){y.exp(X1,Y,h1)},X1,Y))) - 0.75*Yi) /
       (Reduce("+",Map(function(f){f(X1i)},mapply(function(X1){epanechnikov(X1,h1)},X1))) - 0.75)
       )^2
    },X1,Y))
  },c(min,max))$minimum)[3]
  cat(h_c,"\t",h_r,fill=TRUE)
  time <- c(t_c,t_r)
  names(time) <- NULL
  return(time)
}

計算時間の短縮

プログラム間違ってるんじゃないかと思ったので、Cでのクロスバリデーションを最小にするhとRでのhとを出力させている。5回の計算時間の平均を求めるとこんな感じに。

> apply(mapply(mean.h,rep(100,5)),1,mean)
0.1893555 	 0.1893555
0.2037649 	 0.2037649
0.2211508 	 0.2211508
0.1818008 	 0.1818008
0.2865864 	 0.2865864
[1] 0.0040 5.5206
> apply(mapply(mean.h,rep(200,5)),1,mean)
0.1615935 	 0.1615935
0.1903541 	 0.1903541
0.3099866 	 0.3099866
0.3015350 	 0.3015350
0.2411889 	 0.2411889
[1]  0.0134 22.2692
> apply(mapply(mean.h,rep(400,5)),1,mean)
0.2658077 	 0.2658077
0.2255306 	 0.2255306
0.2060307 	 0.2060307
0.2256248 	 0.2256248
0.2548975 	 0.2548975
[1]   0.0572 106.5370
> apply(mapply(mean.h,rep(1000,5)),1,mean)
0.2054734 	 0.2054734
0.1681855 	 0.1681855
0.1848514 	 0.1848514
0.1892796 	 0.1892796
0.1982217 	 0.1982217
[1]   0.3398 769.6374
サンプル数 n=100 n=200 n=400 n=1000
Cを使った場合の計算時間 0.0040 0.0134 0.0572 0.3398
Rを使った場合の計算時間 5.5206 22.2692 106.5370 769.6374
計算時間比(R/C) 1380倍 1661倍 1862倍 2264倍

もちろん同じデータを使っているので、最適解も同じになるはずで上のように同じになっている。そんでもって驚くべきは計算時間がアホみたいに早くなっているということだ。1000倍以上早くなっているし、データ数が多くなるほど、比が大きくなっている。

なんかうまく行きすぎてて、たぶんどっかでミスってると思うんだけどどうなのか。最適化をするための関数を自分で用意するのは非常に手間がかかる。実装するとなると、最適化手法は初期点の問題や収束性の問題があるから一つの手法だけではなく、いくつか作って試すべきだ(最初は初期点に依存しないような手法であたりを付けて、それを元に鋭い刃物を使う、というような)。そうなると、勉強しないといけないことが増えまくるし、自分が作ったものはバグが入る可能性が大きい(笑…ちゃいけないけど)。

最適化してあげたい関数の中身だけCで書いてあげて(これは昨日のコードを見ると分かるけど、ものすごーく簡単)、信頼のおけるRの最適化関数にあとは投げられて、しかもRだけの時よりものすごく早い。こんなおいしい話があるときには注意しなさいよ、と小学校の頃から言われているんだけど、どうなのか。コンパイルするとそんなに早くなるものなのか?とりあえず色んな人に聞いてみよう。

*1:leave-one outしてるので、正確には違うのだが