在R中计算一个向量的所有子向量之和

Question

假设给定长度为k的向量x，我想获得一个k x k的矩阵X，其中X[i,j]是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）。

- Theodor

---

尝试使用rollapply()函数，怎么样？ - Gopala

你的for循环并不那么糟糕，因为这些循环中没有任何东西“增长”。它们只是填充一个已经具有最终大小的矩阵。如果你在循环中构建一个矩阵，例如通过rbind，那么for循环会变得很慢。 - mra68

你的问题是循环速度慢还是不够优雅？ - David Arenburg

主要是一个优雅的问题。对于实际目的来说，它足够快速。 - Theodor

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)

也许可以并行运行它。

- J.R.

---

我有同样的想法。但是，“outer”比“for”循环慢，即使对于更大的向量，例如长度为100的向量。 - mra68

由于这个答案只是隐藏了一个循环，所以它变慢了。Vectorize是mapply的包装器。如果您有一个向量化的函数，那么outer会更快。如果需要一个真正快速的解决方案，将OP的代码翻译成Rcpp几乎是微不足道的。 - Roland

感谢你的回答。foreach和Vectorize是我使用较少的函数，我应该更深入地研究它们。 - Theodor

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

- Neal Fultz

---

谢谢。特别是cumsum的使用，这是我没有想到的。 - Theodor

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编辑的函数已更新。

- talat

---

f4()函数速度非常快。我必须承认它相当反直觉和不太透明，但这是R语言在特定使用方式下可以达到的快速的一个很好的例子。谢谢！ - Theodor

Answer 4

您可以尝试以下方法：

x <- 1:10

matrix(apply(expand.grid(1:10, 1:10), 1, function(y) sum(x[y[2]:y[1]])), 10, 10)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    1    3    6   10   15   21   28   36   45    55
 [2,]    3    2    5    9   14   20   27   35   44    54
 [3,]    6    5    3    7   12   18   25   33   42    52
 [4,]   10    9    7    4    9   15   22   30   39    49
 [5,]   15   14   12    9    5   11   18   26   35    45
 [6,]   21   20   18   15   11    6   13   21   30    40
 [7,]   28   27   25   22   18   13    7   15   24    34
 [8,]   36   35   33   30   26   21   15    8   17    27
 [9,]   45   44   42   39   35   30   24   17    9    19
[10,]   55   54   52   49   45   40   34   27   19    10

- DatamineR

---

这比使用“for”循环要慢，也比使用“outer”要慢。 - mra68

Answer 5

这是一个几乎与你的代码字面上相同的Rcpp函数：

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

library(inline)
library(Rcpp)

cppFunction(
  'NumericMatrix allSums(NumericVector x) {
        int n = x.length();
        NumericMatrix X(n, n);
        for (int i = 0; i < n; ++i) {
          for (int j = 0; j < n; ++j) {
             for (int k = i; k <= j; ++k) {
               X(i,j) += x(k);
             }
            X(j,i) = X(i,j);
          }
        }
        return X;
    }')

Y <- allSums(x)
all.equal(X, Y)
#[1] TRUE

当然，这个算法可以改进：

cppFunction(
  'NumericMatrix allSums2(NumericVector x) {
        int n = x.length();
        NumericMatrix X(n, n);
        X(0,0) = x(0);
        for (int j = 0; j < n; ++j) {
          if (j > 0) {
            X(0,j) = X(0, j-1) + x(j);
            X(j,0) = X(0,j);
          }
          for (int i = 1; i < n; ++i) {
            X(i,j) = X(i-1,j) - x(i-1); 
            X(j,i) = X(i,j);
            }
          }
        return X;
    }')

Z <- allSums2(x)
all.equal(X, Z)
#[1] TRUE

一些基准测试结果：

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

f4 <- function(x, n) {
  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
}


microbenchmark(f4(x, n), allSums(x), allSums2(x), times = 25)#
#Unit: microseconds
#       expr      min       lq      mean   median       uq      max neval cld
#   f4(x, n)  933.441  938.061 1121.0901  975.633 1045.232 2635.561    25  b 
# allSums(x) 1385.533 1391.693 1466.4784 1395.080 1408.630 2996.803    25   c
#allSums2(x)  127.499  129.038  198.8475  133.965  139.201 1737.844    25 a  

- Roland

---

