6

コードを最適化しているところですが、いくつか問題が発生しています。R での最大の高速化は、ループを使用する代わりにコードをベクトル化することから得られることを知っています。しかし、リストにデータがあり、コードをベクトル化できるかどうかわかりません。apply関数(lapply、 など)を使用してみましたvapplyが、これらの関数はよりクリーンなコードを記述するためのものであり、実際には内部でループを使用していることを読みました!

私のコードの 3 つの最大のボトルネックを次に示しますが、最初の部分については何もできないと思います。

1) データの読み込み

私は 277x349 の次元の 1000 行列のバッチを扱っています。doMCこれが私のスクリプトの最大のボトルネックですが、関数で複数のコアを利用するパッケージを使用することで、問題を少し緩和しましたforeach。これにより、1000 個の 277x349 マトリックスを含むリストが作成されます。

質問の目的のために、次元が 277 x 349 の 1000 個の行列のリストがあるとします。

# Fake data
l <- list()
for(i in 1:1000) {
  l[[i]] <- matrix(rnorm(277*349), nrow=277, ncol=349)
}

2) ボトルネック #1

(同じ次元の)いくつかの参照マトリックスと比較する必要があります。これにより、リスト内の 1000 の行列を参照行列と比較して、1000 の距離のベクトルを取得します。行列が同じ次元であることがわかっている場合、このステップをベクトル化できますか?

ここにいくつかのコードがあります:

# The reference matrix
r <- matrix(rnorm(277*349), nrow=277, ncol=349)
# The number of non NA values in matrix. Do not need to worry about this...
K <- 277*349

# Make a function to calculate distances
distance <- function(xi, xj, K, na.rm=TRUE) {
  sqrt(sum((xi - xj)^2, na.rm=na.rm)/K)
}

# Get a vector containing all the distances
d <- vapply(l, distance, c(0), xj=r, K=K)

このステップは を使用するとかなり高速ですvapplyが、コードの中で 3 番目に遅い部分です。

3) ボトルネック #2

ここで、参照行列に「最も近い」J 個の行列の加重平均行列を作成したいと考えています。d[1] < d[2] < ... < d[1000](並べ替えの手順がありますが、簡単にするためにそれを想定しています)。J=1,2,...,1000のときの加重平均行列を取得したい

# Get the weighted matrix
weightedMatrix <- function(listOfData, distances, J) {
  # Calculate weights:
  w <- d[1:J]^{-2} / sum(d[1:J]^{-2})

  # Get the weighted average matrix
  # *** I use a loop here ***
  x_bar <- matrix(0, nrow=nrow(listOfData[[1]]), ncol=ncol(listOfData[[1]]))
  for(i in 1:J) {
    x_bar <- x_bar + {listOfData[[i]] * w[i]}
  }

  return(x_bar)
}

# Oh no! Another loop...
res <- list()
for(i in 1:length(l) ) {
  res[[i]] <- weightedMatrix(l, d, J=i)
}

私は少し困惑しています。行列のリストに対する操作をベクトル化する直感的な方法がわかりません。

私が書いているスクリプトはかなりの頻度で呼び出されるため、少し改善するだけでも効果があります。


編集:

RE: 1) データの読み取り

私のデータは特別な形式であることを忘れていたので、R でデータを読み取るには特別なデータ読み取り関数を使用する必要があります。ファイルはnetcdf4形式でありnc_open、パッケージの関数を使用ncdf4してファイルにアクセスしています。 、そして、ncvar_get関数を使用して対象の変数を読み取る必要があります。良い点は、ファイル内のデータをディスクから読み取ることができ、そのデータをメモリに読み取っncvar_getて R で操作できることです。

foreachそうは言っても、行列のサイズと行列の数はわかっていますが、並列計算を可能にする関数が並列化されたループの結果を出力するため、データのリストで質問しました。リスト。foreach関数を使用すると、データの読み取りステップが約 3 倍速くなることがわかりました。

後でデータを 3 次元配列として配置できると思いますが、もしかしたら 3 次元配列の割り当てにかかる時間の方が節約できるよりも時間がかかるのではないでしょうか? 明日試してみる必要があります。


編集2:

ここに、スクリプトのタイミングの一部を示します。

元のスクリプト:

[1] "Reading data to memory"
user  system elapsed 
176.063  44.070  26.611 

[1] "Calculating Distances"
user  system elapsed 
2.312   0.000   2.308 

[1] "Calculating the best 333 weighted matrices"
user  system elapsed 
63.697  28.495   9.092 

これまでに次の改善を行いました: (1) データを読み取る前にリストを事前に割り当てます。(2) Martin Morgan の提案に従って、加重行列の計算を改善しました。

[1] "Reading data to memory"
user  system elapsed 
192.448  38.578  27.872 

[1] "Calculating Distances"
user  system elapsed 
2.324   0.000   2.326 

[1] "Calculating all 1000 weighted matrices"
user  system elapsed 
1.376   0.000   1.374 

いくつかのメモ:

foreachループで 12 個のコアを使用してデータを読み込みます ( registerDoMC(12))。スクリプト全体の実行には、改善の前後で約 40 秒 / 36 秒かかります。

ボトルネック #2 のタイミングはかなり改善されました。以前は、加重行列の上位 3 分の 1 (つまり 333) のみを計算していましたが、スクリプトはすべての加重行列を元の時間の何分の 1 かで計算できるようになりました。

助けてくれてありがとう。後でコードを微調整して、リストの代わりに 3D 配列で動作するようにスクリプトを変更できるかどうかを確認します。計算が確実に機能することを確認するために、しばらく時間をかけて計算を検証します。

4

1 に答える 1