Vectorize Funktion zur Vermeidung von Schleife

Question

Ich versuche, meinen Code zu beschleunigen, weil es sehr lange läuft. Ich habe bereits herausgefunden, wo das Problem liegt. Man betrachte das folgende Beispiel:

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

Ich habe einen Vektor x mit komplexen Werten, hat der Vektor 12^11 Einträge und dann möchte ich die Summe in der dritten Reihe berechnen. (Ich brauche die Funktion mtx.exp, weil es eine komplexe Matrixleistung ist (die Funktion ist im Paket Biodem). Ich fand heraus, dass die% ^% Funktion komplexe Argumente nicht unterstützt.)

Also mein Problem ist, dass wenn ich versuche

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

ich einen Fehler: "Fehler im Topf% *% Topf. nicht-anpassungs Argumente" Also meine Lösung war eine Schleife zu verwenden:

Aber wie gesagt, das dauert sehr lange. Haben Sie Ideen, wie Sie den Code beschleunigen können? Ich habe auch versucht, sapply, aber das dauert genauso lange wie die Schleife.

Ich hoffe, Sie können mir helfen, weil ich diese Funktion etwa 500 Mal ausführen muss und dies dauerte in ersten Versuch mehr als 3 Stunden. Welche ist nicht sehr befriedigend ..

sehr viel Thank u

Answer 1

Der Code kann durch Pre-Zuteilung Ihre Vektor beschleunigt werden,

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

aber ich verstehe nicht wirklich, was Sie versuchen zu berechnen: in dem ersten Beispiel
Sie die Kraft einer nicht quadratischen Matrix zu nehmen versuchen, in der zweiten, nehmen Sie die Kraft einer Diagonalmatrix (die mit ^ getan werden kann).

Die folgende scheint Ihre Berechnungen äquivalent zu sein:

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

EDIT

Wenn P nicht diagonal und C nicht Skalar ist, ich keine einfache Vereinfachung der mtx.exp(P %*% C, 5) sehen.

Man könnte so etwas wie

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

versuchen, aber wenn Ihr Vektor hat wirklich 12^11 Einträge, , die eine unglaublich lange Zeit in Anspruch nehmen wird.

Alternativ kann, da Sie eine sehr große Anzahl von sehr kleinen haben (2 * 2) Matrizen können Sie explizit das Produkt berechnen P %*% C und seine 5. Leistung (einige Computer-Algebra-System: Maxima, Salbei, Yacas Maple, etc.) und verwenden Sie die resultierenden Formeln: dies sind nur (50 Zeilen) einfache Operationen auf Vektoren.

/* 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); 

ich dann kopieren und das Ergebnis in R einfügen:

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)