Modo più rapido di calcolare le medie fuori diagonale nelle grandi matrici

Question

Ho bisogno di calcolare la media di ogni elemento fuori diagonale in una matrice n × n. I triangoli inferiore e superiore sono ridondanti. Ecco il codice che sto usando:

A <- replicate(500, rnorm(500)) 
sapply(1:(nrow(A)-1), function(x) mean(A[row(A) == (col(A) - x)])) 

Quale sembra funzionare ma non si adatta bene con matrici più grandi. Quelli che ho non sono enormi, circa 2-5.000^2, ma anche con il 1000^2 si sta prendendo più di quanto mi piacerebbe:

A <- replicate(1000, rnorm(1000)) 
system.time(sapply(1:(nrow(A)-1), function(x) mean(A[row(A) == (col(A) - x)]))) 
> user system elapsed 
> 26.662 4.846 31.494 

C'è un modo più intelligente di fare questo?

modifica Per chiarire, mi piacerebbe la media di ciascuna diagonale indipendentemente, ad es. per:

1 2 3 4 
1 2 3 4 
1 2 3 4 
1 2 3 4 

Vorrei:

mean(c(1,2,3)) 
mean(c(1,2)) 
mean(1) 

Answer 1

È possibile ottenere significativamente più veloce semplicemente estraendo le diagonali direttamente tramite indirizzamento lineare: superdiag qui estrae il esima superdiagonal da A (i = 1 è la diagonale principale)

superdiag <- function(A,i) { 
    n<-nrow(A); 
    len<-n-i+1; 
    r <- 1:len; 
    c <- i:n; 
    indices<-(c-1)*n+r; 
    A[indices] 
} 

superdiagmeans <- function(A) { 
    sapply(2:nrow(A), function(i){mean(superdiag(A,i))}) 
} 

l'esecuzione di questo su una matrice quadrata 1K dà un ~ 800x aumento di velocità:

> A <- replicate(1000, rnorm(1000)) 

> system.time(sapply(1:(nrow(A)-1), function(x) mean(A[row(A) == (col(A) - x)]))) 
    user system elapsed 
26.464 3.345 29.793 

> system.time(superdiagmeans(A)) 
    user system elapsed 
    0.033 0.006 0.039 

questo ti dà risultati nello stesso ordine dell'originale.

Answer 2

È possibile utilizzare la seguente funzione:

diagmean <- function(x){ 
    id <- row(x) - col(x) 
    sol <- tapply(x,id,mean) 
    sol[names(sol)!='0'] 
} 

Se controlliamo questo sul vostro matrice, il guadagno di velocità è sostanziale:

> system.time(diagmean(A)) 
    user system elapsed 
    2.58 0.00 2.58 

> system.time(sapply(1:(nrow(A)-1), function(x) mean(A[row(A) == (col(A) - x)]))) 
    user system elapsed 
    38.93 4.01 42.98 

Si noti che questa funzione calcola sia i triangoli superiore che quello inferiore. È possibile calcolare ad es. Solo il triangolare inferiore usando:

diagmean <- function(A){ 
    id <- row(A) - col(A) 
    id[id>=0] <- NA 
    tapply(A,id,mean) 
} 

Ciò si traduce in un altro guadagno di velocità. Si noti che la soluzione sarà invertito rispetto al vostro:

> A <- matrix(rep(c(1,2,3,4),4),ncol=4) 

> sapply(1:(nrow(A)-1), function(x) mean(A[row(A) == (col(A) - x)])) 
[1] 2.0 1.5 1.0 

> diagmean(A) 
-3 -2 -1 
1.0 1.5 2.0