2

R のコードをいくつか継承しましたが、実行速度が非常に遅くなります。ほとんどの時間は、フォームの関数の評価に費やされます (異なる被積分関数 G を持つ約 15 のそのような関数があります)。

TMin <- 0.5

F <- function (t, d) {
    result <- ifelse(((d > 0) & (t > TMin)),
                     mapply(function(t, d) integrate(G, lower=0, upper=t, t, d)$value, t, d),
                     0)

    return(result)

}

テストのために、次のダミー関数を使用していますが、実際のコードでは、G は exp()、log()、dlnorm()、plnorm() などを含むはるかに複雑です。

G <- function(x, t, d) {
    mean(rnorm(1e5))
    x + t - d
}

F は、最悪のケースで約 200 万回計算されます。関数は次の 3 つの異なる方法で呼び出されます:
t が単一の数値で d が数値ベクトルである、または
t が数値ベクトルで d が単一の数値である、または
t が数値ベクトルであり数値ベクトルである

この関数を高速化する (簡単な) 方法はありますか?

これまでのところ、(ifelse ループを取り除くために) 次の行に沿ってバリエーションを試しました。

F2 <- function (t,d) {
    TempRes <- mapply(function(t, d) integrate(G, lower=0, upper=t, t, d)$value, t, d)
    TempRes[(d <= 0) | (t <= TMin)] <- 0
    result <- TempRes

    return(result)
}

と

F3 <- function (t,d) {
    result <- rep(0, max(length(t),length(d)))
    test <- ((d > 0) & (t > TMin))
    result[test] <- mapply(function(t, d) integrate(G, lower=0, upper=t, t, d)$value, t, d)[test]

    return(result)
}

しかし、ほぼ同じ時間がかかります。

4

2 に答える 2

2

多数の独立した統合を実行しています。これらの統合を別々のコアで同時に実行することにより、速度を上げることができます (マルチコア プロセッサが利用可能な場合)。問題は、R がデフォルトでシングル スレッド方式で計算を実行することです。ただし、マルチスレッドのサポートを可能にする利用可能なパッケージが多数あります。私は最近、関連するパッケージと機能に関するいくつかの追加情報とともに、こことここでいくつかの同様の質問に答えました。

さらに、@Mike Dunlavey が既に述べたように、基準に一致しないtとの値の統合を実行しないようにする必要があります。d(現在、これらの値に対して不要な関数評価を実行しており、後で結果を 0 で上書きしています)。

以下に可能な改善を追加しました。Gクラスター ノードで関数を評価するには、関数を含む別のファイルを作成する必要があることに注意してください。以下のコードでは、このファイルが呼び出されると仮定していますfunctionG.R

スニペット:

library(doParallel)
F4 <- function(t,d) {
  results = vector(mode="numeric",max(length=length(t),length(d))) # Zero vector

  logicalVector <- ((d > 0) & (t > TMin))
  relevantT <- t[logicalVector]
  relevantD <- d[logicalVector] # when d is single element, NA values created

  if(length(relevantT) > 1 | length(relevantD) > 1)
  {
    if(length(d)==1) # d is only one element instead of vector --> replicate it
      relevantD <- rep(d,length(relevantT))
    if(length(t)==1) # t is only one element instead of vector --> replicate it
      relevantT <- rep(t,length(relevantD))

    cl <- makeCluster(detectCores()); 
    registerDoParallel(cl)
    clusterEvalQ(cl,eval(parse("functionG.R")))

    integrationResults <- foreach(i=1:length(relevantT),.combine="c") %dopar%
    {
      integrate(G,lower=0,upper=relevantT[i],relevantT[i],relevantD[i])$value;
    }
    stopCluster(cl)
    results[logicalVector] <- integrationResults
  }
  else if(length(relevantT==1)) # Cluster overhead not needd
  {
    results[logicalVector] = integrate(G,lower=0,upper=relevantT,relevantT,relevantD)$value;
  }

  return(results)
}

私の CPU には、ハイパースレッディングが有効な 6 つの物理コア (x2) が含まれています。結果は次のとおりです。

> t = -5000:20000
> d = -5000:20000
> 
> start = Sys.time()
> testF3 = F3(t,d)
> timeNeededF3 = Sys.time()-start
> 
> start = Sys.time()
> testF4 = F4(t,d)
> timeNeededF4 = Sys.time()-start;

> timeNeededF3
Time difference of 3.452825 mins
> timeNeededF4
Time difference of 29.52558 secs
> identical(testF3,testF4)
[1] TRUE

このコードの実行中、コアは常に使用されているようです。ただし、コアの周りでデータをより効率的に事前に分割し、その後、個別のコアで型の適用関数を利用することで、このコードをさらに最適化できる可能性があります。

さらに最適化が必要な場合は、関数をさらに詳しく調べることもできintegrateます。設定をいじって、より厳密でない数値近似を許可することで、パフォーマンスを向上させることができます。別の方法として、独自の単純なバージョンの適応シンプソン求積法を実装し、個別のステップサイズをいじることができます。ほとんどの場合、このようにパフォーマンスが大幅に向上する可能性があります (近似でより多くのエラーを許容できる/許容できる場合)。

編集: すべてのシナリオで機能するようにコードを更新:dおよび/またはt有効/無効な数値またはベクトル

コメント への返信 @mawir: その通りです。は、 test が に評価される行にifelse(test, yes, no)対応する値を返し、 に評価される行に対応する値を返します。ただし、 のベクトルを作成するには、最初に式を評価する必要があります。このコードは、これを示しています。yesTRUEnotestFALSEyesyeslength(test)

> t = -5000:5
> d = -5000:5
> 
> start = Sys.time()
> testF1 = F(t,d)
> timeNeededF1 = Sys.time()-start
> timeNeededF1
Time difference of 43.31346 secs
> 
> start = Sys.time()
> testF4 = F4(t,d)
> timeNeededF4 = Sys.time()-start
> timeNeededF4
Time difference of 2.284134 secs

tこのシナリオでは、との最後の 5 つの値のみdが関連しています。ただし、F1関数内では、ベクトルを作成するためにon allとvalues が最初にifelse評価されます。これが、関数の実行に非常に時間がかかる理由です。次に、条件を満たす要素を選択するか、そうでない場合は 0 を選択します。関数はこの問題を回避します。mapplydtyesF4

tさらに、 anddがベクトルではないシナリオで高速化が得られるという。ただし、このシナリオでは並列化は使用されません。t通常、 /の一方または両方dがベクトルであるシナリオで最大のスピードアップを得る必要があります。

ANOTHER EDIT、ローランドのコメントに応えて：別の関数ファイルを作成したくない場合は、潜在的に置き換えることができclusterEvalQ(cl,eval(parse("functionG.R")))ますclusterExport(cl,"G")。

于 2015-06-20T21:46:21.540 に答える