Rの配列/リストのランダムアクセスのパフォーマンスを調べてみた

よくRでは配列/リストに対してfor文+インデックスでアクセスしてはいけない、ということが言われますが、そうせざるを得ない場面にも時々遭遇します。少し気になってインデックスでアクセスした時のパフォーマンスの特性を調べて見ました。

通常、データ構造の意味で配列といえばランダムアクセスは定数時間、つまり O(1) で、リストのランダムアクセスは O(n) です。ところが、Rの場合は必ずしもそうなっていないようです

実験条件

2つのアクセスパターンで実験をしました。

■ アクセスパターン1
配列とリストの i 番目の成分に100万回 1 を加える。i=1,2000,4000,6000,…,20000 に対して実行時間を調べる。

■ アクセスパターン2
配列とリストの 1~i 番目の中のランダムに選んだ添字の成分に 1 を加える操作を100万回行う。i=1,2000,4000,6000,…,20000 に対して実行時間を調べる。

結果

2つのアクセスパターンに対する結果を表したのが次の図(使ったスクリプトは記事の一番最後):

縦軸は100万回のランダムアクセスに要した時間、横軸は上の実験条件の所で述べた i を示しています。

考察

はっきり言ってうまく実験がうまく制御されているかどうかについては自信がありませんが(定数分のオーバーヘッドがどこに起因しているのかよくわからない)、この実験の結果を真とすればアクセスパターン1では配列/リストともに定数時間であり、一方、アクセスパターン2では配列は定数時間ですが、リストの方は長さが増加するとアクセス時間が増加することになります。
この実験をする前は

  • リストと名のつくものは無条件で所要時間が線形に増えるだろう

と思っていたので少し新鮮な結果でした。

疑問・まとめ

  • リストの同じ添字にアクセスするのは定数時間なのは、最近アクセスした要素を何らかの方法でキャッシュしているから?
  • 配列に関して、パターン1よりもパターン2のほうが所要時間が短い理由はなぜ?

まあ、100万回程度のアクセスにこんなに時間がかかるならやっぱりできるだけ使わないほうがよいということでしょう(100万次元のベクトルの足し算などはこの実験よりもはるかに高速に行うことができる)。どうしても使う場合は、特にリストの場合はキャッシュが効くように同じ添字になるべく連続してアクセスすべし、ということを頭の片隅にいれておくことにします。

実験に使ったスクリプト

ntry <- 1000000 # Number of repeat

ary.len <- 20000

# Defining index for test
idx <- seq(0,ary.len,ary.len/10)
idx[1] <- 1

# Array access time --- Pattern 1
ary <- rep(0,ary.len)
durations.ary <- matrix(0,length(idx),5)
j <- 1
for(i in idx){
  n <- sample(1:i,ntry,T)
  start <- proc.time()
  for(k in 1:ntry)
    ary[i] <- ary[i] + 1
  durations.ary[j,] <- as.numeric(proc.time() - start)
  j <- j+1
  gc();gc();gc();
}

# List access time --- Pattern 1
lst <- list(0)
for(i in 1:ary.len)
  lst[[i]] <- 0
durations.lst <- matrix(0,length(idx),5)
j <- 1
for(i in idx){
  n <- sample(1:i,ntry,T)
  start <- proc.time()
  for(k in 1:ntry)
    lst[[i]] <- lst[[i]] + 1
  durations.lst[j,] <- as.numeric(proc.time() - start)
  j <- j+1
  gc();gc();gc();
}

tmp <- c(durations.ary[,3],durations.lst[,3])
par(mfrow=c(1,2))
plot(idx,durations.ary[,3],type="n",ylim=c(7,8),xlab="Index of array/list",ylab="Excution time")
lines(idx,durations.ary[,3],type="o",lw=2,col="steelblue")
lines(idx,durations.lst[,3],type="o",lw=2,col="pink")
legend("bottomright",legend=c("Array access time","List access time"),col=c("steelblue","pink"),pch=1,lw=2)
title("Access pattern 1")


