Masi Masi - 1 month ago 10
R Question

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

I want to do a simple column (Nx1) times row (1xM) multiplication, resulting in (NxM) matrix.
Code where I create a row by sequence, and column by transposing a similar sequence

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


Output which indicates that R thinks matrices more like columns

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


Expected output: NxM matrix.

OS: Debian 8.5

Linux kernel: 4.6 backports

Hardware: Asus Zenbook UX303UA

Answer

In this case using outer would be a more natural choice

outer(1:6, 1:6)

In general for two numerical vectors x and y, the matrix rank-1 operation can be computed as

outer(x, y)

If you want to resort to real matrix multiplication routines, use tcrossprod:

tcrossprod(x, y)

If either of your x and y is a matrix with dimension, use as.numeric to cast it as a vector first.

It is highly not recommended to use general matrix multiplication operation "%*%" for this. But if you want, make sure you get the dimension comformable, so that x is a single-column matrix and y is a single-row matrix: x %*% y.


Can you say anything about efficiency?

Matrix rank-1 operation is known to be memory-bound. So make sure we use gc() for garbage collection to tell R to release memory from heap after every replicate (otherwise your system will stall):

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 

In terms of performance, they are all very alike.

Comments