如何对矩阵 NxM 与 Nx1 1xM 进行 R 乘法？

Question

我想做一个简单的列 (Nx1) 乘以行 (1xM) 乘法，得到 (NxM) 矩阵。我按顺序创建行，并通过转置类似序列

创建列的代码

row1 <- seq(1:6) 
col1 <- t(seq(1:6))      
col1 * row1

表明 R 认为矩阵更像列的输出

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    4    9   16   25   36

预期输出：NxM 矩阵。

OS：Debian 8.5
Linux 内核：4.6 向后移植
硬件：华硕 Zenbook UX303UA

Answer 1

在这种情况下，使用 outer 将是更自然的选择

outer(1:6, 1:6)

一般来说，对于两个 数值向量 x 和 y，矩阵 rank-1 运算可以计算为

outer(x, y)

如果您想求助于实矩阵乘法例程，请使用 tcrossprod:

tcrossprod(x, y)

如果您的 x 和 y 中的任何一个是具有维度的矩阵，请先使用 as.numeric 将其转换为向量。

不建议为此使用一般的矩阵乘法运算"%*%"。但是如果你想要，确保你得到合适的维度：x 是一个单列矩阵，y 是一个单行矩阵，所以 x %*% y.

Can you say anything about efficiency?

已知矩阵 rank-1 操作受内存限制。因此，请确保我们使用 gc() 进行垃圾收集，以告诉 R 在每次复制后从堆中释放内存（否则您的系统将停止）：

x <- runif(500)
y <- runif(500)
xx <- matrix(x, ncol = 1)
yy <- matrix(y, nrow = 1)

system.time(replicate(200, {outer(x,y); gc();}))
#   user  system elapsed 
#  4.484   0.324   4.837 

system.time(replicate(200, {tcrossprod(x,y); gc();}))
#   user  system elapsed 
#  4.320   0.324   4.653 

system.time(replicate(200, {xx %*% yy; gc();}))
#   user  system elapsed 
#  4.372   0.324   4.708

在性能方面，它们都非常相似。

跟进

当我回来时，我看到了另一个基准不同的答案。好吧，问题是，这取决于问题的大小。如果您只是尝试一个小示例，则无法消除所有三个函数的函数解释/调用开销。如果你这样做

x <- y <- runif(500)
microbenchmark(tcrossprod(x,y), x %*% t(y), outer(x,y), times = 200)

您将再次看到大致相同的性能。

#Unit: milliseconds
#             expr     min      lq     mean  median      uq      max neval cld
# tcrossprod(x, y) 2.09644 2.42466 3.402483 2.60424 3.94238 35.52176   200   a
#       x %*% t(y) 2.22520 2.55678 3.707261 2.66722 4.05046 37.11660   200   a
#      outer(x, y) 2.08496 2.55424 3.695660 2.69512 4.08938 35.41044   200   a

Answer 2

下面是当使用的向量长度为 100 时三种方法的执行速度比较。最快的是 tcrossprod，x%*%t(y) 需要 17% 的时间，outer(x,y) 花费 45% 的时间（中位数时间）。在table中，neval是函数被评估的次数以获得基准分数。

> x <- runif(100,0,100)
> y <- runif(100,0,100)
> microbenchmark(tcrossprod(x,y), x%*%t(y), outer(x,y), times=5000)
Unit: microseconds
             expr    min      lq     mean  median      uq       max neval
 tcrossprod(x, y) 11.404 16.6140 50.42392 17.7300 18.7555  5590.103  5000
       x %*% t(y) 13.878 19.4315 48.80170 20.5405 21.7310  4459.517  5000
      outer(x, y) 19.238 24.0810 72.05250 25.3595 26.8920 89861.855  5000

要得到下图，有

library("ggplot2")
bench <- microbenchmark(tcrossprod(x,y), x%*%t(y), outer(x,y), times=5000)
autplot(bench)

编辑：性能取决于x和y的大小，当然还有机器运行的代码。我最初使用长度为 100 的向量进行基准测试，因为这是 Masi 所询问的。然而，这三种方法似乎对较大的向量具有非常相似的性能。对于长度为 1000 的向量，三种方法的中位数时间在我的机器上彼此相差在 5% 以内。

> x <- runif(1000)
> y <- runif(1000)
> microbenchmark(tcrossprod(x,y),x%*%t(y),outer(x,y),times=2000)
Unit: milliseconds
             expr      min       lq     mean   median       uq       max neval
 tcrossprod(x, y) 1.870282 2.030541 4.721175 2.916133 4.482346  75.77459  2000
       x %*% t(y) 1.861947 2.067908 4.921061 3.067670 4.527197 105.60500  2000
      outer(x, y) 1.886348 2.078958 5.114886 3.033927 4.556067  93.93450  2000

Answer 3

查看此内容的一种简单方法是将向量转换为矩阵

row1.mat = matrix(row1)
col1.mat = matrix(col1)

然后使用dim查看矩阵的维数：

dim(row1.mat)
dim(col1.mat)

如果您希望产品为此工作，您需要一个 6*1 矩阵，乘以一个 1*6 矩阵。所以你需要使用 t(col1.mat).

转置 col1.mat

你可能知道矩阵乘积是 %*%

row1.mat %*% t(col1.mat)

此方法与其他方法的比较

library("microbenchmark")
x <- runif(1000)
y <- runif(1000)
xx = matrix(x)
yy = matrix(y)
microbenchmark(tcrossprod(x,y),x%*%t(y),outer(x,y), xx %*% t(yy), times=2000)

Unit: milliseconds
             expr      min       lq     mean   median       uq      max neval
 tcrossprod(x, y) 2.829099 3.243785 6.015880 4.801640 5.040636 77.87932  2000
       x %*% t(y) 2.847175 3.251414 5.942841 4.810261 5.049474 86.53374  2000
      outer(x, y) 2.886059 3.277811 5.983455 4.788054 5.074997 96.12442  2000
     xx %*% t(yy) 2.868185 3.255833 6.126183 4.699884 5.056234 87.80024  2000

如何对矩阵 NxM 与 Nx1 1xM 进行 R 乘法？

How to do R multiplication with Nx1 1xM for Matrix NxM?

statistics

r

matrix

matrix-multiplication