# Array access time --- Pattern 2
ary <- rep(0,ary.len)
durations.ary2 <- matrix(0,length(idx),5)
j <- 1
for(i in idx){
  n <- sample(1:i,ntry,T)
  start <- proc.time()
  for(m in n){
    ary[m] <- ary[m] + 1
  }
  durations.ary2[j,] <- as.numeric(proc.time() - start)
  j <- j+1
  gc();gc();gc();
}

# List access time --- Pattern 2
lst <- list(0)
for(i in 1:ary.len)
  lst[[i]] <- 0
durations.lst2 <- matrix(0,length(idx),5)
j <- 1
for(i in idx){
  n <- sample(1:i,ntry,T)
  start <- proc.time()
  for(m in n){
    lst[[m]] <- lst[[m]] + 1
  }
  durations.lst2[j,] <- as.numeric(proc.time() - start)
  j <- j+1
  gc();gc();gc();
}


tmp <- c(durations.ary2[,3],durations.lst2[,3])
plot(idx,durations.ary2[,3],type="n",ylim=c(7,8),xlab="Index of array/list",ylab="Excution time")
lines(idx,durations.ary2[,3],type="o",lw=2,col="steelblue")
lines(idx,durations.lst2[,3],type="o",lw=2,col="pink")
legend("bottomright",legend=c("Array access time","List access time"),col=c("steelblue","pink"),pch=1,lw=2)
title("Access pattern 2")

Rで認証プロキシを乗り越える

Rではパッケージのインストールなどはネットワーク経由(CRAN)で行われます。認証プロキシが必要なネットワークの内部でRのパッケージをインストールする必要がある場合、それを乗り越えるためのいくつかの方法についてメモ。

–internet2 を使う

Rの実行時に –internet2 というオプションを指定する方法。この指定をすると windows の場合はOSのネットワークプロキシの設定を参照する。これはR+プロキシで検索するとよく出てくる方法だけど、認証プロキシには対応していないのでなにか方法を考えないといけない。

以前紹介したRubyのlocalproxyを使って、OSのプロキシの設定をlocalhostにしてしまえば –internet2 でもプロキシを通過可能。

この方法は非常に汎用性が高く、rubyのgemとか、その他ネットワーク越しのパッケージインストールでは威力を発揮する方法ですが、rubyを導入したり、OS起動時にスクリプトを起動するなど、初期設定はめんどくさいので、とりあえずRだけ認証プロキシを通過すればよい場合は次の方法がおすすめです。

http_proxy, http_proxy_user を指定

Rの実行時に http_proxy=http://proxy.hogehoge.com:PORT http_proxy_user=ask を指定すれば(指定の方法は –internet2 と同じ方法)パッケージインストールなどの際にプロキシユーザーとパスワードを入力するダイアログが出現します。これに入力するだけ。特に難しいことはないはずです。

Chinese Restaurant Process を発生させてみた

タイトルの通りです。単なる遊びです。Blackwell-McQueen Polya’s urn scheme とか表記されたりもするみたいです。

Chinese Restaurant Process (CRP)

よく(名前にもあるように)中国料理屋の例が出ますが、個人的には「初期の壺に連続的なスペクトルの色(虹色)を持つボールが入っているバージョンのポリアの壺」と考えるほうがわかりやすかったです。つまり、以前に記事にしたポリアの壺の色数が(可算とは限らない)無限集合に拡張されたもの、ということです。兎に角、CRP では以下のスキームで無限の長さの確率変数列を発生させます。

  1. \(G_0\) を適当な確率分布、\(\alpha\) を正の実数とします
  2. \(k\leftarrow1\)とする
  3. \(G_0\) から乱数(色)をひとつサンプリング
  4. その乱数を記録し、同じ色のボールを空の壺にいれる
  5. \(k\leftarrow k+1\)とする
  6. \(\alpha/(\alpha+k-1)\) の確率で \(G_0\) からサンプリングし、それ以外のときは壺からボールを一つ取り出してそれをサンプルとする
  7. サンプルと同じ色のボールを壺に入れる
  8. ステップ5に戻る

実際には無限に数列をつくることはできないので、10000個サンプルした段階でやめることにします。長さ10000の系列を15回生成して、その分布をプロットしてみます。\(G_0\) は標準正規分布、\(\alpha=3\) としました。

