R matrix/array 和积变换

R matrix/array sumproduct transformation

我正在尝试对两个矩阵进行矩阵和积计算。和积是一个矩阵中的行与另一个矩阵中的列。

set.seed(123)
y <- matrix(sample(1:6,6,FALSE), nrow=3, ncol=2) 
x <- matrix(sample(1:5,15,TRUE), nrow=5, ncol=3)

> x
     [,1] [,2] [,3]
[1,]    5    4    3
[2,]    2    3    2
[3,]    2    3    3
[4,]    4    3    1
[5,]    3    5    5

> y
     [,1] [,2]
[1,]    3    4
[2,]    6    5
[3,]    2    1

结果中 [1,1] 中的单元格值将是 x[ 1 , ] 和 y[ , 1] 的和积。 [1,2] 中的单元格值是 x[ 1 , ] 和 y[ , 2] 的和积。 [2,1] 中的单元格值是 x[ 2 , ] 和 y[ , 1] 的和积。等等。最终结果应该是这样的

> result

         [,1] [,2] 
    [1,]   45   43 
    [2,]   28   25 
    [3,]   30   26 
    [4,]   32   32 
    [5,]   49   42 

我可以使用循环来完成此操作,但如果有一个函数可以自动执行此操作,那就容易多了。

我建议采用下一种方法。有一个名为 crossprod() 的函数,但适用于定义的矩阵维度。你可以检查一下。下一个代码可能有用:

set.seed(123)
#Data
y <- matrix(sample(1:6,6,FALSE), nrow=3, ncol=2) 
x <- matrix(sample(1:5,15,TRUE), nrow=5, ncol=3)

矩阵:

x
     [,1] [,2] [,3]
[1,]    5    5    1
[2,]    4    3    1
[3,]    1    3    5
[4,]    2    1    3
[5,]    3    4    2

y
     [,1] [,2]
[1,]    3    4
[2,]    6    5
[3,]    2    1

代码:

#Code
z <-   t(sapply(1:nrow(x), function(i){
    x[i, ] %*% sapply(1:ncol(y), function(j) {y[,j]})}))

输出:

z
     [,1] [,2]
[1,]   47   46
[2,]   32   32
[3,]   31   24
[4,]   18   16
[5,]   37   34

这是一种使用 outer 对矩阵索引进行运算的方法。 Vectorize 允许我们在 outer 内合并结果。如果您需要性能,我真的认为使用 循环会相对容易。

set.seed(123)
y <- matrix(sample(1:6,6,FALSE), nrow=3, ncol=2) 
x <- matrix(sample(1:5,15,TRUE), nrow=5, ncol=3)

outer(seq_len(nrow(x)),
      seq_len(ncol(y)),
      FUN  = Vectorize(function(i, j) sum(x[i, ] * y[, j])))
#>      [,1] [,2]
#> [1,]   47   46
#> [2,]   32   32
#> [3,]   31   24
#> [4,]   18   16
#> [5,]   37   34

出于好奇,花了大约 5 分钟的时间想出了一个 解决方案:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerMatrix custom_crossprod(IntegerMatrix x, IntegerMatrix y) {
  int nrow = x.rows();
  int ncol = y.cols();
  IntegerMatrix ans(nrow, ncol);
  
  for (int i = 0; i < nrow; i++) {
    IntegerVector x_tmp = x(i, _);
    for (int j = 0; j < ncol; j++) {
      ans(i, j) = sum(x_tmp * y(_, j));
    }
  }
  return(ans);
}

基准点::

bench::mark(
  use_outer = 
outer(seq_len(nrow(x)),
      seq_len(ncol(y)),
      FUN  = Vectorize(function(i, j) sum(x[i, ] * y[, j])))
, custom_crossprod(x, y)
)
# A tibble: 2 x 13
##  expression                 min  median `itr/sec` mem_alloc `gc/sec`
##  <bch:expr>             <bch:t> <bch:t>     <dbl> <bch:byt>    <dbl>
##1 use_outer              179.8us 191.2us     4394.   11.73KB     4.37
##2 custom_crossprod(x, y)   4.9us   6.6us   135309.    2.49KB     0  

@user20650 提供了最合适的答案。

看来 x %*% y 就足够了

set.seed(123)
y <- matrix(sample(1:6,6,FALSE), nrow=3, ncol=2) 
x <- matrix(sample(1:5,15,TRUE), nrow=5, ncol=3)
x %*% y

> x %*% y
     [,1] [,2]
[1,]   47   46
[2,]   32   32
[3,]   31   24
[4,]   18   16
[5,]   37   34

感谢大家的回复。您的编码让我学到了很多关于 R 的工作原理以及如何为某些功能编写脚本的知识。谢谢