在循环中优化 lm() 函数
Optimizing lm() function in a loop
我在循环中使用 R 内置 lm() 函数来估计自定义统计信息：
for(i in 1:10000)
{
x<-rnorm(n)
reg2<-lm(x~data$Y)
Max[i]<-max(abs(rstudent(reg2)))
}
当增加循环计数器（通常我们想测试超过 10^6 或 10^9 的迭代值以解决精度问题）和 Y 的大小时，这真的很慢。
阅读以下内容后 Stack topic，第一个尝试是尝试使用并行回归优化整体（使用 calm()）：
cls = makeCluster(4)
distribsplit(cls, "test")
distribsplit(cls, "x")
for(i in 1:10000)
{
x<-rnorm(n)
reg2 <- calm(cls, "x ~ test$Y, data = test")
Max[i]<-max(abs(reg2$residuals / sd(reg2$residuals)))
}
与原始的、无与伦比的循环相比，这个版本以更慢的版本结束（6 倍）。我的假设是我们要求在每个循环迭代中创建/销毁线程，并且在 R.
中会大大减慢进程
第二次尝试是根据 this Stack topic 使用 lm.fit():
for(i in 1:10000)
{
x<- rnorm(n)
reg2<- .lm.fit(as.matrix(x), data$Y)
Max[i]<-max(abs(reg2$residuals / sd(reg2$residuals)))
}
与初始版本和原始版本相比，它的处理速度要快得多。这样我们现在有：lm.fit() < lm() < calm()，说到整体处理时间。
但是，我们仍在寻找提高此代码效率（在处理时间方面）的选项。有哪些可能的选择？我假设使循环并行会节省一些处理时间？
编辑：最小示例
这是一个最小的例子：
#Import data
sample <- read.csv("sample.txt")
#Preallocation
Max <- vector(mode = "numeric", length = 100)
n <- length(sample$AGE)
x <- matrix(rnorm(100 * n), 100)
for(i in 1 : 100)
{
reg <- lm(x ~ data$AGE)
Max[i] <- max(abs(rstudent(reg)))
}
使用以下数据集 'sample.txt'：
AGE
51
22
46
52
54
43
61
20
66
27
从这里开始，我们进行了几次测试并记录了以下内容：
- 根据@Karo 的贡献，我们在循环外生成了正常样本的矩阵，以节省一些执行时间。我们预计会产生明显的影响，但 运行 测试表明这样做会产生意想不到的相反结果（ 即 更长的执行时间）。也许当增加模拟次数时效果相反。
- 根据@BenBolker 的建议，我们还测试了
fastlm()，它减少了执行时间，但与典型的 lm() 相比，结果似乎有所不同（0.05 倍）
我们仍在努力有效地减少执行时间。按照@Karo 的建议，我们将尝试直接将向量传递给 lm() 并研究并行化（但由于未知原因 calm() 失败）。
因为我还不能评论:
尽量避免在 R 中出现循环。出于某种原因，您每次迭代都会重新计算这些随机数。你可以在没有循环的情况下做到这一点：
duration_loop <- system.time({
for(i in 1:10000000)
{
x <- rnorm(10)
}
})
duration <- system.time({
m <- matrix(rnorm(10000000*10), 10000000)
})
这两种方法都应该为每个 iteration/matrix 行创建 10 个随机值，并且 iterations/rows 的数量相同。尽管这两种方式似乎都是线性扩展的，但您应该会看到执行时间的差异，循环可能会受到 CPU 限制，而“矢量化”方式可能会受到内存限制。
考虑到这一点，您可能应该并且很可能可以完全避免循环，例如，您可以 pass a vector into the lm-function。如果在那之后你仍然需要更快，你肯定可以并行化多种方式，建议如何使用 data.
的工作示例会更容易
上面的评论范围很广，但我会尽量回答一些比较狭隘的问题。
- 我似乎得到了与
.lm.fit 和 fitLmPure 相同的结果（即 all.equal() 是 TRUE），如果我对随机数种子很小心的话:
library(Rcpp)
library(RcppEigen)
library(microbenchmark)
nsim <- 1e3
n <- 1e5
set.seed(101)
dd <- data.frame(Y=rnorm(n))
testfun <- function(fitFn=.lm.fit, seed=NULL) {
if (!is.null(seed)) set.seed(seed)
x <- rnorm(n)
reg2 <- fitFn(as.matrix(x), dd$Y)$residuals
return(max(abs(reg2) / sd(reg2)))
}
## make sure NOT to use seed=101 - also used to pick y -
## if we have y==x then results are unstable (resids approx. 0)
all.equal(testfun(seed=102), testfun(fastLmPure,seed=102)) ## TRUE
fastLmPure 最快（但不是很多）：
(bm1 <- microbenchmark(testfun(),
testfun(lm.fit),
testfun(fastLmPure),
times=1000))
Unit: milliseconds
expr min lq mean median uq max
testfun() 6.603822 8.234967 8.782436 8.332270 8.745622 82.54284
testfun(lm.fit) 7.666047 9.334848 10.201158 9.503538 10.742987 99.15058
testfun(fastLmPure) 5.964700 7.358141 7.818624 7.471030 7.782182 86.47498
如果你想拟合许多独立的响应，而不是许多独立的预测变量（即如果你在回归中改变 Y 而不是 X），你可以提供一个矩阵Y in .lm.fit，而不是循环大量的回归，这可能是一个巨大的胜利。如果您只关心“随机回归的残差”，那可能值得一试。 （不幸的是，提供一个组合的矩阵可能会分离 X 向量 运行s a multiple regression，没有很多 univariate regressions ...）
并行化是值得的，但只会（最多）根据您可用的内核数量进行扩展。只做一个 运行 而不是一组基准测试，因为我很懒惰 ...
运行 5000 次连续复制对我来说大约需要 40 秒（现代 Linux 笔记本电脑）。
system.time(replicate(5000,testfun(fastLmPure), simplify=FALSE))
## user system elapsed
## 38.953 0.072 39.028
运行 在 5 个内核上并行运行大约需要 13 秒，因此 5 个内核的速度提高了 3 倍。如果单个作业更大，这可能会好一点，但显然永远不会比核心数量更好......（8 个核心并没有做得更好）。
library(parallel)
system.time(mclapply(1:5000, function(x) testfun(fastLmPure),
mc.cores=5))
## user system elapsed
## 43.225 0.627 12.970
我认为 higher/coarser 级别的并行化（跨 运行 秒而不是在 lm 适合范围内）会执行得更好。
我想知道您是否可以根据 t 分布的 order statistics 使用分析结果...？
我在循环中使用 R 内置 lm() 函数来估计自定义统计信息：
for(i in 1:10000)
{
x<-rnorm(n)
reg2<-lm(x~data$Y)
Max[i]<-max(abs(rstudent(reg2)))
}
当增加循环计数器（通常我们想测试超过 10^6 或 10^9 的迭代值以解决精度问题）和 Y 的大小时，这真的很慢。
阅读以下内容后 Stack topic，第一个尝试是尝试使用并行回归优化整体（使用 calm()）：
cls = makeCluster(4)
distribsplit(cls, "test")
distribsplit(cls, "x")
for(i in 1:10000)
{
x<-rnorm(n)
reg2 <- calm(cls, "x ~ test$Y, data = test")
Max[i]<-max(abs(reg2$residuals / sd(reg2$residuals)))
}
与原始的、无与伦比的循环相比，这个版本以更慢的版本结束（6 倍）。我的假设是我们要求在每个循环迭代中创建/销毁线程，并且在 R.
中会大大减慢进程第二次尝试是根据 this Stack topic 使用 lm.fit():
for(i in 1:10000)
{
x<- rnorm(n)
reg2<- .lm.fit(as.matrix(x), data$Y)
Max[i]<-max(abs(reg2$residuals / sd(reg2$residuals)))
}
与初始版本和原始版本相比，它的处理速度要快得多。这样我们现在有：lm.fit() < lm() < calm()，说到整体处理时间。
但是，我们仍在寻找提高此代码效率（在处理时间方面）的选项。有哪些可能的选择？我假设使循环并行会节省一些处理时间？
编辑：最小示例
这是一个最小的例子：
#Import data
sample <- read.csv("sample.txt")
#Preallocation
Max <- vector(mode = "numeric", length = 100)
n <- length(sample$AGE)
x <- matrix(rnorm(100 * n), 100)
for(i in 1 : 100)
{
reg <- lm(x ~ data$AGE)
Max[i] <- max(abs(rstudent(reg)))
}
使用以下数据集 'sample.txt'：
AGE
51
22
46
52
54
43
61
20
66
27
从这里开始，我们进行了几次测试并记录了以下内容：
- 根据@Karo 的贡献，我们在循环外生成了正常样本的矩阵，以节省一些执行时间。我们预计会产生明显的影响，但 运行 测试表明这样做会产生意想不到的相反结果（ 即 更长的执行时间）。也许当增加模拟次数时效果相反。
- 根据@BenBolker 的建议，我们还测试了
fastlm()，它减少了执行时间，但与典型的lm() 相比，结果似乎有所不同（0.05 倍）
我们仍在努力有效地减少执行时间。按照@Karo 的建议，我们将尝试直接将向量传递给 lm() 并研究并行化（但由于未知原因 calm() 失败）。
因为我还不能评论:
尽量避免在 R 中出现循环。出于某种原因，您每次迭代都会重新计算这些随机数。你可以在没有循环的情况下做到这一点：
duration_loop <- system.time({
for(i in 1:10000000)
{
x <- rnorm(10)
}
})
duration <- system.time({
m <- matrix(rnorm(10000000*10), 10000000)
})
这两种方法都应该为每个 iteration/matrix 行创建 10 个随机值，并且 iterations/rows 的数量相同。尽管这两种方式似乎都是线性扩展的，但您应该会看到执行时间的差异，循环可能会受到 CPU 限制，而“矢量化”方式可能会受到内存限制。
考虑到这一点，您可能应该并且很可能可以完全避免循环，例如，您可以 pass a vector into the lm-function。如果在那之后你仍然需要更快，你肯定可以并行化多种方式，建议如何使用 data.
上面的评论范围很广，但我会尽量回答一些比较狭隘的问题。
- 我似乎得到了与
.lm.fit和fitLmPure相同的结果（即all.equal()是TRUE），如果我对随机数种子很小心的话:
library(Rcpp)
library(RcppEigen)
library(microbenchmark)
nsim <- 1e3
n <- 1e5
set.seed(101)
dd <- data.frame(Y=rnorm(n))
testfun <- function(fitFn=.lm.fit, seed=NULL) {
if (!is.null(seed)) set.seed(seed)
x <- rnorm(n)
reg2 <- fitFn(as.matrix(x), dd$Y)$residuals
return(max(abs(reg2) / sd(reg2)))
}
## make sure NOT to use seed=101 - also used to pick y -
## if we have y==x then results are unstable (resids approx. 0)
all.equal(testfun(seed=102), testfun(fastLmPure,seed=102)) ## TRUE
fastLmPure 最快（但不是很多）：
(bm1 <- microbenchmark(testfun(),
testfun(lm.fit),
testfun(fastLmPure),
times=1000))
Unit: milliseconds
expr min lq mean median uq max
testfun() 6.603822 8.234967 8.782436 8.332270 8.745622 82.54284
testfun(lm.fit) 7.666047 9.334848 10.201158 9.503538 10.742987 99.15058
testfun(fastLmPure) 5.964700 7.358141 7.818624 7.471030 7.782182 86.47498
如果你想拟合许多独立的响应，而不是许多独立的预测变量（即如果你在回归中改变 Y 而不是 X），你可以提供一个矩阵Y in .lm.fit，而不是循环大量的回归，这可能是一个巨大的胜利。如果您只关心“随机回归的残差”，那可能值得一试。 （不幸的是，提供一个组合的矩阵可能会分离 X 向量 运行s a multiple regression，没有很多 univariate regressions ...）
并行化是值得的，但只会（最多）根据您可用的内核数量进行扩展。只做一个 运行 而不是一组基准测试，因为我很懒惰 ...
运行 5000 次连续复制对我来说大约需要 40 秒（现代 Linux 笔记本电脑）。
system.time(replicate(5000,testfun(fastLmPure), simplify=FALSE))
## user system elapsed
## 38.953 0.072 39.028
运行 在 5 个内核上并行运行大约需要 13 秒，因此 5 个内核的速度提高了 3 倍。如果单个作业更大，这可能会好一点，但显然永远不会比核心数量更好......（8 个核心并没有做得更好）。
library(parallel)
system.time(mclapply(1:5000, function(x) testfun(fastLmPure),
mc.cores=5))
## user system elapsed
## 43.225 0.627 12.970
我认为 higher/coarser 级别的并行化（跨 运行 秒而不是在 lm 适合范围内）会执行得更好。
我想知道您是否可以根据 t 分布的 order statistics 使用分析结果...？