Ayuda a acelerar un bucle en R

Question

básicamente quiero realizar un promedio diagonal en R. A continuación se muestra un código adaptado del paquete simsalabim para hacer el promedio diagonal. Solo que esto es lento.

¿Alguna sugerencia para vectorizar esto en lugar de usar sapply?

reconSSA <- function(S,v,group=1){ 
### S : matrix 
### v : vector 

    N <- length(v) 
    L <- nrow(S) 
    K <- N-L+1 
    XX <- matrix(0,nrow=L,ncol=K) 
    IND <- row(XX)+col(XX)-1 
    XX <- matrix(v[row(XX)+col(XX)-1],nrow=L,ncol=K) 
    XX <- S[,group] %*% t(t(XX) %*% S[,group]) 

    ##Diagonal Averaging 
    .intFun <- function(i,x,ind) mean(x[ind==i]) 

    RC <- sapply(1:N,.intFun,x=XX,ind=IND) 
    return(RC) 
} 

Para los datos que usted podría utilizar el siguiente

data(AirPassengers) 
v <- AirPassengers 
L <- 30 
T <- length(v) 
K <- T-L+1 

x.b <- matrix(nrow=L,ncol=K) 
x.b <- matrix(v[row(x.b)+col(x.b)-1],nrow=L,ncol=K) 
S <- eigen(x.b %*% t(x.b))[["vectors"]] 
out <- reconSSA(S, v, 1:10) 

Answer 1

Puede acelerar el cálculo en casi 10 veces con la ayuda de un truco muy especializado con rowsum:

reconSSA_1 <- function(S,v,group=1){ 
### S : matrix 
### v : vector 
    N <- length(v) 
    L <- nrow(S) 
    K <- N-L+1 
    XX <- matrix(0,nrow=L,ncol=K) 
    IND <- row(XX)+col(XX)-1 
    XX <- matrix(v[row(XX)+col(XX)-1],nrow=L,ncol=K) 
    XX <- S[,group] %*% t(t(XX) %*% S[,group]) 
    ##Diagonal Averaging 
    SUMS <- rowsum.default(c(XX), c(IND)) 
    counts <- if(L <= K) c(1:L, rep(L, K-L-1), L:1) 
    else c(1:K, rep(K, L-K-1), K:1) 
    c(SUMS/counts) 
} 

all.equal(reconSSA(S, v, 1:10), reconSSA_1(S, v, 1:10)) 
[1] TRUE 

library(rbenchmark) 

benchmark(SSA = reconSSA(S, v, 1:10), 
      SSA_1 = reconSSA_1(S, v, 1:10), 
      columns = c("test", "elapsed", "relative"), 
      order = "relative") 


    test elapsed relative 
2 SSA_1 0.23 1.0000 
1 SSA 2.08 9.0435 

[Actualización: Como sugiere Joshua podría ser acelerar aún más mediante el uso de el quid del código rowsum:

reconSSA_2 <- function(S,v,group=1){ 
### S : matrix 
### v : vector 
    N <- length(v) 
    L <- nrow(S) 
    K <- N-L+1 
    XX <- matrix(0,nrow=L,ncol=K) 
    IND <- c(row(XX)+col(XX)-1L) 
    XX <- matrix(v[row(XX)+col(XX)-1],nrow=L,ncol=K) 
    XX <- c(S[,group] %*% t(t(XX) %*% S[,group])) 
    ##Diagonal Averaging 
    SUMS <- .Call("Rrowsum_matrix", XX, 1L, IND, 1:N, 
        TRUE, PACKAGE = "base") 
    counts <- if(L <= K) c(1:L, rep(L, K-L-1), L:1) 
    else c(1:K, rep(K, L-K-1), K:1) 
    c(SUMS/counts) 
} 

    test elapsed relative 
3 SSA_2 0.156 1.000000 
2 SSA_1 0.559 3.583333 
1 SSA 5.389 34.544872 

Un aumento de velocidad de x34.5 comparando con el código original !!

]

Answer 2

no puede conseguir su ejemplo para producir resultados razonables. Creo que hay varios errores en tu función.

1. XX se utiliza en sapply, pero no está definido en la función
2. sapply trabajos sobre 1:N, donde N=144 en su ejemplo, pero x.b sólo tiene 115 columnas
3. reconSSA simplemente devuelve x

De todos modos, creo que quiere:

data(AirPassengers) 
x <- AirPassengers 
rowMeans(embed(x,30)) 

---

ACTUALIZACIÓN: He vuelto a trabajar y perfilé la función. La mayor parte del tiempo se gasta en mean, por lo que puede ser difícil obtenerlo mucho más rápido utilizando el código R. Dicho esto, puedes acelerar el 20% usando sum en su lugar.

reconSSA <- function(S,v,group=1){ 

    N <- length(v) 
    L <- nrow(S) 
    K <- N-L+1 
    XX <- matrix(0,nrow=L,ncol=K) 
    IND <- row(XX)+col(XX)-1 
    XX <- matrix(v[row(XX)+col(XX)-1],nrow=L,ncol=K) 
    XX <- S[,group] %*% t(t(XX) %*% S[,group]) 

    ##Diagonal Averaging 
    .intFun <- function(i,x,ind) { 
     I <- ind==i 
     sum(x[I])/sum(I) 
    } 

    RC <- sapply(1:N,.intFun,x=XX,ind=IND) 
    return(RC) 
}