2

次の値を計算しようとしています。

1/N * sum[i=0 to N-1]( log(abs(r_i - 2 * r_i * x_i)) )

ここで x_i は次のように再帰的に計算されます:

x_{i+1} = r_i * x_i * (1 - x_i)

すべてのr_is が指定されている場合 (ただし、 で変わりますi)、 andx_0が指定されています。（私が知る限り、この計算を非反復式に単純化してそのように高速化するためのトリッキーな数学的な方法はありません）。

私の問題は、それが非常に遅いことです.外部の視点が私をスピードアップするのに役立つのではないかと思います.

# x0: a scalar. rs: a numeric vector, length N
# N: typically ~5000
f <- function (x0, rs, N) {
    lambda <- 0                                                                 
    x <- x0                                                                     
    for (i in 1:N) {                                                            
        r <- rs[i]                                                              
        rx <- r * x                                                             
        lambda <- lambda + log(abs(r - 2 * rx))                                 
        # calculate the next x value
        x <- rx - rx * x                                                        
    }                                                                           
    return(lambda / N)                                                          
}

現在、この関数自体はかなり高速ですがrs、それぞれ異なるベクトルを使用して、〜 4,000,000 回 (2000 x 2000 マトリックスの各セルに対して 1 回) 呼び出したいと考えています。

しかし、たった 2500 回 (N=1000 の場合) 呼び出しても、次のプロファイルでは約 25 秒かかります。

      self.time self.pct total.time total.pct
"f"       19.98    81.22      24.60    100.00
"*"        2.00     8.13       2.00      8.13
"-"        1.32     5.37       1.32      5.37
"+"        0.70     2.85       0.70      2.85
"abs"      0.56     2.28       0.56      2.28
":"        0.04     0.16       0.04      0.16

これをスピードアップする方法を知っている人はいますか？乗算には時間がかかるように見えますが、繰り返される乗算は既に事前にキャッシュされています。

また、 andの呼び出しを減らすのsum( log(stuff(i)) )と同じ利点を利用しようとし ましたが、これは長さ(数千単位) のベクトルと少なくとも 1 の典型的な値であるため実行不可能であることが判明したため、最終的に Rになりました。log(prod(stuff(i))logabsstuffNprod(stuff)Inf

4

1 に答える 1

5

私の意見では、ボトルネックはfor関数のループです。

Rcpp で次のように書き換えます。

# x0: a scalar. rs: a numeric vector, length N
# N: typically ~5000
x0 <- runif(1)
N <- 5000
rs <- rnorm(5000)
f <- function (x0, rs, N) {
    lambda <- 0                                                                 
    x <- x0                                                                     
    for (i in 1:N) {                                                            
        r <- rs[i]                                                              
        rx <- r * x                                                             
        lambda <- lambda + log(abs(r - 2 * rx))                                 
        # calculate the next x value
        x <- rx - rx * x                                                        
    }                                                                           
    return(lambda / N)                                                          
}

library(inline)
library(Rcpp)
f1 <- cxxfunction(sig=c(Rx0="numeric", Rrs="numeric"), plugin="Rcpp", body='
  double x0 = as<double>(Rx0);
  NumericVector rs(Rrs);
  int N = rs.size();
  double lambda = 0, x = x0, r, rx;
  for(int i = 0;i < N;i++) {
    r = rs[i];
    rx = r * x;
    lambda = lambda + log( fabs(r - 2 * rx) );
    x = rx - rx * x;
  }
  lambda /= N;
  return wrap(lambda);
  ')
f(x0, rs, N)
f1(x0, rs)

library(rbenchmark)

benchmark(f(x0, rs, N), f1(x0, rs))

f1f前回のテストよりも 140 倍高速です。

于 2013-01-24T07:25:39.800 に答える