我需要计算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)
通过使用线性寻址直接提取对角线，你可以获得明显更快的速度:< code>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
这将为您提供与原始顺序相同的结果。
您可以使用以下功能：
diagmean <- function(x){
id <- row(x) - col(x)
sol <- tapply(x,id,mean)
sol[names(sol)!='0']
}
如果我们在您的矩阵上检查这一点，速度增益是很大的:
> 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
注意，这个函数同时计算上下三角形。您可以使用以下公式仅计算下三角形:
diagmean <- function(A){
id <- row(A) - col(A)
id[id>=0] <- NA
tapply(A,id,mean)
}
这导致另一个速度提升。请注意，与您的解决方案相比，解决方案将相反：
> 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
