COLOR(magenta){SIZE(20){Rコード最適化のコツと実例集 (RjpWiki 誕生一周年記念企画)}}
// 間瀬：2004.06.05 開始
#contents
~
この tips (集)では R コードの最適化のコツと、実際例を紹介したいと思います。もちろん何が最適かは(実行時間という絶対的基準があるものの)多分に主観的要素が入る余地(実行環境にも依存するかもしれません)がありますが。
私見によれば、個別の問題に特有な工夫を別にしても、R で高速なコードを書く比較的小数のコツがあるような気がします。
R はインタプリタ言語ですから、C 等に比べれば実行速度は原理的に格段に遅くなります。しかし、実は R にはそうした欠点を補うための工夫がなされています。例えば
(1) 時間のかかる計算を内部的に C や FORTRAN サブルーチンを呼び出して高速化する(数値計算関係の関数のほとんどが該当します)、~
(2) 関数のベクトル化。R の関数のほとんどは引数としてベクトル(特別なベクトルとしての行列、配列を含む)をとることができ、返り値も対応するベクトルになります。こうした関数の内、汎用性が高い一部は内部的に C サブルーチンで処理することにより高速化をはかっています。
特に高速ベクトル化関数を意識した R コードを書くコツを身につければ、驚く程の高速化がしばしば可能になります。
なぜ高速化を意識したコードを書くことを心がける必要があるのでしょうか。実際、かなりの計算は高速化を特に意識しなくても間に合います。無理に高速化をはかったコードはしばしば理解しにくいものになり、書いた当人も一月も経つと何をしているの分からなくなりがちです(これは面倒でもコメントを丁寧に付け加えればかなり防げますが)。高速化を普段から心がけておく事が望ましい理由として、次をあげることができるでしょう。
- R の主要な利用目的にランダムシミュレーションがあります。特に流行りのMCMC(マルコフ連鎖モンテカルロ)法では、くり返し回数は時間の許す限り多く繰り返すのが基本です。(メトロポリス法を用いたある物理学のシミュレーション結果の論文の最後に、著者が「我々の結果は数台のワークステーションを一年間動かし続けて得られた」(言外に、悔しかったら2年間動かしてみろ)と誇らしげに書いていました。)こうした場合たとえ2倍の速度の向上でも大きな違いとなります。
- 統計計算では解を optim() 等の最適化関数を用いて数値的に求めるほかない事がしばしば起きます。こうした場合、ある関数を繰り返し評価する必要が起きます。したがって、一回あたりの計算時間を高速化しておかないと、いつまで経っても解が求められないことになります。
- 必要になった時だけ最適化したコードを書くことはできません。常日頃からそれを意識したコードを書く癖を付けておく必要があります。またそれにより、結果としてより洗練された R コードを書けるようになります(なる事もあるというべきですね)。
- R で素朴なプログラムを実行しても、多くの場合何とかなるという事実そのものが、1秒で済む計算を、わざわざ2秒かけて計算することが我慢できないという開発者たちのハッカー精神の賜物なのです。あの優れたメルセンヌツィスタ乱数も実行速度が遅ければ、誰も使わなかったでしょう。ボランティアによる高速発生プログラム(素人には何がなんだかさっぱり分かりませんが、それでも美しいコードであることだけはなんとなく分かります (^_^))の開発があったからこそ、R ユーザーはそれと意識すること無く享受できるわけです。
- これはまた R というプログラム言語の特性を理解することにつながります。
もちろん、問題に素直なコードを書くことが基本であることは当然です。しかし、更にそれを最適化したコードに書き換えて見ることをお勧めします。素直なコードと最適化されたコードは必ずしも矛盾しません。何が素直かとは、結局ユーザーの R プログラミングの理解度・熟達度に依存するからです。また、コードの簡潔さと高速性は必ずしも一致しませんが、(これまた独断的意見ですが)高速なコードはしばしば簡潔にならざるを得ません。
以下は私の個人的経験による R コードの高速化のコツをまとめたものです。一部内容が重複しています。これ以外にお気づきの点(当然あるでしょう反論も含め)があれば付け加えて下さい。
~
~
COLOR(red){注意}：system.time() 関数の返す5種類の時間のうち、最初が実働時間(CPU
タイム)、三番目がインタフェイス・表示等の時間を含んだ経過時間です。また、この時間は、全く同じ計算を実行しても、使用機種が違えば当然、同じ機械でも他の実行中のプログラム、現在使用可能なメモリ量等により相当違いますから、単なる目安と考えて下さい。本当は何回も実行した平均値を計算すべきでしょう。
~
* 意味のあるときは常に、意味の無いときも常にベクトル化を心がけよ [#befcb87d]
ベクトル化されたシステム関数を極力使うベクトル化演算に徹する。関数を書くときは、引数はベクトル、返り値もベクトルである関数を書くようにする。これはコードを簡単にするという大きなメリットもあります。但し、ベクトル化しさえすれば早くなるわけではありません。内部で返り値ベクトルの要素を個別に計算するような関数は、遅くなりがちです(見掛け上のベクトル化、これはこれで大きな価値がある)。速度の向上を伴うベクトル化は、高速なシステムベクトル化関数の「継目の無い連続的使用」によるベクトル化です。つまり、COLOR(red){ベクトル引数 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
* むき出しの繰り返しは避ける [#pb51b837]
(ほとんど)同じことを2回以上していれば改良の余地があります。コード中で同じことを繰り返している部分を高速化すれば、全体としての速度は確実に向上します。特に要注意は for, while, repeat ループです。どうしても避けられない事も多いですが、工夫で避けられることもまた多いです。プロファイル機能でコードの処理時間のボトルネックを探してみると、問題の箇所を発見できます。COLOR(red){ループを使うことをためらうな、但し最後の手段として。}
## 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 # 少し遅い
* 全体を操作せよ [#if173a98]
狭い部屋で作業するよりも、広い場所で仕事をするほうが当然能率的です。作業対象を一度にベクトル、行列、配列にまず放りこんでおいて、それから全体として加工すると見通しが良くなリます。数学でいう集合の考えに通じるところがあります。いっぽう R には特にベクトル、行列、配列を高速に処理する関数がいくつも用意されていますから、それらを組み合わせて(ここでも成分を意識した処理はタブー)処理すれば高速な処理が実現できます。そうした関数としては例えば、ベクトル、行列の四則演算、添字操作、積(%*%)、外積 (%o% 等)、周辺和(colSums(). colMeans()等)があります。COLOR(red){作業対象を全体として処理する事を心がける。}
## 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 倍ほどの速度になります。
* リスト変数は高価と心得る [#n169a2ea]
例外として、リスト処理があります。(私の理解不足かもしれませんが)リスト処理は必ずしも高速化されていないような気がします。もちろん、これは極めて便利なデータ構造ですが、高速化が問題になる処理の途中ではできるだけ避けた方が良いような(?)気がします。
> 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)
* 一つずつ沢山よりも、沢山を一度に [#f20ba136]
R は相当サイズの大きな変数を使ってもめげたりしません。必要と思われる以上(例えば数十倍)のデータを一度に生成し、一旦ベクトル等に格納し、最終的に必要なものだけを使って、後は気前良く捨て去っても処理時間の点で十分割にあうことがあります。例えば、乱数を一個ずつ生成し、順に処理するよりも、一度に生成してその中から、必要なだけ使う方が場合によれば早くなります(関数呼び出しのオーバヘッド等の理由で)。またそうすることにより、直前で述べた「全体として操作する」工夫が使えるようになります。COLOR(red){関数呼び出しはタダではない。}
> 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
* 贅沢は素敵だ(花森安治、若い世代には意味不明？) [#ha026d4f]
上の系になりますが、予めサイズが分かっている変数(ベクトル・配列・リスト)はそのサイズで始めに宣言しておく( x <- numeric(1000) 等)。予めサイズが分からない場合も、絶対これを越えることが無いというサイズで宣言しておくと代入操作が高速になります。R ではベクトル、行列、配列、リストは途中でサイズの拡大ができますが、これは言語的には無理をしているわけで、頻繁にこれをすると目立って効率が悪くなりますCOLOR(red){大名プログラミングに徹する、余っているメモリをけちる小市民プログラミングスタイルを止める。}
## ベクトルを段々伸ばす
> 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
* 論理判断は守銭奴のごとくけちる [#t4f9ae2f]
ループと並ぶ高速化の敵は論理判断です。どんな言語でもこれは例外ではありません。繰り返し操作と、論理判断はほとんどのアルゴリズムの中核ですから、全く無くすことは無理ですが、工夫により最小限におさえることはしばしば可能です。場合によれば全く無くすこともできます。
## 百万個の数のランダムな加減算(工夫で論理判断を掛け算で置き換え)
> 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 関数はベクトル化されていますから、(二者択一)論理判断＋繰り返しを同時に比較的高速に行うことができお勧めです。何よりもコードが簡単になります。COLOR(red){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))})
* ベクトルは添字で操作せず、(論理)添字集合で操作する [#kb44ad0e]
ベクトルの部分ベクトルを x[3] の如く、個々の成分で操作するのは以上の原理のほとんどに反するCOLOR(red){極悪非道}な仕打と考えて下さい。さすがに x[1:5] を得るのにx[i] を i=1,2,3,4,5 で求めるなどということをする方は入門初期の方を除けば少ないと思いますが、案外ベクトル(一般に配列)の一部分を取り出すのに、要素毎に取り出す人が多いような気がします(あなたのことではありません)。COLOR(red){ベクトルは添字で操作せず、(論理)添字集合で操作する}というコツを習得するとコードが簡略化されるだけでなく高速化されます(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
* 返り値 [#g2273a6d]
これは高速化という意味ではなく、コードの簡潔化という意味での最適化です。R のほとんどすべての関数は返り値を持ちます。一見返り値とは無縁に見える演算子も暗黙の返り値をもつように設計されています。この事実を使うと、簡潔なコードを書くことが出来ます。ただし、これは初心者には気持ちが悪く、また調子に乗って使い過ぎると確かにコードは短くなりますが、わかりやすさという意味での簡潔さは失われます。COLOR(red){常に返り値を意識する。}
# 代入し、そして和を求める
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
* 問題を変形する [#c091531e]
R に助けてもらうことばかり考えず、R を助けてあげることも考えるべきです。当初の問題が高速コーディング困難となれば、頭を使い高速なコーディングが可能な(つまりこれまで述べて来たようなテクニックが使いやすい)問題に変形します。本体で A を一生懸命に計算しても、最後に必要なのはそれから導かれる B にすぎないなら、最初から B を計算したら問題が早くできるかもしれません(トランプ問題を参照、ここでは絵柄だけが問題ですから、4種類の絵柄を配ると考えれば高速化のヒントになります。)。逆(特殊な場合よりも一般的な問題を考える)もまた真かも知れません。紙と鉛筆を使った一時間の計算が、R を使った一カ月の計算より早くて正確な結果を生むかもしれません。結局理論的計算が不可能・困難・面倒臭いから R を使うわけです。COLOR(red){R の高速化の究極的なコツは R ではなく、頭を使うこと。}それがいつでもできたら世話は無いわけですが。
* デバッグしやすさを最初から考慮する [#t18a8ffe]
経験によれば次の近似的不等式が成り立ちます：
デバッグに必要な時間 >= (プログラミングに要するトータルの時間)/c where 1 < c <= 2
ですから最初からデバッグしやすいようにコーディングすべきです。
これには決め手はありませんが、いくつかのコツはあります。
-多重ループは混乱しやすい。ループが閉じるところにコメントを置く
-(日本語の)コメントを随所に置く
-分かりやすい、そして他と紛らわしくない変数名を使う。特に局所的な変数名と大局的な変数名を区別しやすくする。いつも x,y,i,j を使っていませんか？局所変数名に .x 等の名前を使っている例を見たことがあります。
-永続付値演算子 <<- はしばしば便利ですが、何をしているのか理解していないとなかなか突き止められないバグの原因になることがあります。
-プログラムを一連の関数に分割する。特に補助的な(しかし本質的な)箇所を副関数化する。小さな単位はデバッグがしやすい。しかし、個々は正しくても、組み合わせて使う際に思わぬエラーが生じる可能性もある。
-プログラムが求める結果を確かに計算していることを確認できるようなテストデータを複数用意する。プログラムが動かないうちは、バグがあることを確信できますが、一旦動き出した後に真の困難が始まります。
COLOR(red){誤つは人の常、デバッグは世の常、そして結果の正しさは神のみぞ知る。}
* 人生という長期的スパンで最適化を考える [#f78fa9ae]
高速化はあくまで手段であって、最終目的ではありません。高速化にある程度なれて来ると高速化が自己目的化し、アクロバット的なコードをついつい書いてしまうようになりますが、他人には全く解読不能、自分でも明日には何をしているのか分からなくなってしまうこともありがちです(私自身の経験から)。コードは再利用できないと勿体ないです。COLOR(red){人生という長期的スパンで見れば、分かりやすさも高速化、つまりトータルの作業時間の節約、の重要な要素です。}
* 違いのわかる(気になる)人向け [#p210557e]
- 頻繁に使う定数は変数に代入し、変数として引用する。COLOR(red){チリも積もれば高速化。}
> 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
- 既存変数を再利用(変数の構造・サイズが変わる場合はどうかはわかりません)。COLOR(red){環境にやさしい資源リサイクル}
> 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
- 永続付値（他の環境中の変数への付値）は手間がかかる。COLOR(red){遠くの x よりも近くの x}
> 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
- 変数初期化のコツ。同一のベクトル、配列変数を全要素 0 で繰返し初期化する(良くあること)場合、既存の同じサイズの変数に 0 を掛ければ済む。x-x でも同じこと(少し早くなるとのコメント)。全要素 1 で初期化したければ、0*x + 1 とする。COLOR(red){0 は無ならず}
> 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
-行列・配列の添字に関するCOLOR(red){ループは順序が大事}、逆順にループすると早いらしい(コメント参照、thanks)
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
* 補遺、それでもまだ遅いと文句の多い人は [#z79a5f55]
- COLOR(blue){金持なら：} CPU 速度の早い機械を買う、メモリを増やす。64ビットワークステーションを買い、64 ビット版 R を構築する。並列計算を検討。商用ソフトの購入を考える。人を雇って高速化させる。
- COLOR(blue){金がなければ：} 線形関数数値ライブラリ BLAS, LAPACK の高速版である ATLAS, 後藤ライブラリを利用する(ものによっては数倍早くなる! 何でも掲示版を参照)。半年待つ。CPU 速度の早い機械の登場で確実に数倍早くなります。R のバージョンアップを待つ。すぐ使えるアドオンパッケージを探しまくる。RCC コンパイラを試す。RjpWiki に質問（ある程度簡単で一般的興味を引くものでないと無理かも）。
- COLOR(blue){環境に恵まれれていれば：} 周りの R wizard に頭を下げる。学生・部下に丸なげ(先生に丸投げなどふらちなことを考えてはいけません)。
- COLOR(blue){Tierney 氏のバイトコンパイル関数}を使ってみる（[[なんでも掲示板]]参照）。ものによっては数倍早くなります。
- COLOR(blue){プログラミングに堪能なら：} C, FORTRAN サブルーチンを書く。
- COLOR(blue){頭に自信があれば：} 計算機を使わず、頭を使う。統計・確率の本をもう一度勉強する。
- COLOR(blue){上のいずれにも該当しなければ：} さっさとあきらめる。
~
以上を参考に R プログラミングに精進して下さい。いつの日か、プログラム中の if 文、for 文、そして x[i] 等の表現を見ると、COLOR(blue){無意識に蝿叩きに手が伸びる}ようになれば、あなたは R プログラマとしてCOLOR(blue){最適化}されたことになります。
* コード最適化の具体例 [#x2ef9ff6]
以下コード最適化の具体例を初心者への参考にいくつか挙げたいと思います。一つの記事が結構長くなりそうなので、例毎に独立したページにし、以下にそれへのリンクを張りたいと思います。中級、上級ユーザーからの投稿を期待します。但し、問題の意味を理解するだけで疲れるようなものは避け、また簡単すぎ無い例を挙げられたらと思いますので、ご協力下さい。RjpWiki への投稿記事から無断引用をするかもしれませんが、お許し下さい。
** (1) r-help 記事より [[Rコードの最適化（初心者用メモ）]] [#ef8e1c29]
** (2) Q&A 記事より [[Rコード最適化例：空間点パターンのメディアン]] [#o809a138]
** (3) 「The R Book」 中の例を元にして [[Rコードの最適化例：トランプ]] [#b475743b]
** (4) 「工学のためのデータサイエンス入門」中の例を元にして [[Rコードの最適化例：クイックソート]] [#v561e8b5]
** (5) [[Rコードの最適化例：ベクトルはベクトルとして操作する]] [#w295054a]
** (6) [[Rコードの最適化例：配列はベクトルとして操作する]] [#rd0abcf9]
** (7) [[Rコードの最適化例：レスリー行列による個体群成長]] Ｑ＆Ａ(初級者コース)投稿記事を例として [#ma190f92]
** (8) [[Rコードの最適化例：混合正規乱数の発生コード]] 舟尾さんのコードを教材に（論理ベクトルのうまい使いかたの例 [#y8077e84]
** (9) [[Rコードの最適化例：行が同じ数かどうかの判定]] r-help 記事より [#uaf7ecf0]
** (10) [[Rコードの最適化例：行列のクロス積]] r-help 記事より(2004.10.07) [#c98c7ad4]
** (11) 行列の成分毎の積 r-help 記事より（2005.04.17） [#f70f9281]
行列と１列行列の成分毎の積(行列としての積でなく)は、行列とベクトルの成分毎の積とした方が早い。
> 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
* コメント、「俺ならこうする」等ありましたら以下にどうぞ [#t0c1895e]
-RjpWiki 一周年記念コンテンツにふさわしい内容ではないでしょうか！R に多少馴染んだあたりで疑問に思うことの答えが余すところ無く書かれています！素晴らしいです！ -- [[舟尾]] &new{2004-06-13 (日) 14:50:58};
-因みに論理判断が遅いのは、判断の結果次第で次に行う処理が変わり、CPUが「次に備える」ことができないから。キーワードは「パイプライン化」。（論理判断命令の実行が遅い訳ではない） -- [[flatline]] &new{2004-06-21 (月) 01:01:13};
-Rの持つ組み込みベクトル化機能を使うと速いのは、「メモリ内のデータ領域の10番目を探して取り、それを処理し、次に11番目を探して取り、それを処理」という方式ではなく「メモリ内で次（隣り）にあるデータを取り、それを処理」という方式が取れるため（と思う）。 -- [[flatline]] &new{2004-06-21 (月) 01:04:04};
-計算機の記憶領域は階層化されている。メモリは小さいがデータの出し入れが早い。ハードディスクは大きいが遅い。頻繁に利用するデータのサイズがメモリ内に収まるようにするのが上級者。メモリのキャッシュに収まるとなおよい。 -- [[flatline]] &new{2004-06-21 (月) 01:56:19};
-R のライブラリは FORTRAN で書かれているのだろうか。~
配列の格納法（メモリ上の配置）から考えると，C とは異なり，列優先のアクセスが速そうな気がする。
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
追試希望。 -- &new{2004-06-22 (火) 19:27:10};
-確かに，y <- 0*x は速いですね。y <- x-x の方がもう少しだけ速いようです。 -- &new{2004-06-23 (水) 10:41:05};
- 1で初期化する場合は0*x+1よりx^0が高速。さらに論理値で良ければx==xがもっと高速 -- [[aaa]] &new{2009-05-01 (金) 06:31:21};
- 行列の行や列に対する演算には、rowMeansなどを使う方法とapply関数を使う方法があるとよく紹介されますが、速度は随分違うみたいです。前者が圧倒的に高速。 -- [[mia]] &new{2010-04-02 (金) 19:49:57};
- そういうことは，オンラインヘルプにちゃんと書いてあります。These functions are equivalent to use of apply with FUN = mean or FUN = sum with appropriate margins, but are a lot faster. -- &new{2010-04-02 (金) 20:11:52};
- sample(1:n,m) やa<-1:n; sample(a,m)より sample.int(n,m)のほうが速い -- &new{2013-09-20 (金) 01:18:14};
- その差は，ノミの金玉くらいだ -- &new{2013-09-20 (金) 11:12:00};
> 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 の引数を適切に設定すべきということ。
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
#comment
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
