Более быстрый способ вычисления недиагональных средних в больших матрицах

Мне нужно вычислить среднее значение каждого недиагонального элемента в 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)

Ответы

Ответ 1

Вы можете получить значительно быстрее, просто извлекая диагонали непосредственно с помощью линейной адресации: superdiag здесь извлекает i-ю сверхдиагональную из A (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))})
}

Выполнение этого на квадратной матрице размером 1К дает ускорение ~ 800x:

> 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

Это дает результаты в том же порядке, что и оригинал.

Ответ 2

Вы можете использовать следующую функцию:

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