谢谢，那实际上是我现在正在使用的。这就是为什么我想知道是否也可以在“纯R”中进行快速实现。但这是一个很好的例子，显示了Rcpp的强大之处。 - Theodor

Answer 6

除了已经提供的优秀答案，这里还有一个超级快速的base R解决方案：

subVecSum <- function(v, s) {
    t <- c(0L, cumsum(v))
    n1 <- s+1L
    m <- matrix(0L,s,s)
    for (i in 4L:n1) {
        m[i-2L,1L:(i-3L)] <- t[i-1L]-t[1L:(i-3L)]
        m[i-2L,i-2L] <- v[i-2L]
        m[i-2L,(i-1L):s] <- t[i:n1]-t[i-2L]
    }
    m[1L,] <- t[-1L]; m[s,] <- t[n1]-t[1L:s]
    m
}

实际上，根据以下基准测试，它是最快的基本R解决方案（@Roland的Rcpp解决方案仍然是最快的）。 相对而言，随着向量大小的增加，它变得更快（我只比较了由@docendo提供的最快的基本R解决方案f4以及@Roland的Rcpp实现。您会注意到，我使用了由@Roland定义的修改后的f4函数）。

## We first compile the functions.. no need to compile the Rcpp
## function as it is already done by calling cppFunction
c.f4 <- compiler::cmpfun(f4)
c.subVS1 <- compiler::cmpfun(subVecSum)

n <- 100
x <- 1:n
microbenchmark(c.f4(x,n), c.subVS1(x,n), allSums2(x), times = 1000, unit = "relative")
Unit: relative
          expr       min        lq     mean    median        uq       max neval cld
    c.f4(x, n) 11.355013 11.262663 9.231756 11.545315 12.074004 1.0819186  1000   c
c.subVS1(x, n)  7.795879  7.592643 5.414135  7.624209  8.080471 0.8490876  1000  b 
   allSums2(x)  1.000000  1.000000 1.000000  1.000000  1.000000 1.0000000  1000 a  

n <- 500
x <- 1:n
microbenchmark(c.f4(x,n), c.subVS1(x,n), allSums2(x), times = 500, unit = "relative")
Unit: relative
          expr      min       lq     mean   median       uq       max neval cld
    c.f4(x, n) 6.231426 6.585118 6.442567 6.438163 6.882862 10.124428   500   c
c.subVS1(x, n) 3.548766 3.271089 3.137887 2.881520 3.604536  8.854241   500  b 
   allSums2(x) 1.000000 1.000000 1.000000 1.000000 1.000000  1.000000   500 a  

n <- 1000
x <- 1:n
microbenchmark(c.f4(x,n), c.subVS1(x,n), allSums2(x), times = 100, unit = "relative")
Unit: relative
          expr      min        lq      mean    median        uq      max neval cld
    c.f4(x, n) 7.779537 16.352334 11.489506 15.529351 14.447210 3.639483   100   c
c.subVS1(x, n) 2.637996  2.951763  2.937385  2.726569  2.692099 1.211545   100  b 
   allSums2(x) 1.000000  1.000000  1.000000  1.000000  1.000000 1.000000   100 a  

identical(c.f4(x,n), c.subVS1(x,n), as.integer(allSums2(x)))  ## gives the same results
[1] TRUE

这个算法仅利用一次计算并从那里利用索引。对于非常大的向量，效率与@Roland提供的Rcpp解决方案相当。 观察：

n <- 5000
x <- 1:n
microbenchmark(c.subVS1(x,n), allSums2(x), times = 10, unit = "relative")
Unit: relative
          expr      min       lq     mean   median       uq      max neval cld
c.subVS1(x, n) 1.900718 1.865304 1.854165 1.865396 1.769996 1.837354    10   b
   allSums2(x) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000    10  a 


n <- 10000
x <- 1:n
microbenchmark(c.subVS1(x,n), allSums2(x), times = 10, unit = "relative")
Unit: relative
          expr      min      lq     mean   median       uq     max neval cld
c.subVS1(x, n) 1.503538 1.53851 1.493883 1.526843 1.496783 1.29196    10   b
   allSums2(x) 1.000000 1.00000 1.000000 1.000000 1.000000 1.00000    10  a 

虽然使用 base R 不错，但是 Rcpp 仍然是当之无愧的王者!!!

- Joseph Wood