¿Cómo puedo vectorizar este código en R? ¿Quizás con la función Aplicar ()?

Question

Realmente estoy luchando por replicar la salida de la función Dist () en el código R sin usar 1 o 2 para bucles. (Si te preguntas por qué estoy haciendo esto, es para que pueda jugar con el cálculo de distancia y también para mejorar mis habilidades R, así que solo soluciones que involucran a R!)

Descripción general: la matriz se pasa a dist (), que calcula la distancia euclidiana en cuanto a fila y genera una matriz de distancia completa de la distancia entre cada fila (por ejemplo, la distancia entre las filas 1 y 50 estará en distancematrix [1, 50] y distancematrix [ 50, 1]). El código rápido se ve así:

distancematrix <- as.matrix(dist(myMatrix, method="euclidean", diag = T, upper = T))

He producido con éxito la misma salida en R usando el siguiente código:

for (i in 1:nrow(myMatrix)) {
  for (j in 1:nrow(myMatrix)) {
    distancematrix[i, j] <- sum(abs(myMatrix[i,] - myMatrix[j,])) 
  }
}

Sin embargo, usar dos bucles anidados es mucho más lento que usar dist (). He leído mucho sobre el uso de aplicar () para optimizar más lento para los bucles, pero no he podido entenderlo hasta ahora. Creo que al menos uno de los bucles For es definitivamente evitable simplemente generando un vector y lidiar con él al final. Sin embargo, por mi vida no puedo resolver cómo eliminar ambos para los bucles.

Alguien tiene alguna opinión?

Answer 1

En primer lugar, debe tenerse en cuenta que el código que publicó en realidad no replica la salida del dist función, porque la línea:

distancematrix[i, j] <- sum(abs(myMatrix[i,] - myMatrix[j,]))

no calcula la distancia euclidiana; debería ser:

distancematrix[i, j] <- sqrt(sum((myMatrix[i,] - myMatrix[j,]) ^ 2))

Aquí hay dos soluciones en las que confían apply. Se simplifican y, en particular, no aprovechan la simetría de la matriz de distancia (que, si se considera, conduciría a una aceleración de 2 veces). Primero, genere algunos datos de prueba:

# Number of data points
N <- 2000
# Dimensionality
d <- 10
# Generate data
myMatrix = matrix(rnorm(N * d), nrow = N)

Por conveniencia, defina:

# Wrapper for the distance function
d_fun <- function(x_1, x_2) sqrt(sum((x_1 - x_2) ^ 2))

El primer enfoque es una combinación de apply y sapply:

system.time(
    D_1 <- 
        apply(myMatrix, 1, function(x_i) 
            sapply(1:nrow(myMatrix), function(j) d_fun(x_i, myMatrix[j, ]))
    )
)

   user  system elapsed 
 14.041   0.100  14.001 

Mientras que el segundo usa solo apply (pero repasando los índices, que se combinan usando expand.grid):

system.time(
    D_2 <- 
        matrix(apply(expand.grid(i = 1:nrow(myMatrix), j = 1:nrow(myMatrix)), 1, function(I) 
            d_fun(myMatrix[I[["i"]], ], myMatrix[I[["j"]], ])
        )
    )
)

   user  system elapsed 
 39.313   0.498  39.561 

Sin embargo, como se esperaba, ambos son mucho más lentos que dist:

system.time(
    distancematrix <- as.matrix(
        dist(myMatrix, method = "euclidean", diag = T, upper = T)
    )
)

   user  system elapsed 
  0.337   0.054   0.388