Эффективно вычислить суммы строк трехмерного массива в R

Рассмотрим массив a:

> a <- array(c(1:9, 1:9), c(3,3,2))
> a
, , 1

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

, , 2

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

Как мы можем эффективно вычислить суммы строк матриц, проиндексированных в третьем измерении, так что результат будет следующим:

     [,1] [,2]
[1,]   12   12
[2,]   15   15
[3,]   18   18

??

Суммы столбцов легко через 'dims' аргумент colSums():

> colSums(a, dims = 1)

но я не могу найти способ использовать rowSums() на массиве для достижения желаемого результата, так как он имеет другую интерпретацию 'dims' к тому из colSums(),

Вы можете легко вычислить желаемые суммы строк, используя:

> apply(a, 3, rowSums)
     [,1] [,2]
[1,]   12   12
[2,]   15   15
[3,]   18   18

но это просто скрывает петлю. Существуют ли другие эффективные, по-настоящему векторизованные способы вычисления требуемых сумм строк?

Источник

4 ответа

Решение

Ответ @ Фойтасека о разделении массива напомнил мне о aperm() функция, которая позволяет переставлять размеры массива. Как colSums() работает, мы можем поменять местами первые два измерения, используя aperm() и беги colSums() на выходе.

> colSums(aperm(a, c(2,1,3)))
     [,1] [,2]
[1,]   12   12
[2,]   15   15
[3,]   18   18

Некоторые моменты сравнения этого и других предложенных основанных на R ответов:

> b <- array(c(1:250000, 1:250000),c(5000,5000,2))
> system.time(rs1 <- apply(b, 3, rowSums))
   user  system elapsed 
  1.831   0.394   2.232 
> system.time(rs2 <- rowSums3d(b))
   user  system elapsed 
  1.134   0.183   1.320 
> system.time(rs3 <- sapply(1:dim(b)[3], function(i) rowSums(b[,,i])))
   user  system elapsed 
  1.556   0.073   1.636
> system.time(rs4 <- colSums(aperm(b, c(2,1,3))))
   user  system elapsed 
  0.860   0.103   0.966

Так что в моей системе aperm() Решение выглядит незначительно быстрее:

> sessionInfo()
R version 2.12.1 Patched (2011-02-06 r54249)
Platform: x86_64-unknown-linux-gnu (64-bit)

Тем не мение, rowSums3d() не дает такие же ответы, как другие решения:

> all.equal(rs1, rs2)
[1] "Mean relative difference: 0.01999992"
> all.equal(rs1, rs3)
[1] TRUE
> all.equal(rs1, rs4)
[1] TRUE
Источник

Вы можете разбить массив на два измерения, вычислить суммы строк по нему, а затем собрать результаты обратно так, как вы этого хотите. Вот так:

rowSums3d <- function(a){
    m <- matrix(a,ncol=ncol(a))
    rs <- rowSums(m)
    matrix(rs,ncol=2)
}

> a <- array(c(1:250000, 1:250000),c(5000,5000,2))
> system.time(rowSums3d(a))
   user  system elapsed 
   1.73    0.17    1.96 
> system.time(apply(a, 3, rowSums))
   user  system elapsed 
   3.09    0.46    3.74
Источник

Я не знаю о наиболее эффективном способе сделать это, но sapply кажется, хорошо

a <- array(c(1:9, 1:9), c(3,3,2))
x1 <- sapply(1:dim(a)[3], function(i) rowSums(a[,,i]))
x1
     [,1] [,2]
[1,]   12   12
[2,]   15   15
[3,]   18   18

x2 <- apply(a, 3, rowSums)
all.equal(x1, x2)
[1] TRUE

Что дает улучшение скорости следующим образом:

> a <- array(c(1:250000, 1:250000),c(5000,5000,2))

> summary(replicate(10, system.time(rowSums3d(a))[3]))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.784   2.799   2.810   2.814   2.821   2.862 

> summary(replicate(10, system.time(apply(a, 3, rowSums))[3]))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.730   2.755   2.766   2.776   2.788   2.839 

> summary(replicate(10, system.time( sapply(1:dim(a)[3], function(i) rowSums(a[,,i])) )[3]))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.840   1.852   1.867   1.872   1.893   1.914

Сроки были сделаны на:

# Ubuntu 10.10
# Kernal Linux 2.6.35-27-generic
> sessionInfo()
R version 2.12.1 (2010-12-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Источник

Если у вас многоядерная система, вы можете написать простую C-функцию и использовать библиотеку параллельных потоков Open MP. Я сделал что-то подобное для моей проблемы, и я получил 8-кратное увеличение в 8-ядерной системе. Код будет по-прежнему работать в однопроцессорной системе и даже компилироваться в системе без OpenMP, возможно, с небольшим количеством #ifdef _OPENMP здесь и там.

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

Источник
Другие вопросы по тегам