トップ 差分 一覧 ソース 検索 ヘルプ RSS ログイン

BugTrack-R備忘録/64

R備忘録 /状態空間モデリング/donlp2/その他のメモ

R備忘録 - 記事一覧

ソースを修正して、(g)lmの計算速度の改善をはかる

  • 投稿者: みゅ
  • カテゴリ: なし
  • 優先度: 普通
  • 状態: 完了
  • 日時: 2010年05月21日 10時45分33秒

内容

  • OS : LINUX (CentOS 5.4)
  • 業務で使用していると大きなデータの計算に多大な時間がかかる.qr分解を使用せずに行列演算だけで計算することもできないことはないが、やはり安定したqr分解を使って計算したい.
    • 説明変数が2個とか3個なら、今のままでも別に困らない.そんな方にはまったく無用の代物.
  • ATLASを使えるようにソースを修正してみる.
  • 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用ルーティンを作成

  • してdqrls2.c(114)を「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-R備忘録/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関数はほぼ2分の1の時間で、計算が終わっている(激速).
    • lm関数のほうは、修正していないので、パフォーマンスは同じ(修正すればもちろん激速).
    • ちなみにp=100にした場合には、計算スピードの改善はほとんど見られない.

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フォートラン関数を使用しているのは以下なので、ここを修正すれば、大きな行列の(g)lm関連に関しては、パフォーマンスの改善がはかれる.
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

  • 元の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
    • 両方同じ

glm - ランク落ちしている場合

  • 正しい答えは出るが、表現方法が異なるので、出力される値・形式は異なる.
  • 元のglm(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版を使わないほうがよいかもしれません.

コメント