Помогите ускорить петлю в R

Question

В основном я хочу выполнить усреднение по диагонали в R. Ниже приведен код, адаптированный из пакета simsalabim для осреднения по диагонали. Только это медленно.

Есть предложения по векторизации вместо использования 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)
}

Для данных вы можете использовать следующее

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

Вы можете ускорить вычисления почти в 10 раз с помощью очень специализированного трюка с 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

[Обновление: как предположил Джошуа, его можно ускорить, используя суть кода строки:

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

Ускорение x34,5 по сравнению с исходным кодом !!

]

Answer 2

Я не могу заставить ваш пример дать ощутимые результаты. Я думаю, что в вашей функции есть несколько ошибок.

1. XX используется в sapply, но не определено в функции
2. sapply работает над 1:N, где N=144 в вашем примере, но x.b имеет только 115 столбцов
3. reconSSA просто возвращает x

Независимо, я думаю, вы хотите:

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

---

ОБНОВЛЕНИЕ: я переделал и профилировал функцию. Большую часть времени тратится на mean, поэтому может быть трудно получить это намного быстрее, используя R-код. Тем не менее, вы можете увеличить скорость на 20%, используя sum.

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