※ 引述《ror (回血375)》之銘言：
: [問題敘述]:
: 請簡略描述你所要做的事情，或是這個程式的目的
: 各位前輩好
: 想請教一個問題
: 目前我手上有大概幾百萬(快到千萬)組資料 一組資料是兩筆內容
: 現在要對每一組資料跑lm()回歸後 取出相關係數與P-value 等等資訊
: 但因為要跑的量實在是太大量了 導致速度很慢 已經使用parallel
: 目前寫法是 loop-> lm() ->summary-> 取出value
: 想請問是否有方式避免重複宣告lm()物件 ...或是其他可以提升速度的方式
: 感謝 Orz
以下可能跟原PO的問題關係小了一點XD
我提供一個利用RcppArmadillo and RcppParallel做
leave-one-out cross validation (LOOCV) 的例子
十萬個樣本做LOOCV:
程式：http://pastebin.com/4JK6VSd7
在我的電腦上跑一次花了173.79秒 (平均一次 0.17 milliseconds)
這個應該遠比用R的迴圈快很多...
記憶體使用的部分也保持在300MB以下，基本上不太會用到太多記憶體做暫存
給一個簡單的比較：
mses = vector('numeric', N)
st = proc.time()
for (i in 1:N){
coef_lm = coef(lm.fit(cbind(1, X[-i,]), y[-i]))
mses[i] = y[i] - c(1, X[i,]) %*% coef_lm
}
proc.time() - st
上面的程式在N = 5000大概花10秒，Rcpp的版本只需要0.3秒
題外話，Rcpp真的值得一學，如果想要直接使用R的BLAS可以考慮RcppArmadillo
如果覺得內建的BLAS太慢可以考用RcppEigen (他本身的BLAS也算是夠快的)
不過最簡單的方式就是用RRO，或是在linux/mac上也可以考慮openblas
最後再搭配RcppParallel就可以得到很大的加速了XDD
PS: 裡面有一段程式是避免RcppParallel跟RRO使用的MKL用到過多的thread
而造成速度的拖累，所以寫了一段去減少MKL使用的thread
在有超執行緒(HT)的電腦上，MKL會用2個，其餘則是使用1個thread
PS2: RRO是Revolution R Open
最後，原PO的問題也並不複雜，相關係數跟p-value可以分別透過
RcppArmadillo:::cor跟Rcpp:::pt做處理，如果有任何問題，歡迎討論XD
Rcpp:::pt用法：
http://stackoverflow.com/questions/20144528/how-use-correctly-rcpppt
補充，如果要用K-fold可以用下面的函數去生成cvIndex:
(剛好有現成的，順便分享XD)
好讀版：http://pastebin.com/RhJx3kLs
uvec cvfold_index_f(const uword n, const uword fold)
{
uword fold_n = n / fold, rem = n - fold_n * fold;
uvec tmp_vec = linspace<uvec>(0, fold - 1, fold),
index = vectorise(repmat(tmp_vec, fold_n, 1));
if (rem > 0)
index = join_cols(index, index.head(rem));
// random permutation
int j;
for (int i = 0; i < n; i++)
{
j = as_scalar(randi<uvec>(1, distr_param(i, n - 1)));
index.swap_rows(i, j);
}
return index;
}