向量化函数以避免循环

Question

我正在尝试加快我的代码速度，因为它运行时间很长。我已经发现问题出在哪里了。考虑以下示例：

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 个条目，然后我想计算第三行的总和。 （我需要函数 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 个小时。这不是很令人满意..

非常感谢你

Answer 1

可以通过预先分配你的向量来加速代码，

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

但我不太了解您要计算的内容： 在第一个例子中，
你正试图利用非方阵的力量， 第二，你正在利用对角矩阵的力量 （可以使用^ 完成）。

以下似乎等同于您的计算：

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 等） 并使用得到的公式： 这些只是（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)