![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
Rコード最適化のコツと実例集 (RjpWiki 誕生一周年記念企画)
この tips (集)では R コードの最適化のコツと、実際例を紹介したいと思います。もちろん何が最適かは(実行時間という絶対的基準があるものの)多分に主観的要素が入る余地(実行環境にも依存するかもしれません)がありますが。 私見によれば、個別の問題に特有な工夫を別にしても、R で高速なコードを書く比較的小数のコツがあるような気がします。
R はインタプリタ言語ですから、C 等に比べれば実行速度は原理的に格段に遅くなります。しかし、実は R にはそうした欠点を補うための工夫がなされています。例えば
(1) 時間のかかる計算を内部的に C や FORTRAN サブルーチンを呼び出して高速化する(数値計算関係の関数のほとんどが該当します)、
(2) 関数のベクトル化。R の関数のほとんどは引数としてベクトル(特別なベクトルとしての行列、配列を含む)をとることができ、返り値も対応するベクトルになります。こうした関数の内、汎用性が高い一部は内部的に C サブルーチンで処理することにより高速化をはかっています。
特に高速ベクトル化関数を意識した R コードを書くコツを身につければ、驚く程の高速化がしばしば可能になります。
なぜ高速化を意識したコードを書くことを心がける必要があるのでしょうか。実際、かなりの計算は高速化を特に意識しなくても間に合います。無理に高速化をはかったコードはしばしば理解しにくいものになり、書いた当人も一月も経つと何をしているの分からなくなりがちです(これは面倒でもコメントを丁寧に付け加えればかなり防げますが)。高速化を普段から心がけておく事が望ましい理由として、次をあげることができるでしょう。
もちろん、問題に素直なコードを書くことが基本であることは当然です。しかし、更にそれを最適化したコードに書き換えて見ることをお勧めします。素直なコードと最適化されたコードは必ずしも矛盾しません。何が素直かとは、結局ユーザーの R プログラミングの理解度・熟達度に依存するからです。また、コードの簡潔さと高速性は必ずしも一致しませんが、(これまた独断的意見ですが)高速なコードはしばしば簡潔にならざるを得ません。
以下は私の個人的経験による R コードの高速化のコツをまとめたものです。一部内容が重複しています。これ以外にお気づきの点(当然あるでしょう反論も含め)があれば付け加えて下さい。
注意：system.time() 関数の返す5種類の時間のうち、最初が実働時間(CPU
タイム)、三番目がインタフェイス・表示等の時間を含んだ経過時間です。また、この時間は、全く同じ計算を実行しても、使用機種が違えば当然、同じ機械でも他の実行中のプログラム、現在使用可能なメモリ量等により相当違いますから、単なる目安と考えて下さい。本当は何回も実行した平均値を計算すべきでしょう。
ベクトル化されたシステム関数を極力使うベクトル化演算に徹する。関数を書くときは、引数はベクトル、返り値もベクトルである関数を書くようにする。これはコードを簡単にするという大きなメリットもあります。但し、ベクトル化しさえすれば早くなるわけではありません。内部で返り値ベクトルの要素を個別に計算するような関数は、遅くなりがちです(見掛け上のベクトル化、これはこれで大きな価値がある)。速度の向上を伴うベクトル化は、高速なシステムベクトル化関数の「継目の無い連続的使用」によるベクトル化です。つまり、ベクトル引数 x がスカラーであるかのようにコーディングせよ。 その心は、x の成分を一々取り出して作業するなということで、線形代数の考えに通じるものがあります。
## exp() 関数はベクトル化されている
> x <- runif(1000000)
> test1 <- function () { # 百万個の数の指数を個別に計算
res <- numeric(1000000)
for (i in 1:1000000) res[i] <- exp(x[i])
res}
> system.time(test1())
[1] 9.27 0.04 18.51 0.00 0.00
> system.time(res <- exp(x)) # ベクトル化演算で一気に計算
[1] 0.15 0.03 0.27 0.00 0.00
(ほとんど)同じことを2回以上していれば改良の余地があります。コード中で同じことを繰り返している部分を高速化すれば、全体としての速度は確実に向上します。特に要注意は for, while, repeat ループです。どうしても避けられない事も多いですが、工夫で避けられることもまた多いです。プロファイル機能でコードの処理時間のボトルネックを探してみると、問題の箇所を発見できます。ループを使うことをためらうな、但し最後の手段として。
## 2 で割る操作を毎回している問題コード
> x <- runif(1000000)
> system.time({c <- 0; for (i in 1:100000) c <- c + x[i]/2; print(c)})
[1] 0.55 0.00 1.04 0.00 0.00
## 最後に割るだけで良い(し早い)
> system.time({c <- 0; for (i in 1:100000) c <- c + x[i]; print(c/2)})
[1] 0.46 0.00 0.88 0.00 0.00
## for loop 無しの例
> system.time(print(sum(x/2)))
[1] 0.07 0.02 0.16 0.00 0.00
## これでも最後に一回割る方が数倍早い
> system.time(print(sum(x)/2))
[1] 0.02 0.00 0.05 0.00 0.00
注意：同じ操作を副関数化することはコードの簡潔化に有効で常に推奨出来ますが、関数呼出しのオーバーヘッドを考えれば、速度の点ではむしろ遅くなります。
> test1 <- function () {
for (i in 1:1000) sum(runif(i))
}
> aux <- function (i) sum(runif(i)) # 副関数化
> test2 <- function () {
for (i in 1:1000) aux(i)
}
> c <- 0; for (i in 1:100) c <- c + system.time(test1()); c/100
[1] 0.1299 0.0000 0.1303 0.0000 0.0000
> c <- 0; for (i in 1:100) c <- c + system.time(test2()); c/100
[1] 0.1319 0.0000 0.1323 0.0000 0.0000 # 少し遅い
狭い部屋で作業するよりも、広い場所で仕事をするほうが当然能率的です。作業対象を一度にベクトル、行列、配列にまず放りこんでおいて、それから全体として加工すると見通しが良くなリます。数学でいう集合の考えに通じるところがあります。いっぽう R には特にベクトル、行列、配列を高速に処理する関数がいくつも用意されていますから、それらを組み合わせて(ここでも成分を意識した処理はタブー)処理すれば高速な処理が実現できます。そうした関数としては例えば、ベクトル、行列の四則演算、添字操作、積(%*%)、外積 (%o% 等)、周辺和(colSums(). colMeans()等)があります。作業対象を全体として処理する事を心がける。
## 1:100 一様乱数10個の平均値を1万個求める
> test1 <- function() {
res <- numeric(10000)
for (i in 1:10000) res[i] <- mean(sample(1:100, 10))
res
}
> system.time(test1())
[1] 1.87 0.03 3.94 0.00 0.00
> test2 <- function() {
y <- matrix(0, nr=10000, nc=10) # 乱数全体をまず行列に格納
for(i in 1:10000) y[i,] <- sample(1:100, 10)
res <- rowSums(y)/10 # 高速な行和関数を使う
}
> system.time(test2())
[1] 0.46 0.00 0.91 0.00 0.00
test1 も test2 も，1〜100までの非復元抽出による乱数ですね。普通の復元抽出による乱数の場合は，さらに高速になります。@rem:2010/05/20
test3 <- function()
{
y <- matrix(sample(1:100, 10000*10, replace=TRUE), 10000, 10) # 非復元抽出による 10 個の乱数を 10000 セット
rowMeans(y) # rowSums()/n は rowMeans()
}
system.time(test3())
test3 は，test2 の 40 倍ほどの速度になります。
例外として、リスト処理があります。(私の理解不足かもしれませんが)リスト処理は必ずしも高速化されていないような気がします。もちろん、これは極めて便利なデータ構造ですが、高速化が問題になる処理の途中ではできるだけ避けた方が良いような(?)気がします。
> x <- list(rep(0,1000000)) # リストとして処理 > system.time(for (i in 1000000) x[[i]] <- 0) [1] 0.22 0.02 0.24 0.00 0.00 # それなりに時間がかかる > x <- numeric(1000000) # ベクトルとして処理 > system.time(for (i in 1000000) x[i] <- 0) [1] 0 0 0 0 0 # 一瞬
list(rep(0,1000000)) により作られるリストに，x［［1］］ はあるが，x［［2］］ 以降はない。x［［1］］ は長さ 1000000
for (i in 1000000) x［［i］］ <- 0 が遅いのは，i >= 2 のとき，x［［i］］ を毎回作っているからに他ならない。@rem:2010/05/20
> n <- 1000000
> x <- list(rep(0,n)) # リストとして処理
> length(x)
[1] 1 # この時点では長さ 1 で，x[[1]][1]〜x[[1]][1000000] が存在する
> system.time(for (i in n) x[[i]] <- 0)
ユーザ システム 経過
0.01 0.00 0.01
> length(x)
[1] 1000000 # この時点では長さが 10000000
# x[[1]][1]〜x[[1]][1000000] と同時に
# x[[2]]〜x[[1000000]] が存在する要素は NULL
やっている内容が違う（しかも，非効率的）ので，時間がかかるのは当たり前。この後出てくる「R ではベクトル、行列、配列、リストは途中でサイズの拡大ができますが、これは言語的には無理をしているわけで、頻繁にこれをすると目立って効率が悪くなります」ということの実例になっています。
1000000 個の要素を持つリストを作るのは vector("list", 1000000) です。
以下の正しいプログラムでは，実行時間はベクトルを使う場合と同じくらい（一瞬）になります。
x <- vector("list", n)
system.time(for (i in n) x[[i]] <- 0)
R は相当サイズの大きな変数を使ってもめげたりしません。必要と思われる以上(例えば数十倍)のデータを一度に生成し、一旦ベクトル等に格納し、最終的に必要なものだけを使って、後は気前良く捨て去っても処理時間の点で十分割にあうことがあります。例えば、乱数を一個ずつ生成し、順に処理するよりも、一度に生成してその中から、必要なだけ使う方が場合によれば早くなります(関数呼び出しのオーバヘッド等の理由で)。またそうすることにより、直前で述べた「全体として操作する」工夫が使えるようになります。関数呼び出しはタダではない。
> x <- numeric(10000) ## 標準正規乱数, １万個を一つずつ > system.time(for (i in 1:10000) x[i] <- rnorm(1)) [1] 0.44 0.00 0.51 0.00 0.00 ## 標準正規乱数, １万個をいっぺんに > system.time(x <- rnorm(10000)) [1] 0.01 0.00 0.01 0.00 0.00
> x <- numeric(100000) ## ベルヌイ試行, 10万個を一つずつ > system.time(for (i in 1:100000) x[i] <- sample(c(0,1), 1, replace=TRUE)) [1] 3.87 0.00 3.89 0.00 0.0 ## ベルヌイ試行, 10万個を一回に > system.time(x <- sample(c(0,1), 100000, replace=TRUE)) [1] 0.01 0.01 0.03 0.00 0.00
## 関数呼出しはタダでないことを示す例
> test <- function() {} # 何もしない関数
> system.time(for(i in 1:100000) test()) # それでも呼出しにはそれなりの時間がかかる
[1] 0.16 0.00 0.27 0.00 0.00
上の系になりますが、予めサイズが分かっている変数(ベクトル・配列・リスト)はそのサイズで始めに宣言しておく( x <- numeric(1000) 等)。予めサイズが分からない場合も、絶対これを越えることが無いというサイズで宣言しておくと代入操作が高速になります。R ではベクトル、行列、配列、リストは途中でサイズの拡大ができますが、これは言語的には無理をしているわけで、頻繁にこれをすると目立って効率が悪くなります大名プログラミングに徹する、余っているメモリをけちる小市民プログラミングスタイルを止める。
## ベクトルを段々伸ばす > x <- numeric(0) > system.time(for (i in 1:10000) x[i] <- 0) [1] 1.68 0.01 1.69 0.00 0.00 ## 最初から必要なだけ入れ物を確保 > x <- numeric(10000) > system.time(for (i in 1:10000) x[i] <- 0) [1] 0.09 0.00 0.10 0.00 0.00
# 一様乱数の和が一万を越えるまでの必要個数
> test1 <- function() {
s <- c <- 1
while(s <= 10000) {s <- s + runif(1); c <- c + 1}
print(c(s,c))
}
> system.time(test1())
[1] 10000.14 20093.00
[1] 0.51 0.00 0.67 0.00 0.00
> test2 <- function() {
x <- cumsum(runif(40000)) # 気前良く4万個まず発生し累積和を計算
print(sum(x <= 10000) + 1) # x <- 10000 である要素の和を計算
}
> system.time(test2())
[1] 20048
[1] 0.02 0.00 0.01 0.00 0.00
## ベクトルを段々伸ばす例(空いているメモリを必死に探し回る R の苦労を想像すべきです)
> system.time({x <- NULL; for (i in 1:10000) x <- c(x,runif(1))})
[1] 1.74 0.00 3.42 0.00 0.00
## 最初に場所を確保
> system.time({x <- numeric(10000); for (i in 1:10000) x[i] <- runif(1)})
[1] 0.28 0.00 0.50 0.00 0.00
# 最初に十分大きなベクトル用意し、最後に余分な部分を捨てる方が早い例
test1 <- function() {
x <- numeric(0) # 空のベクトルから出発
for (i in 1:rnorm(1, m=3000, sd=200)) x[i] <- i
x
}
test2 <- function() {
x <- rep(NA,10000) # NA 値だけの十分大きなベクトルをあらかじめ用意
for (i in 1:rnorm(1, m=3000, sd=200)) x[i] <- i
na.omit(x) # NA 値部分を捨てる
}
system.time(x <- test1())
[1] 0.16 0.01 0.17 0.00 0.00
system.time(x <- test2())
[1] 0.03 0.00 0.03 0.00 0.00
length(x) # 実際必要だったベクトル長
[1] 3152
ループと並ぶ高速化の敵は論理判断です。どんな言語でもこれは例外ではありません。繰り返し操作と、論理判断はほとんどのアルゴリズムの中核ですから、全く無くすことは無理ですが、工夫により最小限におさえることはしばしば可能です。場合によれば全く無くすこともできます。
## 百万個の数のランダムな加減算(工夫で論理判断を掛け算で置き換え) > x <- 1:1000000 > y <- sample(c(TRUE, FALSE), 1000000, replace = TRUE) # ランダム論理ベクトル > s <- 0; system.time(for (i in 1:1000000) if (y[i]) s <- s + x[i] else s <- s - x[i]) [1] 5.99 0.01 6.01 0.00 0.00 > system.time(s <- sum((2*y-1)*x)) # 工夫で論値判断を無くす [1] 0.21 0.05 0.28 0.00 0.00
また abs() 関数等はいわば関数自身が論理判断を内蔵していますから、これを使って赤 裸々な論理判断を無くすことができます。
また if, while 文はベクトル化されていない稀な構文の例ですから、これを繰り返し使うことは高速化とは相容れません。一方、論理判断 =, <, <=, >, >=, そして ifelse 関数はベクトル化されていますから、(二者択一)論理判断＋繰り返しを同時に比較的高速に行うことができお勧めです。何よりもコードが簡単になります。if, while, repeat 文は悪徳消費者金融、最後の選択肢と心得る。
## 混合正規分布乱数一万個を生成
> x <- rnorm(10000) # N(0,1) 乱数一万個
> y <- rnorm(10000, sd = 2, mean= 1) # N(1,2^2) 乱数一万個
> w <- numeric(10000) # 結果をいれるベクトルを用意
> test1 <- function(){ # if 文を使う素朴なやり方
for (i in 1:10000)
if (runif(1)< 0.35) w[i] <- x[i] else w[i] <- y[i]
w
}
> system.time(test1())
[1] 0.51 0.00 0.52 0.00 0.00
> system.time(w <- ifelse(runif(10000)<0.35, x,y)) # ifelse 文を使えば
[1] 0.02 0.00 0.01 0.00 0.00
> system.time({p <- (runif(10000)<0.35); w <- p*x+(1-p)*y}) # 論理ベクトル p を使う
[1] 0.01 0.00 0.01 0.00 0.00
データの順序を問題にしないのなら，パラメータの異なる複数の母集団からのデータを c でつなぐという方法もある。
以下の方法は，上の例の掛け算を利用する方法（三番目）よりは少し遅いが，ifelse を利用する方法（二番目）よりは三倍程度速い。それぞれのデータの個数も，二項乱数で決めている。@rem:2010/05/20
system.time({r <- rbinom(1, n, 0.35); w3 <- c(rnorm(r), rnorm(n-r, sd = 2, mean = 1))})
ベクトルの部分ベクトルを x[3] の如く、個々の成分で操作するのは以上の原理のほとんどに反する極悪非道な仕打と考えて下さい。さすがに x[1:5] を得るのにx[i] を i=1,2,3,4,5 で求めるなどということをする方は入門初期の方を除けば少ないと思いますが、案外ベクトル(一般に配列)の一部分を取り出すのに、要素毎に取り出す人が多いような気がします(あなたのことではありません)。ベクトルは添字で操作せず、(論理)添字集合で操作するというコツを習得するとコードが簡略化されるだけでなく高速化されます(R はそういう使い方をするように設計されているからです。) いくつか簡単な例をあげておきます。本格的な例は Rコードの最適化例：クイックソート を見て下さい。
x[x < 0] # x の負の要素だけからなるベクトル x[0 < x & x < 1] # 区間 (0,1) に入る要素からなるベクトル x[x == 0] # 値が 0 に等しい要素からなるベクトル x[x != NA] # 欠損値を取り除いたベクトル x[x %% 3 == 0] # 3 の倍数からなる要素のベクトル x[y == 0] # y の 0 の要素に対応する x の要素からなるベクトル
残念なことに，欠損値が絡むと，上の例のはすべて（！）うまく動かない。
四番目以外は，要素に NA がなければちゃんと動く。しかし，NA を除くという四番目のものは要素に NA が含まれるのだから，絶対に思ったとおりには動かない。@rem:2010/05/20
> (x <- c(-2, 0, 3, NA, 9, 0.4, 5, 0)) [1] -2.0 0.0 3.0 NA 9.0 0.4 5.0 0.0 # 一番目 > x[x < 0] [1] -2 NA # 二番目 > x[0 < x & x < 1] [1] NA 0.4 # 三番目 > x[x == 0] [1] 0 NA 0 # 四番目 > x[x != NA] [1] NA NA NA NA NA NA NA NA # さすがに予想外の結果だろう # 五番目 > x[x %% 3 == 0] [1] 0 3 NA 9 0 # 六番目 > (y <- c(NA, 1, 0, 1, 0, 1, 0, 1)) [1] NA 1 0 1 0 1 0 1 > x[y == 0] [1] NA 3 9 5
欠損値があっても正しく機能するようにするには，subset 関数を使おう。
> (x <- c(-2, 0, 3, NA, 9, 0.4, 5, 0)) [1] -2.0 0.0 3.0 NA 9.0 0.4 5.0 0.0 > # 一番目 > subset(x, x < 0) [1] -2 > # 二番目 > subset(x, 0 < x & x < 1) [1] 0.4 > # 三番目 > subset(x, x == 0) [1] 0 0 > # 四番目 > subset(x, !is.na(x)) [1] -2.0 0.0 3.0 9.0 0.4 5.0 0.0 > # または > x[!is.na(x)] [1] -2.0 0.0 3.0 9.0 0.4 5.0 0.0 > # 五番目 > subset(x, x %% 3 == 0) [1] 0 3 9 0 > # 六番目 > (y <- c(NA, 1, 0, 1, 0, 1, 0, 1)) [1] NA 1 0 1 0 1 0 1 > subset(x, y == 0) [1] 3 9 5
教訓：欠損値があるときの挙動もちゃんと確認しておこう（実際のデータには欠損値がありがち）
次の性質も覚えておくと役にたつことがあります(参考例クイックソート参照)
x <- runif(10) y <- x[x < 0] # 長さ0のベクトル(空集合)になる all.equal(x, c(x,y)) # 長さ0のベクトルを付け加えても変化しない TRUE
# 論理値 TRUE, FALSE は数が要求される局面では整数 1,0 に強制変換される > TRUE + TRUE +FALSE # 1 + 1 + 0 と解釈される [1] 2 > sum(c(TRUE, TRUE, FALSE, TRUE)) [1] 3 > FALSE * 0.5 # 0 * 0.5 と解釈される [1] 0 > x <- 0:5 > x < 3 # 論理値ベクトルになる [1] TRUE TRUE TRUE FALSE FALSE FALSE > sum(x < 3) # 3 未満の x の要素の数 [1] 3
これは高速化という意味ではなく、コードの簡潔化という意味での最適化です。R のほとんどすべての関数は返り値を持ちます。一見返り値とは無縁に見える演算子も暗黙の返り値をもつように設計されています。この事実を使うと、簡潔なコードを書くことが出来ます。ただし、これは初心者には気持ちが悪く、また調子に乗って使い過ぎると確かにコードは短くなりますが、わかりやすさという意味での簡潔さは失われます。常に返り値を意識する。
# 代入し、そして和を求める x <- y z <- sum(x) # 一行で書けば(x に y を代入し、その暗黙の返り値、つまり代入された値自身の和を求める) z <- sum(x <- y)
y <- z
x <- exp(y)
# 一行で書けば
x <- exp(y <- z)
# 実行時間は少し遅くなるらしい
> system.time(for(i in 1:10000) {z <- runif(1000); y <- z^2; x <- exp(y)})
[1] 6.27 0.00 6.28 0.00 0.00
> system.time(for(i in 1:10000) {x <- exp(y <- (z <- runif(1000))^2)})
[1] 6.46 0.00 6.47 0.00 0.00
ほとんど誤差範囲で，むしろ後者の方が速いのでは？@rem:2010/05/20
> n <- 100000
> system.time(for(i in 1:n) {z <- runif(1000); y <- z^2; x <- exp(y)})
ユーザ システム 経過
9.788 0.815 10.558
> system.time(for(i in 1:n) {x <- exp(y <- (z <- runif(1000))^2)})
ユーザ システム 経過
9.715 0.827 10.472
R に助けてもらうことばかり考えず、R を助けてあげることも考えるべきです。当初の問題が高速コーディング困難となれば、頭を使い高速なコーディングが可能な(つまりこれまで述べて来たようなテクニックが使いやすい)問題に変形します。本体で A を一生懸命に計算しても、最後に必要なのはそれから導かれる B にすぎないなら、最初から B を計算したら問題が早くできるかもしれません(トランプ問題を参照、ここでは絵柄だけが問題ですから、4種類の絵柄を配ると考えれば高速化のヒントになります。)。逆(特殊な場合よりも一般的な問題を考える)もまた真かも知れません。紙と鉛筆を使った一時間の計算が、R を使った一カ月の計算より早くて正確な結果を生むかもしれません。結局理論的計算が不可能・困難・面倒臭いから R を使うわけです。R の高速化の究極的なコツは R ではなく、頭を使うこと。それがいつでもできたら世話は無いわけですが。
経験によれば次の近似的不等式が成り立ちます：
デバッグに必要な時間 >= (プログラミングに要するトータルの時間)/c where 1 < c <= 2
ですから最初からデバッグしやすいようにコーディングすべきです。 これには決め手はありませんが、いくつかのコツはあります。
誤つは人の常、デバッグは世の常、そして結果の正しさは神のみぞ知る。
高速化はあくまで手段であって、最終目的ではありません。高速化にある程度なれて来ると高速化が自己目的化し、アクロバット的なコードをついつい書いてしまうようになりますが、他人には全く解読不能、自分でも明日には何をしているのか分からなくなってしまうこともありがちです(私自身の経験から)。コードは再利用できないと勿体ないです。人生という長期的スパンで見れば、分かりやすさも高速化、つまりトータルの作業時間の節約、の重要な要素です。
> system.time(for (i in 1:1000000) sample(1:52)) [1] 41.95 0.31 84.71 0.00 0.00 > a <- 1:52 > system.time(for (i in 1:1000000) sample(a)) # 少しは早い [1] 37.29 0.00 74.63 0.00 0.00
> test1 <- function() {x <- runif(1000000); y <- x^2} # 変数 y は無駄
> system.time(test1())
[1] 0.62 0.04 0.67 0.00 0.00
> test2 <- function() {x <- runif(1000000); x <- x^2} # 変数 x を再利用
> system.time(test2())
[1] 0.48 0.03 0.50 0.00 0.00> test1 <- function(){x <- runif(10000)} # 自前の環境中の x へ代入
> test2 <- function(){x <<- runif(10000)} # 親環境中の x へ代入
> rm(x); c = 0; for(i in 1:100) c <- c + system.time(test1()); c/100
[1] 0.0015 0.0000 0.0020 0.0000 0.0000
> x
Error: Object "x" not found # 変数 x は関数終了時に消える
> x <- numeric(10000); c = 0; for(i in 1:100) c <- c + system.time(test2()); c/100
[1] 0.0018 0.0000 0.0021 0.0000 0.0000
> x[1:5] # 当然 x は存在
[1] 0.02132324 0.36266310 0.50504673 0.89662908 0.67079898> x=1:5 > system.time(for(i in 1:10^6) y <- numeric(5)) [1] 6.22 0.03 6.34 0.00 0.00 > system.time(for(i in 1:10^6) y <- 0*x) [1] 2.94 0.01 3.00 0.00 0.00 > system.time(for(i in 1:10^6) y <- x-x) [1] 2.58 0.01 2.75 0.00 0.00 > x=matrix(1:16, nc=4, nr=4) > system.time(for(i in 1:10^4) y <- matrix(0, nc=4, nr=4)) [1] 0.35 0.00 0.38 0.00 0.00 > system.time(for(i in 1:10^4) y <- 0*x) [1] 0.04 0.01 0.05 0.00 0.00 > system.time(for(i in 1:10^4) y <- x-x) [1] 0.04 0.00 0.04 0.00 0.00
x <- array(1:10^6, c(100, 100, 100))
ijk <- function() { # 添字の順にアクセス
for (i in 1:100) for (j in 1:100) for (k in 1:100) x[i,j,k]
}
ikj <- function() { # こんなことをする人はいないだろうが
for (i in 1:100) for (j in 1:100) for (k in 1:100) x[i,k,j]
}
kji <- function() { # 添字の逆順にアクセス
for (i in 1:100) for (j in 1:100) for (k in 1:100) x[k,j,i]
}
e <- function() { # 添字を全く使わない例
for (e in x) e
}
> system.time(ijk())
[1] 3.71 0.00 3.72 0.00 0.00
> system.time(ikj())
[1] 3.56 0.00 3.58 0.00 0.00
> system.time(kji())
[1] 3.26 0.00 3.27 0.00 0.00
> system.time(e())
[1] 0.26 0.00 0.28 0.00 0.00
以上を参考に R プログラミングに精進して下さい。いつの日か、プログラム中の if 文、for 文、そして x[i] 等の表現を見ると、無意識に蝿叩きに手が伸びるようになれば、あなたは R プログラマとして最適化されたことになります。
以下コード最適化の具体例を初心者への参考にいくつか挙げたいと思います。一つの記事が結構長くなりそうなので、例毎に独立したページにし、以下にそれへのリンクを張りたいと思います。中級、上級ユーザーからの投稿を期待します。但し、問題の意味を理解するだけで疲れるようなものは避け、また簡単すぎ無い例を挙げられたらと思いますので、ご協力下さい。RjpWiki への投稿記事から無断引用をするかもしれませんが、お許し下さい。
行列と１列行列の成分毎の積(行列としての積でなく)は、行列とベクトルの成分毎の積とした方が早い。
> mm <- matrix(1, 1000, 1000) > ee <- matrix(1:100,nc=1) > system.time(mm*ee[,], TRUE) # 行列 ee を掛ける [1] 0.26 0.02 0.28 NA NA > system.time(mm*c(ee), TRUE) # ee をベクトルとした上で掛ける [1] 0.07 0.00 0.07 NA NA
おもしろい現象がある。R を起動した後最初に実行するときと二回目以降で，実行時間に違いがある。初回は前者が5，6倍遅いが，二回目以降は差がなくなる。
なお，ee のサイズが 100 なのは間違いかな？@rem:2010/05/20
# 一回目
> mm <- matrix(1, 1000, 1000)
> ee <- matrix(1:100,nc=1)
> system.time(mm*ee[,], TRUE) # 行列 ee を掛ける
ユーザ システム 経過
0.033 0.006 0.039
> system.time(mm*c(ee), TRUE) # ee をベクトルとした上で掛ける
ユーザ システム 経過
0.006 0.006 0.013
# 二回目
> mm <- matrix(1, 1000, 1000)
> ee <- matrix(1:100,nc=1)
> system.time(mm*ee[,], TRUE) # 行列 ee を掛ける
ユーザ システム 経過
0.007 0.006 0.012
> system.time(mm*c(ee), TRUE) # ee をベクトルとした上で掛ける
ユーザ システム 経過
0.006 0.007 0.013
なお，このときの R のバージョンがいくつかによりますが（R-2.3.1 前後だろうけど），注釈には ee[,] は行列と書いてありますが，少なくとも R-2.11.0 では，ee[,] も c(ee) もベクトルです。
> class(ee[,]) # ee[,] は，ベクトル！！ [1] "integer" > class(c(ee)) [1] "integer" > class(ee) # ee は行列なので， [1] "matrix" > system.time(mm*ee, TRUE) # これは，エラーになる 以下にエラー mm * ee : 適切な配列ではありません Timing stopped at: 0 0 0
a <- function()
{
x <- matrix(1:100000, nr=20000, nc=5)
for (loop in 1:50) {
for (i in 1:20000) {
for (j in 1:5) {
x[i,j]
}
}
}
}
b <- function()
{
x <- matrix(1:100000, nr=20000, nc=5)
for (loop in 1:50) {
for (j in 1:5) {
for (i in 1:20000) {
x[i,j]
}
}
}
}
system.time(a())
system.time(b())
> system.time(a())
[1] 31.49 3.27 46.92 0.00 0.00
> system.time(b())
[1] 26.09 2.85 38.93 0.00 0.00
追試希望。 -- 2004-06-22 (火) 19:27:10> library(rbenchmark)
> benchmark(sample(1000, 1e6, replace=TRUE), sample.int(1000, 1e6, replace=TRUE))
test replications elapsed relative user.self sys.self
2 sample.int(1000, 1e+06, replace = TRUE) 100 1.146 1.000 1.050 0.101
1 sample(1000, 1e+06, replace = TRUE) 100 1.179 1.029 1.063 0.121
なぜそうなるかは，sample 関数を表示してみると分かる。つまり，sample は sample.int を呼んでいるのだ。
> sample
function (x, size, replace = FALSE, prob = NULL)
{
if (length(x) == 1L && is.numeric(x) && x >= 1) {
if (missing(size))
size <- x
sample.int(x, size, replace, prob)
}
else {
if (missing(size))
size <- length(x)
x[sample.int(length(x), size, replace, prob)]
}
}
else 節はどういうことかというと，第1引数にベクトルを指定したときの話。この場合には遅い。しかし，第1引数に sample(1:n,m) などと無様な使用をするときであって，sample(1:n,m) は sample(n,m) と同じ結果をもたらすので，if 節で処理される。そして，sample.int が使用されるというわけ。つまり，sample.int が速いのではなく，sample の引数を適切に設定すべき（1:n ではなく n）ということ（n； a positive number, the number of items to choose from. See ‘Details.’）。
> library(rbenchmark)
> benchmark(sample(1:1000, 1e6, replace=TRUE), sample(1000, 1e6, replace=TRUE), sample.int(1000, 1e6, replace=TRUE))
test replications elapsed relative user.self sys.self
3 sample.int(1000, 1e+06, replace = TRUE) 100 1.155 1.006 1.054 0.106
1 sample(1:1000, 1e+06, replace = TRUE) 100 2.330 2.030 1.866 0.481
2 sample(1000, 1e+06, replace = TRUE) 100 1.148 1.000 1.051 0.102
