中古マンション売買データを分析してみた(2.3) - Stanの行列演算速度を軽く試す
とりあえずStanになれるというのも兼ねて,軽くマニュアル読みながら*1,どういう書き方が良い書き方なのかを学ぶことにしました.そうしたら,マニュアルのp141 Vectorizationのところで,基本的にvectorやmatrix使った方が計算早いよって書いてあるので,実際にコード書いてみて試してみました.
サンプルコード
ベタな回帰式の書き方
data { int<lower=1> N; int<lower=1, upper=5> d_d[N]; int<lower=1900, upper=2014> d_f[N]; int<lower=1, upper=20> d_r[N]; int<lower=1, upper=20> d_s[N]; int<lower=0, upper=200> d_p[N]; } parameters { real a; real b_d; real b_f; real b_r; real b_s; real b_p; real<lower=0> s; } model { for (i in 1:N) d_p[i] ~ normal(a + b_d*d_d[i] + b_f*d_f[i] + b_r*d_r[i] + b_s*d_s[i], s); }
マトリックス演算
data { int<lower=1> N; int<lower=1> M; matrix[N, M] X; vector[N] Y; } parameters { real alpha; vector[M] beta; real<lower=0> sigma; } model { for(i in 1:N) Y[i] ~ normal(alpha + X[i] * beta, sigma); }
Rコード
キックする側のコードはこちら
# load library library('rstan') # load data d <- read.delim('data/mantions.csv', header=T, sep=',') d = na.omit(d) attach(d) # store elapsed time time_matrix = array() time_linear = array() # matrix operation X = t(rbind(distance, from, room, space)) Y = price d.stan = list(N=nrow(X), M=ncol(X), X=X, Y=Y) ## matrix operation fit.matrix = stan(file="script/mcmc_basic_lm_vector.stan", data=d.stan, iter=10, chains=1) for (i in 1:100) { time_matrix[i] = system.time( basic_lm.fit<-stan(fit=fit.matrix, data=d.stan, iter=100, chains=1, thin=10, warmup=20) ) } # usual calculation using linear expression d.stan = list(N=nrow(d), d_d=d$distance, d_f=d$from, d_r=d$room, d_s=d$space, d_p=d$price) fit.linear<-stan(file="script/mcmc_basic_lm.stan", data=d.stan, iter=10, chains=1) for (i in 1:100) { time_linear[i] = system.time( basic_lm.fit<-stan(fit=fit.linear, data=d.stan, iter=100, chains=1, thin=10, warmup=20) ) }
結果
結論からいうと,全くそんなことはなかったので,どういうことなのか未だに良くわかっていないんですよね... まぁある程度変数が増えると,マトリックス使わないと変数管理が面倒でやってられない感じもするので,基本的にはマトリックスでいった方がいいのかなぁとは思っていますが...
> summary(time_matrix) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.526 1.542 2.222 42.330 40.490 254.600 > summary(time_linear) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.219 0.614 0.933 22.770 46.710 112.900
追記
このエントリについて@beroberoさんから,以下のようなご指摘をいただきました.おそらく,この程度の行列演算だと高速化の恩恵はほとんどなくて,むしろ通常計算時に値の範囲を絞っているほうが効果がある可能性はたしかにありそうかなぁと思われます.
@smrmkt 行列演算はC++に変換された時にfor文とEigenライブラリの勝負になります。Eigenの中でもfor文回しているような場合は速度はほぼ同じと思います。大きな行列の積になると有名なアルゴリズムがあるので差がでそうです。http://t.co/wzu9SZTS7P
— berobero (@berobero11) 2014, 7月 20
そして,何より階層ベイズで個体差を考えた瞬間に,そもそも行列演算はできないので,結局このエントリの趣旨とは何だったのか状態だったり...
@berobero11 なるほど,どういう行列演算をするかによってくるわけですね.とはいえ,個体差を入れた瞬間にfor文使わないといけないことに気づきました.少なくとも階層ベイズ系の処理では行列演算は使えなさそうですね...
— 眼鏡をかけた猫背のネコ (@smrmkt) 2014, 7月 20