大規模データで
重回帰などのいわゆる多変量解析を行う際、
複数のマシンで計算を分担することで処理が早くなるらしい。
分担の際は、各マシンでデータの一部を持っておき、
途中まで計算を行い、最後に足し合わせるらしい。
ただ、なんでそんな事ができるのかが分からなかったので、
調べてみました。
http://machinelearning.wustl.edu/mlpapers/paper_files/NIPS2006_725.pdf
自分の理解が正しければ、考え方は以下の通りです。
段取りのイメージは以下の通りです。
これで一応理屈は分かったのですが、
実際に結果を出してみたくなったので一応プログラムを組んでみました。
自分はシステム周りのことは正直あまりわかっていないので、
一台のマシンでバラバラのテーブルをバラバラのまま計算して、
最後に足し合わせてみました。
つまり、「ひとり分散処理」をやってみました。
まずはサンプルデータの作成
library(MASS) mu<-c(0,0,0,0) Sigma<-matrix(NA,4,4) diag(Sigma)<-1 Sigma[lower.tri(Sigma)]<-c(0.4,0.3,0.1,0.2,0.3,0.6) Sigma[upper.tri(Sigma)]<-c(0.4,0.3,0.1,0.2,0.3,0.6) data<-mvrnorm(n=3000000,mu=mu,Sigma=Sigma,empirical=F) lst<-split(as.data.frame(data),1:3)#データを３つに分割
計算
#lst=リスト形式のデータ,num=従属変数になる列番号
mplm<-function(lst,num){
#デザイン行列を作成（元データから従属変数を除外して切片をつける）
dzmatlst<-split(rep(NA,length(lst)),1:length(lst),1)
for(k in 1:length(dzmatlst)){
dzmatlst[[k]] <- as.matrix( cbind(1,lst[[k]][,-num]) )
}
#XX
wklst<-dzmatlst
for(k in 1:length(wklst)){
wklst[[k]] <- t(wklst[[k]])%*%wklst[[k]]
}
xx <- 0
for(k in 1:length(wklst)){
xx <- xx + wklst[[k]]
}
xx<-solve(xx)
#XY
for(k in 1:length(dzmatlst)){
wklst[[k]] <- t(dzmatlst[[k]])%*%lst[[k]][,num]
}
xy <- 0
for(k in 1:length(wklst)){
xy <- xy + wklst[[k]]
}
#回帰係数
b<-xx%*%xy
#予測値
for(k in 1:length(dzmatlst)){
if(k==1){
dzmat<-dzmatlst[[k]]
}else{
dzmat<-rbind(dzmat,dzmatlst[[k]])
}
}
yhat<-dzmat%*%b
yhat<-c(yhat)
#決定係数
for(k in 1:length(lst)){
if(k==1){
y<-lst[[k]][,num]
}else{
y<-c(y,lst[[k]][,num])
}
}
Rsq<-cor(y,yhat)^2
Rsq<-c(Rsq)
#残差
e<-y-yhat
#回帰係数の標準誤差
Ssq<- (t(e)%*%e) / (nrow(dzmat)-ncol(dzmat))
se_b<-sqrt(diag(as.numeric(Ssq)*xx))
#出力結果
out<-list(t(b),se_b,Rsq,yhat,e)
names(out)<-c("b","se_b","Rsq","yhat","e")
return(out)
}試してみる
start<-proc.time()
res<-mplm(lst,4)
mptime<-proc.time()-start
> res$b
1 V1 V2 V3
[1,] 0.0003611049 -0.181937 0.2525769 0.6037954
> res$se_b
1 V1 V2 V3
0.0004394828 0.0004947975 0.0004817994 0.0004627307
> res$Rsq
[1] 0.4198231
> mptime
ユーザ システム 経過
5.26 0.50 5.83
#lm()と比べてみる
start<-proc.time()
kk<-lm(data[,4]~data[,1]+data[,2]+data[,3])
summary(kk)
time<-proc.time()-start
> summary(kk)
Call:
lm(formula = data[, 4] ~ data[, 1] + data[, 2] + data[, 3])
Residuals:
Min 1Q Median 3Q Max
-3.8623 -0.5132 -0.0006 0.5129 3.9257
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0003611 0.0004395 0.822 0.411
data[, 1] -0.1819370 0.0004948 -367.700 <2e-16 ***
data[, 2] 0.2525769 0.0004818 524.237 <2e-16 ***
data[, 3] 0.6037954 0.0004627 1304.853 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7612 on 2999996 degrees of freedom
Multiple R-squared: 0.4198, Adjusted R-squared: 0.4198
F-statistic: 7.236e+05 on 3 and 2999996 DF, p-value: < 2.2e-16
> time
ユーザ システム 経過
17.21 0.74 18.07パラメータ推定値は一致しました。
lm()より若干処理速度が速いです。
lm()ほど入力方法が親切でないし、最低限の出力しかしないからかな？
と思いました。
あと、一台のマシンですべての処理をしているので、
テーブルを細かく分ければ分けるほど、処理時間は長くなります（笑）
