r - 在大矩阵中计算非对角线平均值的更快方法

Question

我需要计算 n × n 矩阵中每个非对角元素的平均值。下三角和上三角是多余的。这是我目前使用的代码:

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

这似乎有效，但不适用于较大的矩阵。我所拥有的并不大，大约 2-5000^2，但即使有 1000^2，它也需要比我想要的更长的时间:

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  

有没有更聪明的方法来做到这一点？

编辑 为了澄清，我想独立地计算每个对角线的平均值，例如为了:

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

我想:

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

Answer 1

只需使用线性寻址直接提取对角线，您就可以显着加快速度:superdiag这里从 A 中提取第 i 个超对角线(i=1 是主对角线)

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

在 1K 方阵上运行它可以得到约 800 倍的加速:

> 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 

这将为您提供与原始顺序相同的结果。