R中向量的子向量之和答案

【问题标题】：Sum of subvectors of a vector in RR中向量的子向量之和
【发布时间】：2016-01-18 15:47:46
【问题描述】：

给定一个长度为 k 的向量 x，我想获得一个 k 乘 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）。

【问题讨论】：

试试rollapply()可以吗？
你的for-loops 并没有那么糟糕，因为这些循环中没有任何“增长”。他们只是填充一个已经具有最终大小的矩阵。 for-loops 如果你在循环中构建矩阵in 会很慢，即通过rbind。
你的问题是循环慢还是不够优雅？
主要是优雅的问题。出于实际目的，它已经足够快了。

标签： r

【解决方案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)

也许可以并行运行。

【讨论】：

我也有同样的想法。但是“外部”比“for”循环慢，即使对于更大的向量，例如长度 100。
@mra68 比较慢，因为这个答案只是隐藏了一个循环。 Vectorize 是 mapply 的包装器。如果你有一个矢量化函数，outer 会更快。如果需要一个真正快速的解决方案，将 OP 的代码翻译成 Rcpp 几乎是微不足道的。
感谢您的回答。 foreach 和 Vectorize 是我使用较少的函数之一，我应该多研究一下。

【解决方案2】：

这是另一种方法，它似乎比 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 的编辑功能。

【讨论】：

f4() 函数要快得多。我不得不承认它相当违反直觉且不太透明，但它是一个很好的例子，说明如果以某种方式使用 R 可以多快。谢谢！

【解决方案3】：

您不需要在内部循环中反复重新计算总和，相反，您可以使用一个单元格等于其上方的单元格加上右侧对角线上的单元格这一事实，通过下对角线构建矩阵。这应该将算法的阶数从 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])
}

【讨论】：

谢谢。尤其是 cumsum 的使用是我没有想到的。

【解决方案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

【讨论】：

这比for-loops慢，比outer慢。

【解决方案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

【讨论】：

谢谢，这实际上是我现在使用的。这就是为什么我想知道是否也可以在“纯 R”中完成快速实现。但这是 Rcpp 展示其力量的一个很好的例子。

【解决方案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
}

实际上，根据下面的基准测试，它是最快的base R 解决方案（@Roland 的Rcpp 解决方案仍然是最快的）。相对而言，随着向量大小的增加，它也会变得更快（我只比较了f4（由@docendo 提供），因为它是迄今为止最快的base R 解决方案和@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

该算法利用了仅计算一次cumsum(v) 并利用从那里开始的索引。对于非常大的向量，效率与@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 仍然占据主导地位！！！

【讨论】：