为啥这个循环的时间复杂度是非线性的？
Posted
技术标签:
【中文标题】为啥这个循环的时间复杂度是非线性的？【英文标题】：Why is the time complexity of this loop non-linear?为什么这个循环的时间复杂度是非线性的？ 【发布时间】：2016-04-21 17:37:55 【问题描述】：为什么这个循环的时间复杂度是非线性的，为什么这么慢？该循环采用~38s for N=50k, 和~570s for N=200k。有没有更快的方法来做到这一点？ Rprof() 似乎表示写入内存很慢。
df <- data.frame(replicate(5, runif(200000)))
df[,1:3] <- round(df[,1:3])
Rprof(line.profiling = TRUE); timer <- proc.time()
x <- df; N <- nrow(df); i <- 1
ind <- df[1:(N-1),1:3] == df[2:N,1:3];
rind <- which(apply(ind,1,all))
N <- length(rind)
while(i <= N)
x$X4[rind[i]+1] <- x$X4[rind[i]+1] + x$X4[rind[i]]
x$X5[rind[i]+1] <- x$X4[rind[i]+1] * x$X3[rind[i]+1]
x$X5[rind[i]+1] <- trunc(x$X5[rind[i]+1]*10^8)/10^8
x$X1[rind[i]] <- NA
i <- i + 1
;x <- na.omit(x)
proc.time() - timer; Rprof(NULL)
summaryRprof(lines = "show")
此算法的目的是遍历数据框并组合在某些元素上匹配的相邻行。也就是说，它会删除其中一行并将该行的一些值添加到另一行。生成的数据帧应少 n 行，其中 n 是原始数据帧中匹配的相邻行数。每次组合一对行时，源数据帧和新数据帧的索引就会不同步 1，因为从新帧中删除/省略了一行，因此i 会跟踪源数据框，q 跟踪新数据框上的位置。
感谢@joran 的评论，上面的代码得到了更新。性能大幅提升至~5.5s for N=50k 和~88s for N=200k。但是，时间复杂度仍然是非线性的，我无法理解。我需要在 N = 100 万或更多的情况下运行它，所以它的速度仍然不是很快。
【问题讨论】：
好像你在写 c++，r 有相应的包 @rawr heh，你是基于代码还是我的个人资料？我不知道有更多R 这样做的方式。有点违反直觉的是，用于处理数据集的语言会扼杀这么简单的东西。不过，仍然希望我做错了。那些包裹呢？
两者。绝对不是这样做的“r方式”。就像在 c 中做“r 方式”不是最优的。它可能有助于描述输入和期望的结果。或者只是用 c++ 编写它并使用其中一个包来即时编译它
我认为你对设计和使用 R 的人的“典型”数据处理的期望是错误的。这远非典型。也就是说，为什么要遍历每一行？您可以以矢量化方式识别相邻的重复行：ind <- df[1:(N-1),1:3] == df[2:N,1:3]; dup_row_ind <- which(apply(ind,1,all)) + 1。此外，如果您的数据都是数字的，那么使用矩阵而不是数据框可能会更好。
@joran: s/可能会更好/会更好/。使用矩阵而不是 data.frame 将我机器上 OP 示例的运行时间减少了一个数量级，并且它与矩阵中的行数成线性关系（例如 5e4 行需要 ~1s，5e5 行需要 ~10s ）。并将代码放入一个函数并进行字节编译，它似乎进一步缩短了运行时间（在我的机器上减少了一半）。
【参考方案1】：
只有X4 列的更新取决于先前的值，因此循环可以大部分被“向量化”（稍加优化，避免在每次迭代中将 1 添加到rind）
rind1 <- rind + 1L
for (i in seq_len(N))
x$X4[rind1[i]] <- x$X4[rind1[i]] + x$X4[rind[i]]
x$X5[rind1] <- x$X4[rind1] * x$X3[rind1]
x$X5[rind1] <- trunc(x$X5[rind1] * 10^8) / 10^8
x$X1[rind] <- NA
na.omit(x)
X4 是一个数值，可以通过将其更新为向量而不是 data.frame 的列来提高更新效率
X4 <- x$X4
for (i in seq_len(N))
X4[rind1[i]] <- X4[rind1[i]] + X4[rind[i]]
x$X4 <- X4
为了比较，我们有
f0 <- function(nrow)
set.seed(123)
df <- data.frame(replicate(5, runif(nrow)))
df[,1:3] <- round(df[,1:3])
x <- df; N <- nrow(df); i <- 1
ind <- df[1:(N-1),1:3] == df[2:N,1:3];
rind <- which(apply(ind,1,all))
N <- length(rind)
while(i <= N)
x$X4[rind[i]+1] <- x$X4[rind[i]+1] + x$X4[rind[i]]
x$X5[rind[i]+1] <- x$X4[rind[i]+1] * x$X3[rind[i]+1]
x$X5[rind[i]+1] <- trunc(x$X5[rind[i]+1]*10^8)/10^8
x$X1[rind[i]] <- NA
i <- i + 1
na.omit(x)
f1a <- function(nrow)
set.seed(123)
df <- data.frame(replicate(5, runif(nrow)))
df[,1:3] <- round(df[,1:3])
x <- df; N <- nrow(df)
ind <- df[1:(N-1),1:3] == df[2:N,1:3];
rind <- which(apply(ind,1,all))
rind1 <- rind + 1L
for (i in seq_along(rind))
x$X4[rind1[i]] <- x$X4[rind1[i]] + x$X4[rind[i]]
x$X5[rind1] <- x$X4[rind1] * x$X3[rind1]
x$X5[rind1] <- trunc(x$X5[rind1] * 10^8) / 10^8
x$X1[rind] <- NA
na.omit(x)
f4a <- function(nrow)
set.seed(123)
df <- data.frame(replicate(5, runif(nrow)))
df[,1:3] <- round(df[,1:3])
x <- df; N <- nrow(df)
ind <- df[1:(N-1),1:3] == df[2:N,1:3];
rind <- which(apply(ind,1,all))
rind1 <- rind + 1L
X4 <- x$X4
for (i in seq_along(rind))
X4[rind1[i]] <- X4[rind1[i]] + X4[rind[i]]
x$X4 <- X4
x$X1[rind] <- NA
x$X5[rind1] <- X4[rind1] * x$X3[rind1]
x$X5[rind1] <- trunc(x$X5[rind1] * 10^8) / 10^8
na.omit(x)
结果是一样的
> identical(f0(1000), f1a(1000))
[1] TRUE
> identical(f0(1000), f4a(1000))
[1] TRUE
速度提升很大（使用library(microbenchmark)）
> microbenchmark(f0(10000), f1a(10000), f4a(10000), times=10)
Unit: milliseconds
expr min lq mean median uq max neval
f0(10000) 346.35906 354.37637 361.15188 363.71627 366.74944 373.88275 10
f1a(10000) 124.71766 126.43532 127.99166 127.39257 129.51927 133.01573 10
f4a(10000) 41.70401 42.48141 42.90487 43.00584 43.32059 43.83757 10
在启用内存分析的情况下编译 R 时可以看到差异的原因 --
> tracemem(x)
[1] "<0x39d93a8>"
> tracemem(x$X4)
[1] "<0x6586e40>"
> x$X4[1] <- 1
tracemem[0x39d93a8 -> 0x39d9410]:
tracemem[0x6586e40 -> 0x670d870]:
tracemem[0x39d9410 -> 0x39d9478]:
tracemem[0x39d9478 -> 0x39d94e0]: $<-.data.frame $<-
tracemem[0x39d94e0 -> 0x39d9548]: $<-.data.frame $<-
>
每一行表示一个内存副本，因此更新数据帧中的一个单元格会导致外部结构或向量本身的 5 个副本。相比之下，向量可以在没有任何副本的情况下更新。
> tracemem(X4)
[1] "<0xdd44460>"
> X4[1] = 1
tracemem[0xdd44460 -> 0x9d26c10]:
> X4[1] = 2
>
（第一次赋值很昂贵，因为它代表data.frame列的重复；后续更新到X4，只有X4指的是正在更新的向量，向量不需要重复） .
data.frame 实现似乎是非线性扩展的
> microbenchmark(f1a(100), f1a(1000), f1a(10000), f1a(100000), times=10)
Unit: milliseconds
expr min lq mean median uq
f1a(100) 2.372266 2.479458 2.551568 2.524818 2.640244
f1a(1000) 10.831288 11.100009 11.210483 11.194863 11.432533
f1a(10000) 130.011104 138.686445 139.556787 141.138329 141.522686
f1a(1e+05) 4092.439956 4117.818817 4145.809235 4143.634663 4172.282888
max neval
2.727221 10
11.581644 10
147.993499 10
4216.129732 10
原因在上面 tracemem 输出的第二行中很明显——更新一行会触发整个列的副本。因此，该算法的缩放比例为要更新的行数乘以一列中的行数，近似二次。
f4a() 似乎呈线性缩放
> microbenchmark(f4a(100), f4a(1000), f4a(10000), f4a(100000), f4a(1e6), times=10)
Unit: milliseconds
expr min lq mean median uq
f4a(100) 1.741458 1.756095 1.827886 1.773887 1.929943
f4a(1000) 5.286016 5.517491 5.558091 5.569514 5.671840
f4a(10000) 42.906895 43.025385 43.880020 43.928631 44.633684
f4a(1e+05) 467.698285 478.919843 539.696364 552.896109 576.707913
f4a(1e+06) 5385.029968 5521.645185 5614.960871 5573.475270 5794.307470
max neval
2.003700 10
5.764022 10
44.983002 10
644.927832 10
5823.868167 10
人们可以尝试巧妙地对循环进行矢量化，但现在有必要吗？
函数的数据处理部分的优化版本使用负索引（例如，-nrow(df)）从数据框中删除行，rowSums() 代替 apply()，以及 unname() 以便子集操作不'不要携带未使用的名称：
g0 <- function(df)
ind <- df[-nrow(df), 1:3] == df[-1, 1:3]
rind <- unname(which(rowSums(ind) == ncol(ind)))
rind1 <- rind + 1L
X4 <- df$X4
for (i in seq_along(rind))
X4[rind1[i]] <- X4[rind1[i]] + X4[rind[i]]
df$X4 <- X4
df$X1[rind] <- NA
df$X5[rind1] <- trunc(df$X4[rind1] * df$X3[rind1] * 10^8) / 10^8
na.omit(df)
与@Khashaa 建议的 data.table 解决方案相比
g1 <- function(df)
x <- setDT(df)[, r:=rleid(X1, X2, X3),]
x <- x[, .(X1=X1[.N], X2=X2[.N], X3=X3[.N], X4=sum(X4), X5=X5[.N]), by=r]
x <- x[, X5:= trunc(X3 * X4 * 10^8)/10^8]
x
基础 R 版本的性能与时俱进
> n_row <- 200000
> set.seed(123)
> df <- data.frame(replicate(5, runif(n_row)))
> df[,1:3] <- round(df[,1:3])
> system.time(g0res <- g0(df))
user system elapsed
0.247 0.000 0.247
> system.time(g1res <- g1(df))
user system elapsed
0.551 0.000 0.551
（f4a 中的预调版本大约需要 760 毫秒，所以慢了一倍多）。
data.table 实现的结果不正确
> head(g0res)
X1 X2 X3 X4 X5
1 0 1 1 0.4708851 0.8631978
2 1 1 0 0.8977670 0.8311355
3 0 1 0 0.7615472 0.6002179
4 1 1 1 0.6478515 0.5616587
5 1 0 0 0.5329256 0.5805195
6 0 1 1 0.8526255 0.4913130
> head(g1res)
r X1 X2 X3 X4 X5
1: 1 0 1 1 0.4708851 0.4708851
2: 2 1 1 0 0.8977670 0.0000000
3: 3 0 1 0 0.7615472 0.0000000
4: 4 1 1 1 0.6478515 0.6478515
5: 5 1 0 0 0.5329256 0.0000000
6: 6 0 1 1 0.8526255 0.8526255
而我的 data.table 向导（几乎不是 data.table 用户）还不足以知道正确的公式是什么。
编译（仅从 for 循环中受益？）将速度提高约 20%
> g0c <- compiler::cmpfun(g0)
> microbenchmark(g0(df), g0c(df), times=10)
Unit: milliseconds
expr min lq mean median uq max neval
g0(df) 250.0750 262.941 276.1549 276.8848 281.1966 321.3778 10
g0c(df) 214.3132 219.940 228.0784 230.2098 235.4579 242.6636 10
【讨论】：
非常好的解决方案，一如既往。 我在 OP 下的仓促评论并不是一个确切的解决方案:) 正确的版本是x <- setDT(dt)[, r:=rleid(X1, X2, X3),][, s:=.N:1, r];x <- x[x[, .(X4=sum(X4)), by=r], on="r"][, X4 := i.X4,][shift(r)==r, X5:= trunc(X3 * X4 * 10^8) / 10^8,][which(s==1)][,:=(r=NULL, s=NULL, i.X4=NULL),]
@Khashaa R 在:=(r=NULL, ... 抱怨语法错误
:= 应该用反引号关闭。
正确格式为':='(r=NULL, s=NULL, i.X4=NULL),]。【参考方案2】：
以下只是对@Martin Morgan 答案的重写，利用了data.table 的快速子集。它比 data.frame 方法快大约 3 倍。
library(data.table)
library(matrixStats) # for efficient rowAlls function
g01 <- function(df)
setDT(df)
ind <- df[-nrow(df), 1:3, with=FALSE] == df[-1, 1:3, with=FALSE]
rind <- which(rowAlls(ind)) + 1L
X4 <- df$X4
for (i in seq_along(rind))
X4[rind[i]] <- X4[rind[i]] + X4[rind[i] - 1L]
df$X4 <- X4
df$X5[rind] <- trunc(df$X4[rind] * df$X3[rind] * 10^8) / 10^8
df[-rind + 1L,]
g01c <- compiler::cmpfun(g01)
n_row <- 1e6
set.seed(123)
df <- data.frame(replicate(5, runif(n_row)))
df[,1:3] <- round(df[,1:3])
# data.frame
system.time(g0(df))
# user system elapsed
# 1.14 0.00 1.14
system.time(g0c(df))
# user system elapsed
# 0.82 0.03 0.86
# data.table
system.time(g01(df))
# user system elapsed
# 0.40 0.02 0.43
system.time(g01c(df))
# user system elapsed
# 0.12 0.03 0.16
【讨论】：
很酷的方法。data.table 基本上和data.frame 一样只是优化了很多吗？以上是关于为啥这个循环的时间复杂度是非线性的？的主要内容，如果未能解决你的问题，请参考以下文章