R中向量的子向量之和

Question

给定x长度为k 的向量,我想通过k矩阵获得ak,X其中X[i,j]是k 的总和x[i] + ... + x[j].我现在这样做的方式是

set.seed(1)
x <- rnorm(10)

X <- matrix(0,10,10)
for(i in 1:10) 
  for(j in 1:10)
    X[i,j] <- sum(x[i:j])

#             [,1]       [,2]       [,3]      [,4]        [,5]       [,6]        [,7]      [,8]      [,9]      [,10]
# [1,]  -0.6264538 -0.4428105 -1.2784391 0.3168417  0.64634948 -0.1741189  0.31331014 1.0516348 1.6274162  1.3220278
# [2,]  -0.4428105  0.1836433 -0.6519853 0.9432955  1.27280329  0.4523349  0.93976395 1.6780887 2.2538700  1.9484816
# [3,]  -1.2784391 -0.6519853 -0.8356286 0.7596522  1.08915996  0.2686916  0.75612063 1.4944453 2.0702267  1.7648383
# [4,]   0.3168417  0.9432955  0.7596522 1.5952808  1.92478857  1.1043202  1.59174924 2.3300739 2.9058553  2.6004669
# [5,]   0.6463495  1.2728033  1.0891600 1.9247886  0.32950777 -0.4909606 -0.00353156 0.7347931 1.3105745  1.0051861
# [6,]  -0.1741189  0.4523349  0.2686916 1.1043202 -0.49096061 -0.8204684 -0.33303933 0.4052854 0.9810667  0.6756783
# [7,]   0.3133101  0.9397640  0.7561206 1.5917492 -0.00353156 -0.3330393  0.48742905 1.2257538 1.8015351  1.4961467
# [8,]   1.0516348  1.6780887  1.4944453 2.3300739  0.73479315  0.4052854  1.22575376 0.7383247 1.3141061  1.0087177
# [9,]   1.6274162  2.2538700  2.0702267 2.9058553  1.31057450  0.9810667  1.80153511 1.3141061 0.5757814  0.2703930
# [10,]  1.3220278  1.9484816  1.7648383 2.6004669  1.00518611  0.6756783  1.49614672 1.0087177 0.2703930 -0.3053884

但我不禁感觉必须有更优雅的R方式(除了将其翻译成Rcpp).

Answer 1

我们可以使用outer():

mySum <- function(i,j) sum(x[i:j])
outer(1:10, 1:10, Vectorize(mySum))

编辑:您也可以通过foreach以下方式寻求解决方案:

library(foreach)
mySum <- function(j) sum(x[i:j])
mySum <- Vectorize(mySum)
foreach(i = 1:10, .combine = 'rbind') %do% mySum(1:10)

并且可能并行运行它.

Answer 2

您不需要重复重新计算内循环中的总和,而是可以使用单元格等于其上方的单元格以及右侧对角线上的单元格,通过子对角线构建矩阵.这应该将算法的顺序从O(n ^ 3)减少到O(n ^ 2).

这是一个快速而肮脏的实现:

X <- diag(x)

for(i in 1:9) {
    for(j in 1:(10-i)){
        X[j+i,j] <- X[j,j+i] <- X[j+i,j+i] + X[j+i-1,j]
    }  
}

编辑:

正如其他人所指出的那样,通过使用cumsum和向量化内部循环,您可以获得更快的速度和简单性:

n <- length(x)
X <- diag(x)
for(i in 1:n) {
    X[i:n,i] <- X[i,i:n] <- cumsum(x[i:n])
}

Answer 3

这是另一种方法,它似乎明显快于OP的for循环(因子~30),并且比当前存在的其他答案(因子> = 18)更快:

n <- 5
x <- 1:5
z <- lapply(1:n, function(i) cumsum(x[i:n]))
m <- mapply(function(y, l) c(rep(NA, n-l), y), z, lengths(z))
m[upper.tri(m)] <- t(m)[upper.tri(m)]
m

#     [,1] [,2] [,3] [,4] [,5]
#[1,]    1    3    6   10   15
#[2,]    3    2    5    9   14
#[3,]    6    5    3    7   12
#[4,]   10    9    7    4    9
#[5,]   15   14   12    9    5

基准(向下滚动查看结果)

library(microbenchmark)
n <- 100
x <- 1:n

f1 <- function() {
  X <- matrix(0,n,n)
  for(i in 1:n) {
    for(j in 1:n) {
      X[i,j] <- sum(x[i:j])
    }
  }
  X
}

f2 <- function() {
  mySum <- function(i,j) sum(x[i:j])
  outer(1:n, 1:n, Vectorize(mySum))
}

f3 <- function() {
  matrix(apply(expand.grid(1:n, 1:n), 1, function(y) sum(x[y[2]:y[1]])), n, n)
}

f4 <- function() {
  z <- lapply(1:n, function(i) cumsum(x[i:n]))
  m <- mapply(function(y, l) c(rep(NA, n-l), y), z, lengths(z))
  m[upper.tri(m)] <- t(m)[upper.tri(m)]
  m
}

f5 <- function() {
  X <- diag(x)
  for(i in 1:(n-1)) {
    for(j in 1:(n-i)){
      X[j+i,j] <- X[j,j+i] <- X[j+i,j+i] + X[j+i-1,j]
    }  
  }
  X
}

microbenchmark(f1(), f2(), f3(), f4(), f5(), times = 25L, unit = "relative")
#Unit: relative
# expr      min       lq     mean   median       uq      max neval
# f1() 29.90113 29.01193 30.82411 31.15412 32.51668 35.93552    25
# f2() 29.46394 30.93101 31.79682 31.88397 34.52489 28.74846    25
# f3() 56.05807 53.82641 53.63785 55.36704 55.62439 45.94875    25
# f4()  1.00000  1.00000  1.00000  1.00000  1.00000  1.00000    25
# f5() 16.30136 17.46371 18.86259 17.87850 21.19914 23.68106    25

all.equal(f1(), f2())
#[1] TRUE
all.equal(f1(), f3())
#[1] TRUE
all.equal(f1(), f4())
#[1] TRUE
all.equal(f1(), f5())
#[1] TRUE

更新了Neal Fultz编辑的功能.