rpart: 分类变量和连续变量的回归器计算时间对比

8
我目前正在使用 rpart 包来拟合一个回归树,数据具有相对较少的观测值和数千个分类预测变量,这些变量只有两个可能的值。从在较小的数据上测试该包的结果来看,在这种情况下,无论我将回归变量声明为分类变量(即因子)还是保留它们(它们被编码为+/-1),都没有关系。
然而,我仍然想了解为什么将我的解释变量作为因子传递会显著减慢算法的速度(至少因为我很快就会得到新数据,响应将采取3个不同的值,将其视为连续变量将不再是一个选项)。这难道不应该是反过来吗?
以下是模拟我的数据的示例代码:
library(rpart)

x <- as.data.frame(matrix(sample(c(-1, +1), 50 * 3000, replace = T), nrow = 50))
y <- rnorm(50)

x.fac <- as.data.frame(lapply(x, factor))

现在进行比较:

system.time(rpart( y ~ ., data = x, method = 'anova'))

   user  system elapsed 
   1.62    0.21    1.85 

system.time(rpart( y ~ ., data = x.fac, method = 'anova'))

   user  system elapsed 
   246.87  165.91  412.92 

对于每个变量(因子)只处理一种可能的分裂情况比处理整个潜在分裂范围(对于连续变量)更简单和更快,因此我最困惑的是rpart的行为。任何澄清/建议将不胜感激。

2个回答

6
你需要对代码进行剖析才能确定,但如果计时差异不是因为R将每个因子变量转换为两个二进制变量以准备模型矩阵,我会感到惊讶的。

尝试

Rprof("rpartProfile.Rprof")
rpart( y ~ ., data = x.fac, method = 'anova')
Rprof()

summaryRprof("rpartProfile.Rprof")

并查看时间花费的位置。现在我已经做到了:

> summaryRprof("rpartProfile.Rprof")
$by.self
                          self.time self.pct total.time total.pct
"[[<-.data.frame"            786.46    72.45     786.56     72.46
"rpart.matrix"               294.26    27.11    1081.78     99.66
"model.frame.default"          1.04     0.10       3.00      0.28
"terms.formula"                0.96     0.09       0.96      0.09
"as.list.data.frame"           0.46     0.04       0.46      0.04
"makepredictcall.default"      0.46     0.04       0.46      0.04
"rpart"                        0.44     0.04    1085.38     99.99
"[[.data.frame"                0.16     0.01       0.42      0.04
"<Anonymous>"                  0.16     0.01       0.18      0.02
"match"                        0.14     0.01       0.22      0.02
"print"                        0.12     0.01       0.12      0.01
"model.matrix.default"         0.10     0.01       0.44      0.04
....

$by.total
                          total.time total.pct self.time self.pct
"rpart"                      1085.38     99.99      0.44     0.04
"rpart.matrix"               1081.78     99.66    294.26    27.11
"[[<-"                        786.62     72.47      0.06     0.01
"[[<-.data.frame"             786.56     72.46    786.46    72.45
"model.frame.default"           3.00      0.28      1.04     0.10
"eval"                          3.00      0.28      0.04     0.00
"eval.parent"                   3.00      0.28      0.00     0.00
"model.frame"                   3.00      0.28      0.00     0.00
"terms.formula"                 0.96      0.09      0.96     0.09
"terms"                         0.96      0.09      0.00     0.00
"makepredictcall"               0.50      0.05      0.04     0.00
"as.list.data.frame"            0.46      0.04      0.46     0.04
"makepredictcall.default"       0.46      0.04      0.46     0.04
"as.list"                       0.46      0.04      0.00     0.00
"vapply"                        0.46      0.04      0.00     0.00
"model.matrix.default"          0.44      0.04      0.10     0.01
"[["                            0.44      0.04      0.02     0.00
"model.matrix"                  0.44      0.04      0.00     0.00
....

$sample.interval
[1] 0.02

$sampling.time
[1] 1085.5

请注意,上面提到的大部分时间都花在了函数rpart.matrix中:
> rpart:::rpart.matrix
function (frame) 
{
    if (!inherits(frame, "data.frame") || is.null(attr(frame, 
        "terms"))) 
        return(as.matrix(frame))
    for (i in 1:ncol(frame)) {
        if (is.character(frame[[i]])) 
            frame[[i]] <- as.numeric(factor(frame[[i]]))
        else if (!is.numeric(frame[[i]])) 
            frame[[i]] <- as.numeric(frame[[i]])
    }
    X <- model.matrix(attr(frame, "terms"), frame)[, -1L, drop = FALSE]
    colnames(X) <- sub("^`(.*)`", "\\1", colnames(X))
    class(X) <- c("rpart.matrix", class(X))
    X
}

