//このページが表示される方は、URLから「action=SOURCE&」を削除してみてください [[Ｒ備忘録 - 記事一覧]] !!!ソースを修正して、(g)lmの計算速度の改善をはかる *投稿者： みゅ *カテゴリ： なし *優先度： 普通 *状態： 完了 *日時： 2010年05月21日 10時45分33秒 //{{bugstate}} !!内容 *OS ： LINUX （CentOS 5.4） *業務で使用していると大きなデータの計算に多大な時間がかかる．ｑｒ分解を使用せずに行列演算だけで計算することもできないことはないが、やはり安定したｑｒ分解を使って計算したい． **説明変数が２個とか３個なら、今のままでも別に困らない．そんな方にはまったく無用の代物． *ATLASを使えるようにソースを修正してみる． **[LAPACKのストラテジーを理解する必要があるようで|http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=485] **ご意見・ご質問・バグのご連絡などは[別館|http://phoenixx.sakura.ne.jp/mt/R/2010/05/glm.html]にお願いします． **ご使用にあたっては、自己責任でお願いします． **ソースを修正して、nnet（ニューラルネットワーク）の計算速度の改善をはかる【[[BugTrack-Ｒ備忘録/59]]】もどうぞ． *glm関数は中でglm.fit関数を呼んでいる．で、glm.fit関数の中で、以下のようにdqrlsというFORTRANの関数を呼び出しています． fit <- .Fortran("dqrls", qr = x[good, ] * w, n = ngoodobs, p = nvars, y = w * z, ny = as.integer(1), tol = min(1e-07, control$epsilon/1000), coefficients = double(nvars), residuals = double(ngoodobs), effects = double(ngoodobs), rank = integer(1), pivot = 1:nvars, qraux = double(nvars), work = double(2 * nvars), PACKAGE = "base") *このdqrlsフォートラン関数を見ると以下のように、QR分解→最小二乗問題を解く、という手順を取っているが、これはLINPACKの関数． call dqrdc2(x,n,n,p,tol,k,qraux,jpvt,work) 20 call dqrsl(x,n,n,k,qraux,y(1,jj),rsd(1,jj),qty(1,jj), 1 b(1,jj),rsd(1,jj),rsd(1,jj),1110,info) *従って、これを（ATLAS化された）BLAS・LAPACKにすると早くなるはず． **結果445.671秒かかっていた計算が、270.034秒で終了！ !!手順 !ソースを解凍 tar xvf R-2.11.0.tar.gz cd R-2.11.0 !LAPACK用ルーティンを作成 *して{{ref dqrls2.c}}を「src/appl」に置く． **[2010.05.24]ランク落ちした場合に無限ループにはまっていたバグを修正． **[2010.05.26]LAPACKの最小自乗法に対する、ストラテジーを理解．ランク落ちの場合も対応済み． !Makefile.inを変更する *「src/appl/Makefile.in」のCのソースのリストに、上記Cファイルを加える． SOURCES_C = \ bakslv.c binning.c \ cpoly.c cumsum.c \ fft.c fmin.c integrate.c interv.c \ lbfgsb.c \ machar.c maxcol.c \ pretty.c \ rcont.c rowsum.c \ stem.c strsignif.c \ tabulate.c \ uncmin.c \ zeroin.c dqrls2.c ← !glm.fit関数の定義を修正する *「src/library/stats/R/glm.R」にglm.fit関数の定義があるのでこれを修正する． fit <- .Fortran("dqrls", qr = x[good, ] * w, n = ngoodobs, p = nvars, y = w * z, ny = as.integer(1), tol = min(1e-7, control$epsilon/1000), coefficients = double(nvars), residuals = double(ngoodobs), effects = double(ngoodobs), rank = integer(1), pivot = 1:nvars, qraux = double(nvars), work = double(2 * nvars), PACKAGE = "base") *これを以下のように修正する fit <- .C("dqrls4", qr = x[good, ] * w, n = ngoodobs, p = nvars, y = w * z, ny = as.integer(1), tol = min(1e-7, control$epsilon/1000), coefficients = double(nvars), residuals = double(ngoodobs), effects = double(ngoodobs), rank = integer(1), pivot = integer(nvars), qraux = double(nvars), PACKAGE = "base") !Applic.hを修正する *「src/include/R_ext/Applic.h」を修正する． *dqrlsの下あたりに追加する． void F77_NAME(dqrls)(double *x, int *n, int *p, double *y, int *ny, double *tol, double *b, double *rsd, double *qty, int *k, int *jpvt, double *qraux, double *work); void dqrls2(double *x, int *n, int *p, double *y, int *ny, double *tol, double *b, double *rsd, double *qty, int *k, int *jpvt, double *qraux); void dqrls3(double *x, int *n, int *p, double *y, int *ny, double *tol, double *b, double *rsd, double *qty, int *k, int *jpvt, double *qraux); void dqrls4(double *x, int *n, int *p, double *y, int *ny, double *tol, double *b, double *rsd, double *qty, int *k, int *jpvt, double *qraux); !registration.cを修正する *「src/main/registration.c」を以下のように修正する． ・・・ CDEF(R_rowsum), CDEF(stemleaf), #if 0 CDEF(str_signif), ・・・ *↓ ・・・ static R_NativePrimitiveArgType stemleaf_t[] = {REALSXP, INTSXP, REALSXP, INTSXP, REALSXP}; static R_NativePrimitiveArgType dqrls2_t[] = {REALSXP, INTSXP, INTSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP}; static R_NativePrimitiveArgType dqrls3_t[] = {REALSXP, INTSXP, INTSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP}; static R_NativePrimitiveArgType dqrls4_t[] = {REALSXP, INTSXP, INTSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP}; ・・・ CDEF(R_rowsum), CDEF(stemleaf), CDEF(dqrls2), CDEF(dqrls3), CDEF(dqrls4), #if 0 CDEF(str_signif), ・・・ !コンフィグ ./configure --with-blas='-L/usr/local/atlas/lib -lptf77blas -latlas -lpthread' --with-lapack='-L/usr/local/atlas/lib -llapack -lptcblas' --enable-R-shlib *上記ではなく、以下のようにすると、ATLAS化LAPACKを、libR.soとかの作成に使ってくれる． **まぁ、上のコンフィグではmakeが通らないし・・・ **ATLASをRで使う（Linux編）【[[BugTrack-Ｒ備忘録/63]]】で書いたように、「--with-lapack」オプションはRで動的にロードする動的共有ファイルの作成に使われる．したがって、ソースに直接書いたLAPACKの関数は、「エントリーが見つからない」とエラーが出る．したがって以下のように「--with-blas」オプションに記述することで、エラーを回避． ./configure --with-blas='-L/usr/local/atlas/lib -llapack -lptcblas -L/usr/local/atlas/lib -lptf77blas -latlas -lpthread' --with-lapack --enable-R-shlib *コンフィグの結果 R is now configured for i686-pc-linux-gnu Source directory: . Installation directory: /usr/local C compiler: gcc -std=gnu99 -g -O2 Fortran 77 compiler: gfortran -g -O2 C++ compiler: g++ -g -O2 Fortran 90/95 compiler: gfortran -g -O2 Obj-C compiler: Interfaces supported: X11 External libraries: readline, BLAS(ATLAS), LAPACK(in blas) Additional capabilities: PNG, JPEG, TIFF, NLS, cairo Options enabled: shared R library, R profiling, Java Recommended packages: yes !make *makeする make ・・・ gcc -std=gnu99 -I. -I../../src/include -I../../src/include -I/usr/local/include -DHAVE_CONFIG_H -fpic -g -O2 -c dqrls2.c -o dqrls2.o ・・・ ar cr libappl.a bakslv.o binning.o cpoly.o cumsum.o fft.o fmin.o integrate.o interv.o lbfgsb.o machar.o maxcol.o pretty.o rcont.o rowsum.o stem.o strsignif.o tabulate.o uncmin.o zeroin.o 【dqrls2.o】 ch2inv.o chol.o dchdc.o dpbfa.o dpbsl.o dpoco.o dpodi.o dpofa.o dposl.o dqrdc.o dqrdc2.o dqrls.o dqrsl.o dqrutl.o dsvdc.o dtrco.o dtrsl.o eigen.o ranlib libappl.a *このように、作成した「dqrls2.c」をコンパイルして「libappl.a」を作成している． $ nm src/appl/libappl.a | grep dqrls3 00000000 T dqrls3 *「libR.so」作成の様子． gcc -std=gnu99 -shared -L/usr/local/lib -o libR.so CConverters.o CommandLineArgs.o Rdynload.o Renviron.o RNG.o agrep.o apply.o arithmetic.o array.o attrib.o base.o bind.o builtin.o character.o coerce.o colors.o complex.o connections.o context.o cov.o cum.o dcf.o datetime.o debug.o deparse.o deriv.o devices.o dotcode.o dounzip.o dstruct.o duplicate.o engine.o envir.o errors.o eval.o format.o fourier.o gevents.o gram.o gram-ex.o gramRd.o graphics.o grep.o identical.o inlined.o inspect.o internet.o iosupport.o lapack.o list.o localecharset.o logic.o main.o mapply.o match.o memory.o model.o names.o objects.o optim.o optimize.o options.o par.o paste.o platform.o plot.o plot3d.o plotmath.o print.o printarray.o printvector.o printutils.o qsort.o random.o raw.o registration.o relop.o rlocale.o saveload.o scan.o seq.o serialize.o size.o sort.o source.o split.o sprintf.o startup.o subassign.o subscript.o subset.o summary.o sysutils.o unique.o util.o version.o vfonts.o xxxpr.o ../unix/Rembedded.o ../unix/libunix.a 【../appl/libappl.a】 ../nmath/libnmath.a ../extra/zlib/libz.a ../extra/bzip2/libbz2.a ../extrapcre/libpcre.a ../extra/tre/libtre.a ../extra/xz/liblzma.a -L/usr/local/atlas/lib -llapack -lptcblas -L/usr/local/atlas/lib -lptf77blas -latlas -lpthread -lgfortran -lm -lreadline -lncurses -ldl -lm *make check $ LANG=C LC_ALL=C make check **make check は、通らない．が、気にしない．計算手続きを変更したんだから当たり前． **先にノーマルなRでmake checkまで確認しておけば、OK． **LAPACKを使ったことで、rstudentやcooks.distanceの戻り値が微妙に異なるため？ !!パフォーマンスの確認 *確認用ソース # 大きな行列で計算速度の比較 set.seed(1) nn <- 10000 p <- 1000 xx <- matrix(rnorm(nn*(p-1)), nn, p-1) dim(xx) xx1 <- cbind(1,xx) dim(xx1) b <- runif(p) yy <- cbind(1,xx) %*% b + rnorm(nn,0,0.2) dat <- data.frame(yy=yy, xx) dim(dat) system.time(res0 <- glm(yy~., data=dat)) *元のR > system.time(res0 <- glm(yy~., data=dat)) ユーザ システム 経過 76.083 1.744 77.555 > system.time(res0 <- glm(yy~., data=dat)) ユーザ システム 経過 74.843 1.979 76.541 > system.time(res0 <- glm(yy~., data=dat)) ユーザ システム 経過 74.983 1.591 76.306 > system.time(res0 <- lm(yy~., data=dat)) # ←ちなみにlm ユーザ システム 経過 38.052 1.455 39.505 *修正したR > system.time(res0 <- glm(yy~., data=dat)) ユーザ システム 経過 41.534 2.005 34.705 > system.time(res0 <- lm(yy~., data=dat)) # ←ちなみにlm．計算速度は変わらない． ユーザ システム 経過 37.899 1.490 39.386 **glm関数はほぼ２分の１の時間で、計算が終わっている（激速）． **lm関数のほうは、修正していないので、パフォーマンスは同じ（修正すればもちろん激速）． **ちなみにｐ＝１００にした場合には、計算スピードの改善はほとんど見られない． !FORTRAN・Cだけの部分のパフォーマンス *元のLINPACKの関数を使用した版 > system.time( + fit <- .Fortran("dqrls", qr=cbind(1,xx), n=as.integer(nn), + p = as.integer(p), y = yy, ny = 1L, tol = 1e-07, + coefficients = double(p), + residuals = double(nn), effects = double(nn), + rank = integer(1), pivot = 1L:p, qraux = double(p), + work = double(2 * p), PACKAGE = "base") + ) ユーザ システム 経過 35.558 0.284 35.841 *LAPACK対応版 > system.time( + res5 <- .C("dqrls4", qr = cbind(1,xx), n = as.integer(nn), + p = as.integer(p), y = yy, ny = 1L, tol = 1e-07, + coefficients = double(p), + residuals = double(nn), effects = double(nn), + rank = integer(1), pivot = integer(p), qraux = double(p)) + ) ユーザ システム 経過 18.431 0.338 14.596 **ランク落ちがなければもっと早くなるんだが・・・ > system.time( + res4 <- .C("dqrls4", qr = cbind(1,xx), n = as.integer(nn), + p = as.integer(p), y = yy, ny = 1L, tol = 1e-07, + coefficients = double(p), + residuals = double(nn), effects = double(nn), + rank = integer(1), pivot = 1L:p, qraux = double(p)) + ) ユーザ システム 経過 13.511 0.344 5.637 **これは、今のところ未対応です．あしからず． *このdqrlsフォートラン関数を使用しているのは以下なので、ここを修正すれば、大きな行列の(ｇ)ｌｍ関連に関しては、パフォーマンスの改善がはかれる． src/library/stats/R/glm.R:232: fit <- .Fortran("dqrls", src/library/stats/R/lm.R:115: z <- .Fortran("dqrls", src/library/stats/R/lm.R:197: z <- .Fortran("dqrls", src/library/stats/R/lsfit.R:82: z <- .Fortran("dqrls", !!make checkが通らない件 *先にノーマルなRでmake checkまで確認しておけば、OK． *テスト用スクリプト ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14) trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69) group <- gl(2,10,20, labels=c("Ctl","Trt")) weight <- c(ctl, trt) (lm.D9 <- glm(weight ~ group)) *修正する前 > rstudent(lm.D9) 1 2 3 4 5 6 -1.33259060 0.82197468 0.21801317 1.71787892 -0.79706453 -0.62792250 7 8 9 10 11 12 0.20324545 -0.75058167 0.44085879 0.15898682 0.21949039 -0.73360797 13 14 15 16 17 18 -0.37071215 -1.70480945 1.97125395 -1.27995362 2.30780309 0.33799054 19 20 -0.50536729 0.04266149 > cooks.distance(lm.D9) 1 2 3 4 5 6 0.0945790812 0.0382244879 0.0027880720 0.1479169968 0.0360249862 0.0226675959 7 8 9 10 11 12 0.0024240341 0.0320765752 0.0113035039 0.0014846636 0.0028258760 0.0306862305 13 14 15 16 17 18 0.0080191438 0.1460022339 0.1860514935 0.0878987313 0.2385544132 0.0066750039 19 20 0.0148009406 0.0001070475 *修正後、このように微妙に異なる． > rstudent(lm.D9) 1 2 3 4 5 6 -1.33971090 0.83273377 0.21811584 1.71882586 -0.79745286 -0.62822424 7 8 9 10 11 12 0.20334114 -0.75094587 0.44106820 0.15906160 0.22306112 -0.74589420 13 14 15 16 17 18 -0.37677535 -1.73732737 2.01076013 -1.30277065 2.35734660 0.34351089 19 20 -0.51369116 0.04335358 > cooks.distance(lm.D9) 1 2 3 4 5 6 0.104750903 0.049114618 0.002816910 0.149446941 0.036397602 0.022902053 7 8 9 10 11 12 0.002449107 0.032408352 0.011420419 0.001500020 0.003872741 0.042054159 13 14 15 16 17 18 0.010989892 0.200089783 0.254975571 0.120461431 0.326928565 0.009147806 19 20 0.020284053 0.000146704 !そもそも、LINPACKとLAPACKでQR分解の値自体が異なる． *確認用Rソース ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14) trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69) group <- gl(2,10,20, labels=c("Ctl","Trt")) weight <- c(ctl, trt) x <- model.matrix(weight ~ group) *LINPACK > qr(x) $qr (Intercept) groupTrt 1 -4.4721360 -2.2360680 2 0.2236068 2.2360680 3 0.2236068 0.1827440 4 0.2236068 0.1827440 5 0.2236068 0.1827440 6 0.2236068 0.1827440 7 0.2236068 0.1827440 8 0.2236068 0.1827440 9 0.2236068 0.1827440 10 0.2236068 0.1827440 11 0.2236068 -0.2644696 12 0.2236068 -0.2644696 13 0.2236068 -0.2644696 14 0.2236068 -0.2644696 15 0.2236068 -0.2644696 16 0.2236068 -0.2644696 17 0.2236068 -0.2644696 18 0.2236068 -0.2644696 19 0.2236068 -0.2644696 20 0.2236068 -0.2644696 attr(,"assign") [1] 0 1 attr(,"contrasts") attr(,"contrasts")$group [1] "contr.treatment" $rank [1] 2 $qraux [1] 1.223607 1.182744 $pivot [1] 1 2 attr(,"class") *LAPACK > qr(x, LAPACK=T) $qr (Intercept) groupTrt 1 4.4721360 2.2360680 2 -0.2880072 2.2360680 3 -0.2880072 0.2236068 4 -0.2880072 0.2236068 5 -0.2880072 0.2236068 6 -0.2880072 0.2236068 7 -0.2880072 0.2236068 8 -0.2880072 0.2236068 9 -0.2880072 0.2236068 10 -0.2880072 0.2236068 11 -0.2880072 -0.1236068 12 -0.2880072 -0.1236068 13 -0.2880072 -0.1236068 14 -0.2880072 -0.1236068 15 -0.2880072 -0.1236068 16 -0.2880072 -0.1236068 17 -0.2880072 -0.1236068 18 -0.2880072 -0.1236068 19 -0.2880072 -0.1236068 20 -0.2880072 -0.1236068 attr(,"assign") [1] 0 1 attr(,"contrasts") attr(,"contrasts")$group [1] "contr.treatment" $rank [1] 2 $qraux [1] 0.7763932 1.2880072 $pivot [1] 1 2 attr(,"useLAPACK") [1] TRUE attr(,"class") [1] "qr" *でもこれは分解の方法が異なるだけで最小二乗解として求まるパラメータの値は同じ． > qr.coef(qr(x), weight) (Intercept) groupTrt 5.032 -0.371 > qr.coef(qr(x, LAPACK=T), weight) [1] 5.032 -0.371 !glm *元のｇｌｍ（LINPACK） > (lm.D9 <- glm(weight ~ group)) Call: glm(formula = weight ~ group) Coefficients: (Intercept) groupTrt 5.032 -0.371 Degrees of Freedom: 19 Total (i.e. Null); 18 Residual Null Deviance: 9.417 Residual Deviance: 8.729 AIC: 46.18 *LAPACK版 > (lm.D9 <- glm(weight ~ group)) Call: glm(formula = weight ~ group) Coefficients: (Intercept) groupTrt 5.032 -0.371 Degrees of Freedom: 19 Total (i.e. Null); 18 Residual Null Deviance: 9.417 Residual Deviance: 8.729 AIC: 46.18 **両方同じ !ｇｌｍ - ランク落ちしている場合 *正しい答えは出るが、表現方法が異なるので、出力される値・形式は異なる． *元のｇｌｍ（LINPACK） > (lm.D9 <- glm(weight ~ group+rev(group))) Call: glm(formula = weight ~ group + rev(group)) Coefficients: (Intercept) groupTrt rev(group)Trt 5.032 -0.371 NA Degrees of Freedom: 19 Total (i.e. Null); 18 Residual Null Deviance: 9.417 Residual Deviance: 8.729 AIC: 46.18 > predict(lm.D9) 1 2 3 4 5 6 7 8 9 10 11 12 13 5.032 5.032 5.032 5.032 5.032 5.032 5.032 5.032 5.032 5.032 4.661 4.661 4.661 14 15 16 17 18 19 20 4.661 4.661 4.661 4.661 4.661 4.661 4.661 *LAPACK版 > (lm.D9 <- glm(weight ~ group+rev(group))) Call: glm(formula = weight ~ group + rev(group)) Coefficients: (Intercept) groupTrt rev(group)Trt 4.661 NA 0.371 Degrees of Freedom: 19 Total (i.e. Null); 18 Residual Null Deviance: 9.417 Residual Deviance: 8.729 AIC: 46.18 > predict(lm.D9) 1 2 3 4 5 6 7 8 9 10 11 12 13 5.032 5.032 5.032 5.032 5.032 5.032 5.032 5.032 5.032 5.032 4.661 4.661 4.661 14 15 16 17 18 19 20 4.661 4.661 4.661 4.661 4.661 4.661 4.661 *LINPACKとLAPACKでpivotingの仕方が違うため．これがわからない人は、LAPACK版を使わないほうがよいかもしれません． !!コメント //{{comment}}