このコードをRでベクトル化するにはどうすればよいですか？たぶんApply（）関数がありますか？

Question

ループに1または2を使用せずに、rコードのdist（）関数の出力を複製するのに本当に苦労しています。 （なぜ私がこれをしているのか疑問に思っているなら、それは私が距離計算で遊ぶことができるように、そしてRスキルを向上させることができるように - だから、Rを含む解決策のみをお願いします！）

概要：マトリックスはdist（）に渡されます。これにより、ユークリッド距離が列の距離を計算し、各行間の距離の完全距離マトリックスが出力されます（行1と50の間の距離は、distancematrix [1、50]とdistancematrix [ 50、1]）。高速コードは次のようになります：

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

次のコードを使用して、Rで同じ出力を正常に生成しました。

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

ただし、2つのネストされたループを使用すると、dist（）を使用するよりもはるかに遅くなります。私はApply（）を使用してループの最適化を最適化することについて多くを読みましたが、これまでのところ頭を周回することはできませんでした。私は、ベクトルを出力するだけで、最後にそれを扱うことで、少なくとも1つのループの1つは間違いなく回避可能であると信じています。しかし、私の人生のために、ループのために両方を削除する方法を解決することはできません。

誰かが考えもありますか？

Answer 1

まず第一に、あなたが投稿したコードが実際にの出力を実際に複製しないことに注意する必要があります dist 関数、行：

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

ユークリッド距離を計算しません。そのはず：

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

依存する2つのソリューションがあります apply. 。それらは単純化されており、特に距離マトリックスの対称性を利用しません（考慮した場合、2倍のスピードアップにつながります）。まず、いくつかのテストデータを生成します。

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

便宜上、定義してください。

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

最初のアプローチはの組み合わせです apply と 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 

2番目は使用しますが apply （ただし、使用してペアリングされるインデックスを介して進む 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 

ただし、予想どおり、両方ともはるかに遅いです dist:

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

   user  system elapsed 
  0.337   0.054   0.388