計算速度を上げるためにRコードを変更するには

Question

ここに私のRコードがあります。計算のスピードを加速できるように助言してください:)

まず、関数myfun()は複素数を生成します。

第2に、myfun()を使用して行列Mの要素を計算します。

myfun<-function(a,b,nq,ul,uk) 
{ 
    m<-seq(1,(nq/2)+1,length=(nq/2)+1); 
    k<-m; 

    D<-matrix(NA,nrow = length(k),ncol = length(k)); 

    for(i in 1:length(k)) # row 
    for(j in 1:length(m)) # column 
    { 
     D[i,j]<-(2/nq)*cos(((j-1)*(i-1)*pi)/(nq*0.5)) 
    } 

    D[,1]<-D[,1]*0.5; 
    D[,ncol(D)]<-D[,ncol(D)]*0.5; 

    # compute the vector v 
    vseq<-seq(2,nq-2,by=2); 
    vr<-2/(1-vseq^2); 
    vr<-c(1,vr,1/(1-nq*nq)); 
    v<-matrix(vr,ncol=1); # v is a N by 1 matrix 

    # compute the vector w, length(w)=nq/2+1 

    h<-function(x,ul,uk) 
    { 
     ((b-a)/2)*(exp((b-a)/2*x+(a+b)/2)+1)^(1i*uk)*cos(((b-a)/2*x+(a+b)/2-a)*ul) 
    } 

    w<-matrix(rep(NA,length(v)),ncol=1); 

    for(i in 1:length(w)) 
     { 
     w[i]<-h((cos((i-1)*pi/nq)),ul,uk)+h((-cos((i-1)*pi/nq)),ul,uk) 
     } 

     res<-t(t(D)%*%v)%*%w; # each element of matrix M 

     return(res) 
     } 

次に、行列Mの各要素を計算します.N番目の列とN番目の行はゼロです。

matrix.M<-matrix(0,ncol = N,nrow = N); 

    for(i in 1:N-1) 
    for(j in 1:N-1) 
    { 
    matrix.M[i,j]<-myfun(a,b,nq,i-1,j-1) 
    } 

私たちは、私はあなたの助けに感謝

a<--173.2; 
    b<-78; 
    alpha<-0.24; 
    Dt<-0.1; 
    M<-1000; 
    N<-150; 
    u<-seq(1,150,by=1)*pi/(b-a); 
    nq<-3000; 

ようにパラメータを設定することができます！

Answer 1

機能を高速化するためのいくつかの提案があります。ダブルループ

使用最終行列積

でmyfunのための隠された宝石crossprodするのではなく、可能

使用などの多くの機能としてouter機能を

1. ベクタライズ：私は三つの「トリック」を使用します-function（a、b、nq、ul、uk）{ m < -seq（1、（nq/2）+1、length =（nq/2）+1）; k < -m;

   ## Use outer to compute the elements of the matrix 
   D <- outer(1:length(k), 1:length(m), function(i, j) {(2/nq)*cos(((j-1)*(i-1)*pi)/(nq*0.5))}) 
   
   D[,1]<-D[,1]*0.5; 
   D[,ncol(D)]<-D[,ncol(D)]*0.5; 
   
               # compute the vector v                    
   vseq<-seq(2,nq-2,by=2); 
   vr<-2/(1-vseq^2); 
   vr<-c(1,vr,1/(1-nq*nq)); 
   v<-matrix(vr,ncol=1); # v is a N by 1 matrix                       
   
   h<-function(x,ul,uk) { 
       ((b-a)/2)*(exp((b-a)/2*x+(a+b)/2)+1)^(1i*uk)*cos(((b-a)/2*x+(a+b)/2-a)*ul) 
   } 
   
   ## Compute the full w vector in one go  
   vect <- seq_along(v)-1 
   w <- h((cos(vect*pi/nq)),ul,uk) + h((-cos(vect*pi/nq)),ul,uk) 
   
   ## Compute the cross products. 
   res <- crossprod(crossprod(D, v), w) 
   
   return(res) 
   

   }

私は、これが本来の機能に比べて時間の約80％を保存するべきだと思います。時間は、Dの最初の計算でした。お役に立てれば。