2

コードの実行時間が非常に長いため、コードを高速化しようとしています。問題がどこにあるかはすでにわかりました。次の例を検討してください。

x<-c((2+2i),(3+1i),(4+1i),(5+3i),(6+2i),(7+2i))
P<-matrix(c(2,0,0,3),nrow=2)
out<-sum(c(0.5,0.5)%*%mtx.exp(P%*%(matrix(c(x,0,0,x),nrow=2)),5))

複雑な値を持つベクトル x があり、ベクトルには 12^11 のエントリがあり、3 行目の合計を計算したいと考えています。(関数 mtx.exp が必要なのは、これが複素行列ベキだからです (関数はパッケージ Biodem にあります)。%^% 関数は複素数引数をサポートしていないことがわかりました。)

だから私の問題は、私が試してみると

sum(c(0.5,0.5)%*%mtx.exp(P%*%(matrix(c(x,0,0,x),nrow=2)),5))

エラーが表示されます:「pot %*% pot のエラー: 適合しない引数です。」したがって、私の解決策はループを使用することでした：

tmp<-NULL
for (i in 1:length(x)){
  tmp[length(tmp)+1]<-sum(c(0.5,0.5)%*%mtx.exp(P%*%matrix(c(x[i],0,0,x[i]),nrow=2),5))
}

しかし、前述のとおり、これには非常に時間がかかります。コードを高速化する方法はありますか? 私もsapplyを試しましたが、ループと同じくらい時間がかかります。

この関数を約 500 回実行する必要があり、最初の試行で 3 時間以上かかったので、助けていただければ幸いです。これはあまり満足のいくものではありません..

ありがとうございます

4

1 に答える 1

1

ベクトルを事前に割り当てることで、コードを高速化できます。

tmp <- rep(NA,length(x))

しかし、私はあなたが計算しようとしていることを本当に理解していません.最初の例で
は、非正方行列のべき乗を取得しようとしています.2番目の例では、対角行列のべき乗を取得しています(これは実行できます)と^）。

以下はあなたの計算と同等のようです：

sum(P^5/2) * x^5

編集

Pが対角線でもスカラーでもない場合C、 の簡単な単純化は見当たりませんmtx.exp( P %*% C, 5 )。

次のようなものを試すことができます

y <- sapply(x, function(u) 
  sum( 
    c(0.5,0.5) 
    %*% 
    mtx.exp( P %*% matrix(c(u,0,0,u),nrow=2), 5 )
  )
)

しかし、ベクトルに実際に 12^11 のエントリがある場合、非常に長い時間がかかります。

あるいは、非常に小さい (2*2) 行列が非常に多数あるため、P %*% C (Maxima、Sage、Yacas、Maple などのコンピューター代数システムを使用して) 積とその 5 乗を明示的に計算し、結果の式: これらはベクトルに対する単純な操作 (50 行) です。

/* Maxima code */ 
p: matrix([p11,p12], [p21,p22]);
c: matrix([c1,0],[0,c2]);
display2d: false;
factor(p.c . p.c . p.c . p.c . p.c);

次に、結果をコピーして R に貼り付けます。

c1 <- dnorm(abs(x),0,1); # C is still a diagonal matrix
c2 <- dnorm(abs(x),1,3);
p11 <- P[1,1]
p12 <- P[1,2]
p21 <- P[2,1]
p22 <- P[2,2]
# Result of the Maxima computations: 
# I just add all the elements of the resulting 2*2 matrix,
# but you may want to do something slightly different with them.

          c1*(c2^4*p12*p21*p22^3+2*c1*c2^3*p11*p12*p21*p22^2
                                +2*c1*c2^3*p12^2*p21^2*p22
                                +3*c1^2*c2^2*p11^2*p12*p21*p22
                                +3*c1^2*c2^2*p11*p12^2*p21^2
                                +4*c1^3*c2*p11^3*p12*p21+c1^4*p11^5)
          +
          c2*p12
            *(c2^4*p22^4+c1*c2^3*p11*p22^3+3*c1*c2^3*p12*p21*p22^2
                        +c1^2*c2^2*p11^2*p22^2+4*c1^2*c2^2*p11*p12*p21*p22
                        +c1^3*c2*p11^3*p22+c1^2*c2^2*p12^2*p21^2
                        +3*c1^3*c2*p11^2*p12*p21+c1^4*p11^4)
         +
         c1*p21
            *(c2^4*p22^4+c1*c2^3*p11*p22^3+3*c1*c2^3*p12*p21*p22^2
                        +c1^2*c2^2*p11^2*p22^2+4*c1^2*c2^2*p11*p12*p21*p22
                        +c1^3*c2*p11^3*p22+c1^2*c2^2*p12^2*p21^2
                        +3*c1^3*c2*p11^2*p12*p21+c1^4*p11^4)
         +
         c2*(c2^4*p22^5+4*c1*c2^3*p12*p21*p22^3
                        +3*c1^2*c2^2*p11*p12*p21*p22^2
                        +3*c1^2*c2^2*p12^2*p21^2*p22
                        +2*c1^3*c2*p11^2*p12*p21*p22
                        +2*c1^3*c2*p11*p12^2*p21^2+c1^4*p11^3*p12*p21)
于 2012-04-04T13:19:10.930 に答える