加快嵌套 for 循环答案

【问题标题】：Speed up nested for loop加快嵌套 for 循环
【发布时间】：2016-03-04 07:31:07
【问题描述】：

我试图编写一个带有双循环的函数，因为我必须模拟一个双和。这段代码我的代码可以运行，但速度很慢，如何加快速度？

a & y 是向量，x 是矩阵。它们都是相同的长度，即 100。X 为 100x4（100 行，4 列）

X1 <- matrix(rnorm(4*100), ncol=4)
y1 <- sign(X1[,1] + X1[,2] > 0)*2 - 1

fn <- function (a,x,y){
  dsum <-0
  for(i in 1:length(y)){
    for(j in 1:length(y)){
      dsum <- dsum + a[j]*a[i]*y[j]*y[i]*(t(x)[,j])%*%x[i,]
    }
  }
  res <- sum(a)-.5*dsum
  return (res)
}

我尝试过 sum(a)-.5%*%sum(a%o%X%o%y*a%o%t(X)%o%y) 但我肯定错了。

【问题讨论】：

想提供更多细节？ a 和 y 是否共享相同的长度？ x 的维度是多少？您能解释一下您要的金额是多少吗？
它们的长度都是一样的，即100。X是100x4（100行，4列）
How to make a great R reproducible example?
请在您的帖子中添加更多详细信息，而不是在 cmets 中。
t(x)[,j]中的元素与t[j,]中的元素相同

标签： r

【解决方案1】：

试试这个：

fn2<-function(a,x,y) {
   one<-outer(a*y,a*y)
   two<-rowSums(x[rep(1:nrow(x),nrow(x)),]*x[rep(1:nrow(x),each=nrow(x)),])
   sum(a)-.5*sum(one*two)
}

一个例子：

set.seed(1)
a<-runif(100)
y<-runif(100)
x<-matrix(runif(400),nrow=100)
fn(a,x,y)
#[1,] -303.6947
fn2(a,x,y)
#[1] -303.6947

【讨论】：

【解决方案2】：

dsum <- sum(tcrossprod(x*(a*y), x*(a*y)))

给出与循环中的计算相同的结果。简短的版本是：

dsum <- sum(tcrossprod(x*(a*y)))

我想这个解决方案很快。如果你想把它放在你的函数中，你可以这样做：

fn <- function (a,x,y) sum(a) - 0.5*sum(tcrossprod(x*(a*y)))

例子：

x <- matrix(1, 6, 2)
a <- 1:6
y <- 11:16

dsum <- 0
for(i in 1:length(y)) {
  for(j in 1:length(y)) {
    dsum <- dsum + a[j]*a[i]*y[j]*y[i]*(t(x)[,j])%*%x[i,]
  }
}
dsum

sum(tcrossprod(x*(a*y), x*(a*y)))
sum(tcrossprod(x*(a*y)))

【讨论】：

使用 MRO 或任何与良好 BLAS 相关的发行版，这会飞。根据我笔记本电脑上的microbenchmark，几乎比原来的答案快两个数量级，一个比另一个答案快。
我在某种意义上“重写”了函数第二行中的矩阵乘积，后来意识到一个简单的乘积就足够了。 OP应该接受这个答案而不是我的。当然 +1。
@nicola 感谢您的评论。我还需要几分钟来分析代数。