加快所有行对的点积并添加截距项
Speed up taking dot product of all pairs of rows and adding intercept term
我正在对 83 个数据点应用 100 个不同的模型模拟,并检查每个数据点的每个模型模拟的估计范围。
每个计算本身都是 342 个变量乘以 342 个系数的乘积,再加上一个截距。我在下面编写的代码可以正确计算这些值,但是速度慢得令人发指。有什么办法可以提高处理速度吗?
spec = read.csv(spectra)
coef = read.csv(coefficents)
shell = matrix(data=NA,ncol=101,nrow=nrow(spec))
shell = as.data.frame(shell)
heading = paste("Sim_",seq(1,100,1),sep="")
names(shell)[1] = "Filename"
names(shell)[2:101] = heading
shell[1] = spec[1]
for (i in 1:nrow(spec))
{
for (j in 1:100)
{
shell[i,j+1] = sum(spec[i,2:341]*coef[j,3:342]) + coef[j,2]
}
}
实际上,您正在执行 spec
的矩阵乘法与 coef
的转置,然后向每一列添加一个常数。您应该通过使用内置矩阵乘法函数 %*%
和向量化操作来缩放列来获得加速:
out <- cbind(spec[,1], t(t(spec[,2:341] %*% t(coef[1:100,3:342])) + coef[1:100,2]))
此 1-liner 的输出与原始 post 中代码的输出相同(稍作修改以获取和输出矩阵而不设置任何名称):
OP <- function(spec, coef) {
shell = matrix(data=NA,ncol=101,nrow=nrow(spec))
shell[,1] <- spec[,1]
for (i in 1:nrow(spec)) {
for (j in 1:100) {
shell[i,j+1] = sum(spec[i,2:341]*coef[j,3:342]) + coef[j,2]
}
}
shell
}
all.equal(out, OP(spec, coef))
# [1] TRUE
就运行时间而言,矢量化操作产生了显着的好处 (38x),即使对于这个相对较小的示例(spec
中有 1000 行)也是如此:
system.time(cbind(spec[,1], t(t(spec[,2:341] %*% t(coef[1:100,3:342])) + coef[1:100,2])))
# user system elapsed
# 0.028 0.001 0.030
system.time(OP(spec, coef))
# user system elapsed
# 0.927 0.224 1.161
数据:
set.seed(144)
spec <- matrix(rnorm(1000*341), nrow=1000)
coef <- matrix(rnorm(100*342), nrow=100)
我正在对 83 个数据点应用 100 个不同的模型模拟,并检查每个数据点的每个模型模拟的估计范围。
每个计算本身都是 342 个变量乘以 342 个系数的乘积,再加上一个截距。我在下面编写的代码可以正确计算这些值,但是速度慢得令人发指。有什么办法可以提高处理速度吗?
spec = read.csv(spectra)
coef = read.csv(coefficents)
shell = matrix(data=NA,ncol=101,nrow=nrow(spec))
shell = as.data.frame(shell)
heading = paste("Sim_",seq(1,100,1),sep="")
names(shell)[1] = "Filename"
names(shell)[2:101] = heading
shell[1] = spec[1]
for (i in 1:nrow(spec))
{
for (j in 1:100)
{
shell[i,j+1] = sum(spec[i,2:341]*coef[j,3:342]) + coef[j,2]
}
}
实际上,您正在执行 spec
的矩阵乘法与 coef
的转置,然后向每一列添加一个常数。您应该通过使用内置矩阵乘法函数 %*%
和向量化操作来缩放列来获得加速:
out <- cbind(spec[,1], t(t(spec[,2:341] %*% t(coef[1:100,3:342])) + coef[1:100,2]))
此 1-liner 的输出与原始 post 中代码的输出相同(稍作修改以获取和输出矩阵而不设置任何名称):
OP <- function(spec, coef) {
shell = matrix(data=NA,ncol=101,nrow=nrow(spec))
shell[,1] <- spec[,1]
for (i in 1:nrow(spec)) {
for (j in 1:100) {
shell[i,j+1] = sum(spec[i,2:341]*coef[j,3:342]) + coef[j,2]
}
}
shell
}
all.equal(out, OP(spec, coef))
# [1] TRUE
就运行时间而言,矢量化操作产生了显着的好处 (38x),即使对于这个相对较小的示例(spec
中有 1000 行)也是如此:
system.time(cbind(spec[,1], t(t(spec[,2:341] %*% t(coef[1:100,3:342])) + coef[1:100,2])))
# user system elapsed
# 0.028 0.001 0.030
system.time(OP(spec, coef))
# user system elapsed
# 0.927 0.224 1.161
数据:
set.seed(144)
spec <- matrix(rnorm(1000*341), nrow=1000)
coef <- matrix(rnorm(100*342), nrow=100)