计算R中两个矩阵的所有行组合的函数答案

【问题标题】：Calculate function for all row combinations of two matrices in R计算R中两个矩阵的所有行组合的函数
【发布时间】：2012-05-25 18:24:28
【问题描述】：

我想计算两个矩阵/数据框之间所有行组合的距离度量。

结果将是一个矩阵，其中单元格 i,j 对应于函数给出的结果，该函数应用于第一个矩阵的第 i 行和第二个矩阵的第 j 行。这是一个示例，说明我想用一个示例函数来处理 for 循环。

x<-matrix(rnorm(30),10,3)  ## Example data
y<-matrix(rnorm(12),4,3)

results<-matrix(NA,nrow(x),nrow(y))

for (i in 1:nrow(x)){
  for (j in 1:nrow(y)){
    r1<-x[i,]
    r2<-y[j,]
    results[i,j]<-sum(r1*r2)  ## Example function
  }
}

在现实生活中，我的第一个矩阵有几十万行，第二个矩阵有几百行，我要计算的函数不是点积（我意识到我可能选择了一个函数似乎我想做的只是矩阵乘法）。事实上，我想替换一些函数，所以我想找到一个可推广到不同函数的解决方案。一种思考方式是我想劫持矩阵乘法来执行其他功能。用 for 循环计算这个需要很长时间，这是不切实际的。如果有任何关于如何在没有 for 循环的情况下执行此操作的提示，我将不胜感激。

【问题讨论】：

标签： r

【解决方案1】：

outer(1:nrow(x), 1:nrow(y), Vectorize(function(i, j) sum(x[i, ] * y[j, ])))

【讨论】：

看起来不错。不过，可能不会比 for 循环快，而且可能更慢。
@DWin 是对的。我尝试了这种方法，令我惊讶的是，“外部”加“矢量化”方法所用的时间至少与 for 循环一样长（尽管我不知道在它完成之前我将它切断了多长时间）。最后，我发现我可以将我的函数分解为可以使用应用和矩阵代数计算的部分，并且与约 1 天相比，这可以在几秒钟内计算出一些大型数据帧。但是，我担心我无法使用其他功能做到这一点。

【解决方案2】：

我知道你很久以前就问过这个问题，但我想我可能会与你分享一个解决方案，与for 循环相比，当你拥有的行数变得非常大时，它会变得更有效率。在少量行中，速度差异可以忽略不计（for 循环甚至可能更快）。这仅依赖于子集和 rowSums 的使用，非常简单：

## For reproducibility
set.seed( 35471 )

## Example data - bigger than the original to get and idea of difference in speed
x<-matrix(rnorm(60),20,3)
y<-matrix(rnorm(300),100,3)

# My function which uses grid.expand to get all combinations of row indices, then rowSums to operate on them
rs <- function( x , y ){
rows <- expand.grid( 1:nrow(x) , 1:nrow(y) )
results <- matrix( rowSums( x[ rows[,1] , ] * y[ rows[,2] , ] ) , nrow(x) , nrow(y) )
return(results)
}

# Your orignal function
flp <- function(x ,y){
results<-matrix(NA,nrow(x),nrow(y))
for (i in 1:nrow(x)){
  for (j in 1:nrow(y)){
    r1<-x[i,]
    r2<-y[j,]
    results[i,j]<-sum(r1*r2)  ## Example function
  }
}
return(results)
}


## Benchmark timings:
library(microbenchmark)
microbenchmark( rs( x, y ) , flp( x ,y ) , times = 100L )
#Unit: microseconds
#     expr      min       lq     median        uq      max neval
#  rs(x, y)  487.500  527.396   558.5425   620.486   679.98   100
# flp(x, y) 9253.385 9656.193 10008.0820 10430.663 11511.70   100

## And a subset of the results returned from each function to confirm they return the same thing!
flp(x,y)[1:3,1:3]
#          [,1]       [,2]       [,3]
#[1,] -0.5528311  0.1095852  0.4461507
#[2,] -1.9495687  1.7814502 -0.3769874
#[3,]  1.8753978 -3.0908057  2.2341414

rs(x,y)[1:3,1:3]
#          [,1]       [,2]       [,3]
#[1,] -0.5528311  0.1095852  0.4461507
#[2,] -1.9495687  1.7814502 -0.3769874
#[3,]  1.8753978 -3.0908057  2.2341414

所以你可以看到，通过使用rowSums 和子集，当行组合的数量只有 2000 时，我们可以比 for 循环快 20 倍。如果你有更多，速度上的差异会更大。

HTH。

【讨论】：