CODEWITHSUNDEEP
r
Abhirup Datta
Abhirup Datta

Reputation: 103

Outer not working in R

I have 3 vectors x, y and z. I want to create a matrix whose (i,j)^th entry will be sum((x[i]*z+y[j])^2)

I'm trying to use outer for this but I'm getting an error. Here is an example:

x=y=1:10/10
z=1:10
outer(x,y,FUN=function(x,y,z) sum((x*z+y)^2),z)

Error in outer(x, y, FUN = function(x, y, z) sum((x * z + y)^2), z) : 
  dims [product 100] do not match the length of object [1]

Upvotes: 0

Views: 274

Answers (1)

josliber
josliber

Reputation: 44330

You have the following summation:

\sum_{k=1}^n [(x_i*z_k + y_j)^2]
\sum_{k=1}^n [(x_i*z_k)^2 + 2*x_i*z_k*y_j + y_j^2]
\sum_{k=1}^n [(x_i*z_k)^2] + 2*x_i*y_j*\sum_{k=1}^n [z_k] + n*y_j^2

So an efficient way to proceed would be to compute the first term with an outer product between x and z and to compute the second term with an outer product between x and y:

outer.sol <- function(x, y, z) {
  n <- length(x)
  xmult <- rowSums(outer(x, z, function(x, z) (x*z)^2))
  matrix(rep(xmult, n), nrow=n) +
    outer(x, y, function(x, y) 2*sum(z)*x*y) +
    matrix(rep(n*y^2, each=n), nrow=n)
}
outer.sol(x, y, z)
#         [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]
#  [1,]   5.05   6.45   8.05   9.85  11.85  14.05  16.45  19.05  21.85  24.85
#  [2,]  17.70  20.20  22.90  25.80  28.90  32.20  35.70  39.40  43.30  47.40
#  [3,]  38.05  41.65  45.45  49.45  53.65  58.05  62.65  67.45  72.45  77.65
#  [4,]  66.10  70.80  75.70  80.80  86.10  91.60  97.30 103.20 109.30 115.60
#  [5,] 101.85 107.65 113.65 119.85 126.25 132.85 139.65 146.65 153.85 161.25
#  [6,] 145.30 152.20 159.30 166.60 174.10 181.80 189.70 197.80 206.10 214.60
#  [7,] 196.45 204.45 212.65 221.05 229.65 238.45 247.45 256.65 266.05 275.65
#  [8,] 255.30 264.40 273.70 283.20 292.90 302.80 312.90 323.20 333.70 344.40
#  [9,] 321.85 332.05 342.45 353.05 363.85 374.85 386.05 397.45 409.05 420.85
# [10,] 396.10 407.40 418.90 430.60 442.50 454.60 466.90 479.40 492.10 505.00

While you could use a looping solution like mapply to compute each element individually, which would make your code a good deal easier to read, note that this will be much slower than the fully vectorized solution with outer:

mapply.sol <- function(x, y, z) {
  n <- length(x)
  matrix(mapply(function(i,j) sum((i*z+j)^2), rep(x, n), rep(y, each=n)), nrow=n)
}

# 1000 x 1000 case
x <- y <- z <- 1:1000
all.equal(outer.sol(x, y, z), mapply.sol(x, y, z))
# [1] TRUE
system.time(outer.sol(x, y, z))
#    user  system elapsed 
#   0.060   0.001   0.060 
system.time(mapply.sol(x, y, z))
#    user  system elapsed 
#  18.522   1.986  20.432

From the comments it appears the true function is more complicated than the quadratic one included in the question, so the provided mapply.sol with your function is the (rather slow) way to go.

Upvotes: 1

Related Questions