Help porting broadcasting code to R

[pdeffebach]

I have the following code in Julia

julia> x = [1 2];

julia> y = reshape([3,4], (2, 1));

julia> x .+ y
2×2 Array{Int64,2}:
 4  5
 5  6

In julia, the . is broadcasting and is equivalent to

julia> x_expanded = vcat(x, x)
2×2 Array{Int64,2}:
 1  2
 1  2

julia> y_expanded = hcat(y, y)
2×2 Array{Int64,2}:
 3  3
 4  4

julia> x_expanded + y_expanded
2×2 Array{Int64,2}:
 4  5
 5  6

How would I do the same in R? After some benchmarking in Julia the broadcasted summation looks like the fastest option for my code. I am wondering if there is some sort of highly optimized equivalent in R.

r$> x = matrix(c(1, 2), nrow = 1, ncol = 2)                                                               

r$> y = matrix(c(3, 4), nrow = 2, ncol = 1)                                                               

r$> x + y                                                                                                 
Error in x + y : non-conformable arrays

Can someone point me in the right direction?

[FJCC]

I don't know if this is highly optimized.

x = matrix(c(1, 2), nrow = 1, ncol = 2)                                                               

y = matrix(c(3, 4), nrow = 2, ncol = 1)

sapply(x, function(N) N + y)
#>      [,1] [,2]
#> [1,]    4    5
#> [2,]    5    6

^{Created on 2020-11-07 by the reprex package (v0.3.0)}

[xiaodai]

x = t(matrix(c(1, 2), nrow = 2, ncol = 2))

y = matrix(c(3, 4), nrow = 2, ncol = 2)                                                               

x + y  

[pdeffebach]

Thanks.

However this will allocate the full matrix. In my use case if I were to expand all the dimensions I would be working with 370 x 370 x 370 x 5 size matrices, which are very big.

Additionally, I need to broadcast along the 3rd dimension. Does broadcasting only work along the first dimension in R?

[edavidaja]

Take a look at the {rray} vignette on broadcasting, which might be helpful.

[pdeffebach]

Unfortunately, it seems like it's not available for R 3.6.3, which is what I'm using.

It seems to have been removed from CRAN. Does anyone know the reason for this? Should I just use devtools to install it?

[edavidaja]

you should be fine to install it from GitHub.

[jmcvw]

In general, I believe, R is well optimized for working with matrices, and has oodles of functions for doing so.
Maybe outer(1:2, 3:4, '+') is what you need, no extra packages needed

[pdeffebach]

Hello everyone, thank you for all your replies.

rray is the winner! It works really well. I am very happy with performance. It's within a factor of 2 of Julia after replacing naive "manually expand matrices and multiply" approach.

Here is a MWE showing an example of how it's used

library(rray)
mwe = function() {
  n = 10
  k = 5
  e = array(1:100, dim = c(n, n, 1, 1))
  p = array(1:50, dim = c(1, 1, n, k))

  # need output of dimension n n n k 
  full_dims = c(n, n, n, k)
  # Current approach
  e_full = array(NA, dim = full_dims)
  p_full = array(NA, dim = full_dims)
  for(ni in 1:n){
    for(ki in 1:k) {
      e_full[ , , ni, ki] = e
    }
  }

  for(ni in 1:n){
    for(nj in 1:n) {
      p_full[ni, nj, , ] = p
    }
  }

  correct_out = p_full + e_full


  t = sapply(e, function(x) x + p)
  # This is likely performant, but gives the wrong answer! Not clear
  # how to fil in the array the right way. 
  fast_out = array(t, dim = c(n, n, n, k))

  e_b = rray_broadcast(e, dim = full_dims)
  p_b = rray_broadcast(p, dim = full_dims)

  rray_out = e_b + p_b
}