如何使用 3 个具有复杂索引进度的 for 循环来加速此计算?
How to speed up this calculation using 3 for loops with a complex progress of indices?
给定以下数据框 df
和一个包含一个值的数值向量 p
:
df <- data.frame(id = c(rep(1, 110), rep(2, 290)),
m = c(seq(1, 110), seq(1:290)),
m1 = c(rep(108, 110), rep(288, 290)),
m2 = c(rep(3, 400)),
f1 = c(rep(-100, 110), rep(-50, 290)),
f2 = c(rep(22, 110), rep(15, 290)),
f3 = c(rep(5, 110), rep(0, 290)),
u = c(c(0.12, 0.16, 0.10), rep(0, 107), c(0.085, 0.09, 0.11), rep(0, 287)),
v = c(rep(0.175, 3), rep(0, 107), rep(0.115, 3), rep(0, 287)),
y = rep(0, 400))
df$s <- sqrt(df$m/(df$m1 + df$m2 - 1))/40
p <- 0.01
这是一个片段:
> head(df)
id m m1 m2 f1 f2 f3 u v y s
1 1 1 108 3 -100 22 5 0.12 0.175 0 0.002383656
2 1 2 108 3 -100 22 5 0.16 0.175 0 0.003370999
3 1 3 108 3 -100 22 5 0.10 0.175 0 0.004128614
4 1 4 108 3 -100 22 5 0.00 0.000 0 0.004767313
5 1 5 108 3 -100 22 5 0.00 0.000 0 0.005330018
6 1 6 108 3 -100 22 5 0.00 0.000 0 0.005838742
以下是有关数据的一些事实:
- 变量
id
和m
唯一标识每一行(主键)。
- 变量
m
表示'month'。因此,数据集是一个时间序列。
- 变量
f1
、f2
、f3
、m1
和 m2
对于每个值都是常量 id
个。这些不依赖于变量 m
.
- 对于
id
的每个值,变量 s
、u
和 v
不是常量 ,因此确实取决于m
.
id
的每个值的行数等于m1 + m2 - 1。或等效:id
的每个值的m
的最大值等于m1 + m2 - 1.
目标是使用以下公式计算 y
:
![y_{k}=(f_{1}+f_{2}+f_{3})\cdot \sum_{i=k}^{min(m_{1}+k-1, m_{1}+m_{2}-1)}(\frac{1+p}{1+s_{i}})^{i/12} \cdot \sum_{j=min(k,m_{2})}^{max(1,i-m_{1}+1)}u_{j}\cdot v_{i-j+1}](https://latex.codecogs.com/gif.latex?y_%7Bk%7D=(f_%7B1%7D+f_%7B2%7D+f_%7B3%7D)%5Ccdot&space;%5Csum_%7Bi=k%7D%5E%7Bmin(m_%7B1%7D+k-1,&space;m_%7B1%7D+m_%7B2%7D-1)%7D(%5Cfrac%7B1+p%7D%7B1+s_%7Bi%7D%7D)%5E%7Bi/12%7D&space;%5Ccdot&space;%5Csum_%7Bj=min(k,m_%7B2%7D)%7D%5E%7Bmax(1,i-m_%7B1%7D+1)%7Du_%7Bj%7D%5Ccdot&space;v_%7Bi-j+1%7D)
![\forall k \in \left \{1,...,m_{1}+m_{2}-1 \right \}](https://latex.codecogs.com/gif.latex?%5Cforall&space;k&space;%5Cin&space;%5Cleft&space;%5C%7B1,...,m_%7B1%7D+m_%7B2%7D-1&space;%5Cright&space;%5C%7D)
我已经创建了一个解决方案来做到这一点:
counter <- 0
start <- proc.time()
for(n in 1:nrow(df)){
#index k holds the current value for m
k <- df$m[n]
counter <- counter + 1
#read the current value for m1 and m2
m1 <- df$m1[n]
m2 <- df$m2[n]
counter <- counter + 2
#calculate the sum of f1, f2 and f3.
sum_of_fs <- df$f1[n] + df$f2[n] + df$f3[n]
counter <- counter + 1
#initialize y. Set it to zero.
y <- 0
counter <- counter + 1
for(i in k:min(m1 + k - 1, m1 + m2 - 1)){
#Initialize the sumproduct of u and v. Set it to zero.
sumprod_uv <- 0
counter <- counter + 1
for(j in min(k, m2):max(1, i - m1 + 1)){
sumprod_uv <- sumprod_uv + df$u[j] + df$v[i - j + 1]
counter <- counter + 1
}
z <- ((1 + p)/(1 + df$s[i]))^(i / 12)
y <- y + sumprod_uv * z
counter <- counter + 2
}
y <- y * sum_of_fs
df$y[n] <- y
counter <- counter + 2
}
counter
proc.time() - start
在这段代码中,我添加了 2 个额外的东西:
- 一个名为
counter
的计数器,用于计算执行的语句数。
- 测量脚本持续时间的计时器。
现在的问题是脚本花费的时间太长 运行。对于这个玩具示例,它花费了大约 2 秒(注释掉了计数器语句),这是可以接受的:
user system elapsed
1.829 0.002 1.872
此持续时间对应的语句数为 290,188(脚本完成时 counter
的值 运行ning)
在现实生活中,我有一个包含超过 90k 条记录的数据集。除此之外,真实的数据集稍微复杂一些(组成 id 的是 7 个变量,而不是一个)。我 运行 使用该数据集的脚本,它 运行 大约 17 分钟。
问题是:我怎样才能加速这个算法?应该有一种更简洁的方法来做到这一点。
最简单的改进应该是在循环之前将列重新定义为向量:(+ 在第一个循环中计算 v1
并删除 sum_of_fs
计算,因为它没有在任何地方使用)
# redefine df columns as vectors
dfm <- df$m
dfm1 <- df$m1
dfm2 <- df$m2
u <- df$u
v <- df$v
s <- df$s
start <- proc.time()
for (n in 1:nrow(df)) {
k <- dfm[n]
m1 <- dfm1[n]
m2 <- dfm2[n]
v1 <- min(k, m2)
# sum_of_fs <- df$f1[n] + df$f2[n] + df$f3[n] # not used anywhere !!
y <- 0
for (i in k:min(m1 + k - 1, m1 + m2 - 1)) {
sumprod_uv <- 0
for (j in v1:max(1, i - m1 + 1)) {
sumprod_uv <- sumprod_uv + u[j] + v[i - j + 1]
}
z <- ((1 + p)/(1 + s[i]))^(i / 12)
y <- y + sumprod_uv * z
}
df$y[n] <- y
}
proc.time() - start
对我来说,这会在 0.39
秒内运行(初始方法为 1.03
秒)。
我建议为速度测试创建更复杂的数据集。
她有一个 C++ 变体,它可能比 R.
更快
library(Rcpp)
sourceCpp(code = "#include <Rcpp.h>
#include <vector>
#include <algorithm>
using namespace Rcpp;
// [[Rcpp::export]]
std::vector<double> fun(double &p
, std::vector<int> &dfm
, std::vector<int> &dfm1
, std::vector<int> &dfm2
, std::vector<double> &u
, std::vector<double> &v
, std::vector<double> &s
) {
std::vector<double> yy(s.size());
for(size_t n=0; n<s.size(); ++n) {
int k = dfm[n];
int m1 = dfm1[n];
int m2 = dfm2[n];
int v1 = std::min(k, m2);
double y = 0.;
int ii = std::min(m1 + k - 1, m1 + m2 - 1);
for(int i=std::min(k,ii); i<=std::max(k,ii); ++i) {
double sumprod_uv = 0.;
int jj = std::max(1, i - m1 + 1);
for (int j=std::min(v1, jj); j<=std::max(v1, jj); ++j) {
sumprod_uv += u[j-1] + v[i - j];
}
y += sumprod_uv * std::pow(((1. + p)/(1. + s[i-1])), (i / 12.));
}
yy[n] = y;
}
return yy;
}")
system.time(df$y <- fun(p, df$m, df$m1, df$m2, df$u, df$v, df$s))
# user system elapsed
# 0.005 0.000 0.004
包含 f1、f2 和 f3 的问题更新后:
df$y <- fun(p, df$m, df$m1, df$m2, df$u, df$v, df$s) * (df$f1 + df$f2 + df$f3)
为了比较时间我电脑上的时间:
#Your code
# user system elapsed
# 0.358 0.004 0.362
#@minem
# user system elapsed
# 0.090 0.003 0.093
给定以下数据框 df
和一个包含一个值的数值向量 p
:
df <- data.frame(id = c(rep(1, 110), rep(2, 290)),
m = c(seq(1, 110), seq(1:290)),
m1 = c(rep(108, 110), rep(288, 290)),
m2 = c(rep(3, 400)),
f1 = c(rep(-100, 110), rep(-50, 290)),
f2 = c(rep(22, 110), rep(15, 290)),
f3 = c(rep(5, 110), rep(0, 290)),
u = c(c(0.12, 0.16, 0.10), rep(0, 107), c(0.085, 0.09, 0.11), rep(0, 287)),
v = c(rep(0.175, 3), rep(0, 107), rep(0.115, 3), rep(0, 287)),
y = rep(0, 400))
df$s <- sqrt(df$m/(df$m1 + df$m2 - 1))/40
p <- 0.01
这是一个片段:
> head(df)
id m m1 m2 f1 f2 f3 u v y s
1 1 1 108 3 -100 22 5 0.12 0.175 0 0.002383656
2 1 2 108 3 -100 22 5 0.16 0.175 0 0.003370999
3 1 3 108 3 -100 22 5 0.10 0.175 0 0.004128614
4 1 4 108 3 -100 22 5 0.00 0.000 0 0.004767313
5 1 5 108 3 -100 22 5 0.00 0.000 0 0.005330018
6 1 6 108 3 -100 22 5 0.00 0.000 0 0.005838742
以下是有关数据的一些事实:
- 变量
id
和m
唯一标识每一行(主键)。 - 变量
m
表示'month'。因此,数据集是一个时间序列。 - 变量
f1
、f2
、f3
、m1
和m2
对于每个值都是常量id
个。这些不依赖于变量m
. - 对于
id
的每个值,变量s
、u
和v
不是常量 ,因此确实取决于m
. id
的每个值的行数等于m1 + m2 - 1。或等效:id
的每个值的m
的最大值等于m1 + m2 - 1.
目标是使用以下公式计算 y
:
我已经创建了一个解决方案来做到这一点:
counter <- 0
start <- proc.time()
for(n in 1:nrow(df)){
#index k holds the current value for m
k <- df$m[n]
counter <- counter + 1
#read the current value for m1 and m2
m1 <- df$m1[n]
m2 <- df$m2[n]
counter <- counter + 2
#calculate the sum of f1, f2 and f3.
sum_of_fs <- df$f1[n] + df$f2[n] + df$f3[n]
counter <- counter + 1
#initialize y. Set it to zero.
y <- 0
counter <- counter + 1
for(i in k:min(m1 + k - 1, m1 + m2 - 1)){
#Initialize the sumproduct of u and v. Set it to zero.
sumprod_uv <- 0
counter <- counter + 1
for(j in min(k, m2):max(1, i - m1 + 1)){
sumprod_uv <- sumprod_uv + df$u[j] + df$v[i - j + 1]
counter <- counter + 1
}
z <- ((1 + p)/(1 + df$s[i]))^(i / 12)
y <- y + sumprod_uv * z
counter <- counter + 2
}
y <- y * sum_of_fs
df$y[n] <- y
counter <- counter + 2
}
counter
proc.time() - start
在这段代码中,我添加了 2 个额外的东西:
- 一个名为
counter
的计数器,用于计算执行的语句数。 - 测量脚本持续时间的计时器。
现在的问题是脚本花费的时间太长 运行。对于这个玩具示例,它花费了大约 2 秒(注释掉了计数器语句),这是可以接受的:
user system elapsed
1.829 0.002 1.872
此持续时间对应的语句数为 290,188(脚本完成时 counter
的值 运行ning)
在现实生活中,我有一个包含超过 90k 条记录的数据集。除此之外,真实的数据集稍微复杂一些(组成 id 的是 7 个变量,而不是一个)。我 运行 使用该数据集的脚本,它 运行 大约 17 分钟。
问题是:我怎样才能加速这个算法?应该有一种更简洁的方法来做到这一点。
最简单的改进应该是在循环之前将列重新定义为向量:(+ 在第一个循环中计算 v1
并删除 sum_of_fs
计算,因为它没有在任何地方使用)
# redefine df columns as vectors
dfm <- df$m
dfm1 <- df$m1
dfm2 <- df$m2
u <- df$u
v <- df$v
s <- df$s
start <- proc.time()
for (n in 1:nrow(df)) {
k <- dfm[n]
m1 <- dfm1[n]
m2 <- dfm2[n]
v1 <- min(k, m2)
# sum_of_fs <- df$f1[n] + df$f2[n] + df$f3[n] # not used anywhere !!
y <- 0
for (i in k:min(m1 + k - 1, m1 + m2 - 1)) {
sumprod_uv <- 0
for (j in v1:max(1, i - m1 + 1)) {
sumprod_uv <- sumprod_uv + u[j] + v[i - j + 1]
}
z <- ((1 + p)/(1 + s[i]))^(i / 12)
y <- y + sumprod_uv * z
}
df$y[n] <- y
}
proc.time() - start
对我来说,这会在 0.39
秒内运行(初始方法为 1.03
秒)。
我建议为速度测试创建更复杂的数据集。
她有一个 C++ 变体,它可能比 R.
更快library(Rcpp)
sourceCpp(code = "#include <Rcpp.h>
#include <vector>
#include <algorithm>
using namespace Rcpp;
// [[Rcpp::export]]
std::vector<double> fun(double &p
, std::vector<int> &dfm
, std::vector<int> &dfm1
, std::vector<int> &dfm2
, std::vector<double> &u
, std::vector<double> &v
, std::vector<double> &s
) {
std::vector<double> yy(s.size());
for(size_t n=0; n<s.size(); ++n) {
int k = dfm[n];
int m1 = dfm1[n];
int m2 = dfm2[n];
int v1 = std::min(k, m2);
double y = 0.;
int ii = std::min(m1 + k - 1, m1 + m2 - 1);
for(int i=std::min(k,ii); i<=std::max(k,ii); ++i) {
double sumprod_uv = 0.;
int jj = std::max(1, i - m1 + 1);
for (int j=std::min(v1, jj); j<=std::max(v1, jj); ++j) {
sumprod_uv += u[j-1] + v[i - j];
}
y += sumprod_uv * std::pow(((1. + p)/(1. + s[i-1])), (i / 12.));
}
yy[n] = y;
}
return yy;
}")
system.time(df$y <- fun(p, df$m, df$m1, df$m2, df$u, df$v, df$s))
# user system elapsed
# 0.005 0.000 0.004
包含 f1、f2 和 f3 的问题更新后:
df$y <- fun(p, df$m, df$m1, df$m2, df$u, df$v, df$s) * (df$f1 + df$f2 + df$f3)
为了比较时间我电脑上的时间:
#Your code
# user system elapsed
# 0.358 0.004 0.362
#@minem
# user system elapsed
# 0.090 0.003 0.093