R備忘録 - 記事一覧
- 投稿者: みゅ
- カテゴリ: なし
- 優先度: 普通
- 状態: 完了
- 日時: 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
- してdqrls2.c(114)を「src/appl」に置く.
- [2010.05.24]ランク落ちした場合に無限ループにはまっていたバグを修正.
- [2010.05.26]LAPACKの最小自乗法に対する、ストラテジーを理解.ランク落ちの場合も対応済み.
- 「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 ←
- 「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")
- 「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);
- 「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
・・・
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
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
$ 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))
> 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
> 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にした場合には、計算スピードの改善はほとんど見られない.
> 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
> 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",
- 先にノーマルな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
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)
> 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")
> 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
> (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
> (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(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
> (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版を使わないほうがよいかもしれません.
R備忘録 /状態空間モデリング/donlp2/その他のメモ