# Generating n random variables via CRP(alpha, G0)
rcrp<-function(n,alpha,G0){
  res<-numeric(n); res[1]<-G0(1);
  for(i in 2:n){
    p<-alpha/(alpha+i-1)
    res[i]<-ifelse(runif(1)<p,G0(1),sample(res[1:(i-1)],1))
  }
  class(res)<-"crp"
  res
}

# Simple simulation
par(mfrow=c(5,3)); par(mar=rep(1.5,4)); par(oma=c(1,1,1,1))

N <- 10000

for(i in 1:15){
  plot(table(rcrp(N,3,rnorm))/N,xlim=c(-3,3),ylim=c(0,0.4),ylab="",xaxt="n")
  axis(side=1,at=-2:2)
  curve(dnorm,add=T)
}

ここからわかること:

  • \(G_0\) は連続値をとる密度関数だが、得られる分布は離散的(実はこれは10000個で生成をやめずに、無限個生成しても離散的)
  • 10000個のサンプルのほとんどは数~十数のポイントに落ちる

いろいろ調べると(証明すると)わかること:

  • 実は得られた15個の離散分布はディリクレ過程(←ここでは説明していない)に従う「ランダム確率分布」である
  • (これは、もちろん有限で止めているという意味では正確な表現ではない)
  • 本当のディリクレ過程からのサンプルされる分布は(離散だけども)サポートは無限集合になる
  • CRPによって得られる確率変数列は過去に依存しているので、iidには見えない。しかし、実際にはこれらは得られた離散分布からのiid とみなすことができる

共役事前分布について

定義など少し混乱気味になって、調べてみたので書いてみました。

定義

パラメトリックな確率分布の族 \(\{p(\theta); \theta\in\Theta\}\) を考える。確率変数 \(X\) の従う分布が \(p(x|\theta)\) のとき(もしくはそのように仮定するとき)、ベイズの定理より事後分布は

\begin{align*}
p(\theta|x)\propto p(x|\theta)p(\theta)
\end{align*}

となるが、このとき事前分布 \(p(\theta)\) と事後分布 \(p(\theta|x)\) が同じ分布族に属するときに \(p(\theta)\) は \(p(x|\theta)\) の共役事前分布であるという。同じ分布族に属するので、事前パラメータ \(\theta_0\) とデータ \(x\) が決定されると事後パラメータ \(\tilde\theta\) は \(\theta_0, x\) の関数として\(\tilde\theta(\theta_0, x)\) と書ける。
“共役事前分布について”の続きを読む

ディリクレ分布に関するメモ (6)

ディリクレ分布に関するメモ (3) ではガンマ乱数を正規化することでディリクレ分布に従う乱数を発生させました。

今回は別の方法を試してみます。

ポリアの壺

ポリアの壺(Polya’s urn)というのは以下のような手続きのことです。

  1. n種類の色のボールが入った壺を考える
  2. 壺から一つボールを取り出す
  3. 取り出したボールを壺に戻し、同じ色のボールを壺に追加する
  4. ステップ2に戻る

この手続きによって、壺の中のボールは増加し続けます。その色のシェアはランダムに変動しますが、やがて収束します。この収束した時のシェアを \(p^*\) とします。この極限のシェア \(p^*\) はボールを取り出した経路に依存し、これ自体が確率変数です。これは、例えば壺に赤と青のボールが一つずつ入っていたとして、最初の数回でたまたま赤が連続して出ると、その後も赤が出やすくなり、極限では赤のシェアが90%ということも起こりえるわけです。

さて、n色のポリアの壺を考えます。最初に \(\alpha=(\alpha_1,\cdots,\alpha_n)\) 個のボールが壺に入っていたとすると、このポリアの壺の極限分布が \(\text{Dirichlet}(\alpha)\) に従います(ちなみに、これは \(\alpha\) が自然数でなくても成立します)。
“ディリクレ分布に関するメモ (6)”の続きを読む

数学ビデオ “Dimensions” を見た


知人から数学に関する Dimensions という DVD をもらったので見た。クリエイティブ・コモンズなので、いろいろな動画サイトでも見れるようですが。

幾何学と代数学の関係に関して、CGを用いて説明したおよそ二時間の動画。

9つのチャプターから構成されていて、1つが13分程度。
“数学ビデオ “Dimensions” を見た”の続きを読む