Why does this apply function return a different value than the loop?

Question

I am attempting to covert an apply function to a nested loop for use in software outside of R. The loop I have written appears similar, but returns slightly different values than the apply function.

I have tried adjusting the indexing throughout, but still cannot return identical values.

require(MASS)


# set up grid
B <- 5
co <- seq(0, 1, length=B)
Z <- cbind(rep(co, each=B), rep(co, times=B))



# long hand euclidean distance
x <- Z[,1]
y <- Z[,2]
d <- matrix(NA, nrow = nrow(Z), ncol = nrow(Z))
N <- nrow(Z)
for(i in 1:N){
  for(j in 1:N){
    d[i,j] <- sqrt((x[i] - x[j])*(x[i] - x[j]) + (y[i] - y[j])*(y[i] - y[j]))
  }
}

# generate a value for each pixel
cov <- MASS::mvrnorm(1, mu = rep(0, nrow(Z)), Sigma = exp(-d/0.25))

# calculate weights of each pixel
# apply
applyW <- apply(d, 1, function(x) {
  w0 <- exp(-x^2 / (2*0.3 * 0.3))
  w0[which.min(x)] <- 0
  w <- w0/sum(w0)
  sum(cov * w)
})


# for loop approach
w0 <- matrix(NA, ncol = ncol(d), nrow = nrow(d))
w <- matrix(NA, ncol = ncol(d), nrow = nrow(d))
loopW <- c()
for(i in 1:nrow(d)){
  for(j in 1:ncol(d)){
    w0[i, j] <- exp(-(d[i,j] * d[i,j]) / (2 * 0.3 * 0.3))
    w0[which.min(w0)] <- 0
  }
  w[i, ] <- w0[i ,]/sum(w0[i,])
  loopW[i] <- sum(cov * w[i,])
}
# check the values
cbind(applyW, loopW)

I am expecting applyW and loopW to be identical values, but I am not finding that. Values are similar, but not the same. Is this just a rounding error or is there an error in the code?

Answer 1

I have not spent the time to figure out exactly what you are trying to accomplish, but these two parts of your code do not match in terms of what they are doing:

In your apply:

w0[which.min(x)] <- 0

Here w0 is a vector (i.e., one row) and you are replacing the minimum value in that vector with a 0. This is done for every row separately as per the apply(MARGIN = 1).

However, in your loops you have:

w0[which.min(w0)] <- 0

Here w0 is a matrix, and the location in the code where you have w0[which.min(w0)] <- 0 continually replaces the minimum value with 0 every time you move to the next matrix element (i.e., not per row as in your apply()).

If you move this line outside the inner loop and make some tweaks to match the apply(), you will get identical values:

require(MASS)

# set up grid
B <- 5
co <- seq(0, 1, length=B)
Z <- cbind(rep(co, each=B), rep(co, times=B))

# long hand euclidean distance
x <- Z[,1]
y <- Z[,2]
d <- matrix(NA, nrow = nrow(Z), ncol = nrow(Z))
N <- nrow(Z)
for(i in 1:N){
  for(j in 1:N){
    d[i,j] <- sqrt((x[i] - x[j])*(x[i] - x[j]) + (y[i] - y[j])*(y[i] - y[j]))
  }
}

# generate a value for each pixel
cov <- MASS::mvrnorm(1, mu = rep(0, nrow(Z)), Sigma = exp(-d/0.25))

# calculate weights of each pixel
# apply
applyW <- apply(d, 1, function(x) {
  w0 <- exp(-x^2 / (2*0.3 * 0.3))
  w0[which.min(x)] <- 0
  w <- w0/sum(w0)
  sum(cov * w)
})


# for loop approach
w0 <- matrix(NA, ncol = ncol(d), nrow = nrow(d))
w <- matrix(NA, ncol = ncol(d), nrow = nrow(d))
loopW <- c()
for(i in 1:nrow(d)){
  for(j in 1:ncol(d)){
    w0[i, j] <- exp(-(d[i,j] * d[i,j]) / (2 * 0.3 * 0.3))
    #w0[which.min(w0)] <- 0 move this outside inner loop
  }
  w0[i, which.min(d[i, ])] <- 0 # and make some tweaks to match the apply
  w[i, ] <- w0[i ,]/sum(w0[i,])
  loopW[i] <- sum(cov * w[i,])
}
# check the values
output <- cbind(applyW, loopW)
output <- cbind(output, "dif" = output[, "applyW"] - output[, "loopW"])
output

Produces:

           applyW       loopW dif
 [1,]  0.62023939  0.62023939   0
 [2,]  0.45082535  0.45082535   0
 [3,]  0.07066739  0.07066739   0
 [4,]  0.01518537  0.01518537   0
 [5,]  0.13383600  0.13383600   0
 [6,]  0.52896297  0.52896297   0
 [7,]  0.38474014  0.38474014   0
 [8,]  0.23321748  0.23321748   0
 [9,]  0.07092092  0.07092092   0
[10,] -0.04388823 -0.04388823   0
[11,]  0.46926149  0.46926149   0
[12,]  0.39166724  0.39166724   0
[13,]  0.12309979  0.12309979   0
[14,] -0.04778011 -0.04778011   0
[15,]  0.10359100  0.10359100   0
[16,]  0.41273927  0.41273927   0
[17,]  0.33122926  0.33122926   0
[18,] -0.10115094 -0.10115094   0
[19,] -0.15051985 -0.15051985   0
[20,] -0.23035539 -0.23035539   0
[21,]  0.37065262  0.37065262   0
[22,] -0.24497145 -0.24497145   0
[23,] -0.04993695 -0.04993695   0
[24,] -0.24990339 -0.24990339   0
[25,] -0.27183745 -0.27183745   0