Векторизация функции, чтобы избежать цикла

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)) 

У меня есть вектор х с комплексными значениями, вектор имеет 12^11 записей, а затем я хочу, чтобы вычислить сумму в третьем ряду. (Мне нужна функция mtx.exp, потому что это сложная матричная мощность (функция находится в пакете Biodem). Я выяснил, что функция% ^% не поддерживает сложные аргументы.)

Так что моя проблема в том, что если Я стараюсь

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

Я получаю сообщение об ошибке: «ошибка в горшочке% *% банке: без облегать аргументы.» Поэтому моим решением было использовать петлю:

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)) 
} 

Но, как сказано, это занимает очень много времени. У вас есть идеи по ускорению кода? Я также пробовал использовать, но это занимает столько же времени, сколько и цикл.

Надеюсь, вы можете мне помочь, потому что я должен выполнять эту функцию примерно 500 раз, и это потребовало в первую очередь более 3 часов. Но это не очень приятно ..

Спасибо и очень

Answer 1

Код может быть ускорен путем предварительного выделения вашего вектора,

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

, но я не очень понимаю, что вы пытаетесь вычислить: в первом примере,
вы пытаетесь взять власть не-квадратной матрицы, во втором, вы принимаете силу диагональной матрицы (что может быть сделано с ^).

Следующая, кажется, что эквивалентно ваших вычислений:

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

EDIT

Если 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 и его 5-ой степени (с помощью системы компьютерной алгебры: Maxima, шалфей, 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)