但是在该函数中,大部分时间花费在 for 循环上,该循环将每列数据进行转换并添加回数据框中。


看起来不是 model.matrix 本身的问题,而是 rpart:::rpart.matrix 需要在 for 循环中扩展模型简写 y~.。这可能也是人们不建议在随机森林中使用公式接口的原因。 - joran
+1 @joran 确实 - 我正在对其进行分析,但在我的电脑上需要一些时间。现在已经添加了这些细节。 - Gavin Simpson
@GavinSimpson 谢谢你,Gavin。是的,Rprof 对我也给出了类似的结果,因此确定了罪魁祸首是 rpart.matrix。 - stas g
@joran 我不确定这是否是一个简单的问题,但如果函数参数需要,是否可以放弃公式语法呢? - stas g
@stasg 我其实有点惊讶,因为在这种情况下似乎没有非公式接口可用于 rpart。我会调查一下下面 Hong 的解决方法。 - joran
那么这个答案确定了问题,但解决方案是什么? - littleO

6

延续@gavin simpson的发现...我决定挑战一下rpart.matrix,看看能否解决执行时间过长的问题。

问题在于使用了for循环。通常情况下,我对for和[sl]apply持中立态度;虽然后者通常被认为更加优雅,但只要for正常工作,我就不会仅仅因为这个原因而替换掉它。特别地,我认为*apply的性能优势有时候被过度吹嘘了;相比S-Plus时代,for在速度和内存使用方面得到了显著改进。

然而,在这种情况下并非如此。仅仅将for替换成lapply就可以使这个示例的运行时间缩短超过两个数量级。期待其他人也能确认这一点。

m <- model.frame(x.fac)


# call rpart.matrix
system.time(mm <- rpart:::rpart.matrix(m))
   user  system elapsed 
 208.25   88.03  296.99 


# exactly the same as rpart.matrix, but with for replaced by lapply
f <- function(frame)
{
    if (!inherits(frame, "data.frame") || is.null(attr(frame, 
        "terms"))) 
        return(as.matrix(frame))
    frame[] <- lapply(frame, function(x) {
        if (is.character(x))
            as.numeric(factor(x))
        else if(!is.numeric(x))
            as.numeric(x)
        else x
    })
    X <- model.matrix(attr(frame, "terms"), frame)[, -1L, drop = FALSE]
    colnames(X) <- sub("^`(.*)`", "\\1", colnames(X))
    class(X) <- c("rpart.matrix", class(X))
    X
}

system.time(mm2 <- f(m))
   user  system elapsed 
   0.65    0.04    0.70 


identical(mm, mm2)
[1] TRUE

谢谢您!我在我的机器上运行了这个程序,得到了类似的结果:原始循环函数为 248.64 170.68 421.89,而重写的 lapply 函数为 0.81 0.20 1.57。我认为,这不仅仅是一个例外,更是证明了 apply 函数族比循环更快,并且利用了 R 语言向量化的特性。但是,我发现调用 model.frame 的时间太长了。而且,由于 rpart.matrix 是一个隐藏函数,我不确定这对我有什么影响! - stas g
@stasg 要小心得出结论,即for循环很慢!这里并不是真的情况。问题在于 for 循环中的代码,而不是循环本身。所有在循环中发生的赋值都是问题所在。lapply的匿名函数可以规避许多导致对象复制的问题,仅此而已。我可能可以构建一个略有不同但更快的 for 循环。 - joran
我敢打赌你会将此作为对软件包维护者的建议改进提交!;) - joran
@stasg 这里的问题不在于循环速度(尽管 lapply() 的循环基础设施运行速度比 for() 快一些 [稍微],但只有当循环内部的代码是微不足道的时才会产生影响,否则该组件就会占据主导地位),而是像 Joran 所提到的那样,循环内部函数调用的计算成本。lapply 调用正在消耗除一个 [[<.data.frame() 之外的所有调用(还剩下一个 R-side 调用 [<-.data.frame())。数据框架在处理上非常缓慢,我认为这个问题很好地突显了这个问题。 - Gavin Simpson
+1 这是一个非常好的修改。在接受 @joran 的挑战之前,请确保代码没有做任何奇怪的事情,已经对 rpart 包中的所有代码示例和测试(如果有)运行过,并且不依赖于 [<-.data.frame 的未记录行为 - 我对 ?[<-.data.frame 的阅读似乎表明,您的用法是正确的,而且更重要的是有文档支持的正确。 - Gavin Simpson

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接