大型矢量距離計算[性能]

Question

嗨，大家好我有一個問題我有一些dimmentions（4-6集羣）中的聚類中心和一個非常大的數據集，我需要將每行分配到最近的集羣。所以它不是一個真正的距離問題，而是性能問題，我的代碼如下：

distances <- matrix(NA, nrow = nrow(ClusterCenters), ncol = nrow(data)) 
calcData <- data[, colnames(ClusterCenters), drop=FALSE] 
for(i in 1:nrow(ClusterCenters)) { 
    distances[i,] <- (rowSums((matrix(unlist(apply(calcData, 1, function(x) {x - ClusterCenters[i,]})), ncol = ncol(calcData), byrow = TRUE))^2))^0.5 
} 
ClusterMemberships <- vector(mode="numeric", length=nrow(calcData)) 
for(i in 1: nrow(calcData)) { 
    ClusterMemberships[i] <- which.min(distances[,i]) 
} 
return(ClusterMemberships) 

有沒有辦法加快速度？我在Windows服務器上工作。

Answer 1

爲了針對六組50 X 1萬個數據行矩陣匹配每50個值，我得到的結果在大約3秒鐘：

vals <- 50 
clusts <- 6 
clusters <- matrix(runif(vals * clusts), nrow=clusts) 

data.count <- 1e6 # large number 
data <- matrix(runif(data.count * vals), nrow=data.count) 

system.time({ 
    dists <- apply(clusters, 1, function(x) rowSums((data - x)^2)^.5) 
    min.dist <- max.col(-dists, ties.method="first") 
}) 
# user system elapsed 
# 2.96 0.47 3.49 

最關鍵的事情是要確保我們限制R函數調用的數量變得很貴。請注意，我如何在羣集上（其中只有六個羣集）而不是數據行（其中有一百萬個）。然後，我使用回收來計算每個集羣與整個集合的距離（注意：與您的數據相比，data已轉置，行數與集羣中的項數一樣多;這對於回收工作是必需的）。

感謝@ user20650提供max.col件。

Answer 2

有幾種方法可以優化R的性能，例如使用BrodieG顯示矢量化。

或者，您可以利用現有計算模式的性能優勢，例如矩陣乘法，排序。

[Michael McCool等，Structured Parallel Programming - Patterns for Efficient Computation]。

而且我們還可以從多核CPU或GPU的這些現有模式的並行庫中獲得進一步的性能優勢。在這種情況下，我們可以用矩陣運算或KNN來表示計算。

1.仿形

通過system.time爲一個大的數據集輸入（1E6），我們可以看到第一環路，其計算兩個向量之間的距離的概要分析這一段代碼，佔95總計算時間的百分比。所以，我們的優化將從這裏開始。

# Original code 
vals <- 50 
clusts <- 6 
ClusterCenters <- matrix(runif(vals * clusts), nrow=clusts) 

data.count <- 1e6 # large number 
calcData <- matrix(runif(data.count * vals), nrow=data.count) 

system.time({ 
    for(i in 1:nrow(ClusterCenters)) { 
    dists[i,] <- (rowSums((matrix(unlist(apply(calcData, 1, function(x) {x  ClusterCenters[i,]})), ncol = ncol(calcData), byrow = TRUE))^2))^0.5 
    } 
}) 
user system elapsed 
71.62 1.13 73.13 

system.time({ 
    for(i in 1: nrow(calcData)) { 
    ClusterMemberships[i] <- which.min(dists[,i]) 
    } 
}) 
user system elapsed 
5.29 0.00 5.31 

2.矢量

矢量是加速R代碼裏面尤其關於「環」，如圖@BrodieG絕對有用的方法。順便說一句，我已經修改了一些在他的解決方案，以獲得正確的結果如下，它可以獲得約3-5X比原來的代碼加速。

#Vectorization: From BrodieG 
dists1 <-matrix(NA, nrow = nrow(ClusterCenters), ncol = nrow(calcData)) system.time({ 
    dists1 <- apply(ClusterCenters, 1, function(x) rowSums(sweep(calcData, 2,x, '-')^2)^.5) 
    min.dist.vec <- max.col(-dists1, ties.method="first") 
}) 
user system elapsed 
16.13 1.42 17.61 
all.equal(ClusterMemberships, min.dist.vec) 
[1] TRUE 

3.矩陣圖案

然後，讓我們來回顧第一環路，它通過總結列計算的距離（calcData [I，] - ClusterCenters [J，]）^ 2。

所以，我們可以將這種操作通過如下擴展該方程爲矩陣：

calcData [I，]^2 - 2 * calcData [I，] * ClusterCenters [J，] + ClusterCenters [J，]^2

因此，對於第一和第三部分中，我們可以做簡單的矩陣乘法掃，如

calcData * calcData

而對於第二項，我們需要矩陣傳遞的特技技能，然後它改變的

ClusterCenters％*％噸的矩陣乘法（calcData）

最後，整個代碼的矩陣運算如下：

# Pattern Representation 1: Matrix 
dists2 <-matrix(NA, nrow = nrow(ClusterCenters), ncol = nrow(calcData)) 
system.time({ 
    data2  <- rowSums(calcData*calcData) 
    clusters2 <- rowSums(ClusterCenters*ClusterCenters) 
    # NVIDIA GPU: nvBLAS can speedup this step 
    # Futher Info on ParallelR.com 
    ClustersXdata <- calcData %*% t(ClusterCenters) 
    # compute distance 
    dists2 <- sweep(data2 - 2 * ClustersXdata, 2, clusters2, '+') ^0.5 
    min.dist.matrix <- max.col(-dists2, ties.method="first") 
}) 
user system elapsed 
1.17 0.09 1.28 
all.equal(ClusterMemberships, min.dist.matrix) 
[1] TRUE 

現在我們可以看到所有這三部分都可以通過矩陣運算來完成。通過這種方法，計算時間幾乎是線性的，從10^3到10^7，並且其約爲1e6 calcData集合的原始代碼更快。

4. KNN

在這種情況下，計算可被視爲找到在簇數據集第一最近鄰所以這是一個種簡單的KNN，其中k = 1。而且我們可以很容易地使用由C編寫的類包中的knn函數。同時它也會非常高效。如下面的示例代碼：

# Pattern Representation 2: KNN 
library("class") 
system.time(
    min.dist.knn <- knn(ClusterCenters, calcData, cl = 1:nrow(ClusterCenters), k = 1) 
) 
user system elapsed 
1.21 0.12 1.35 

all.equal(ClusterMemberships, as.integer(min.dist.knn)) 
[1] TRUE 

KNN的運行時間是作爲我們的1E6矩陣運算代碼相似，但如果我們申請更多的大數據這兩種算法，我們可以看到矩陣算法仍然贏了，矩陣算法是快於KNN（15.9 .vs。29.1）的2X。

終於，在這篇文章中，我只是展示了一些性能調優的想法，我們可以繼續微調這些代碼，包括架構指定的優化和使用c/C++來重寫它。甚至可以在多核CPU和NVIDIA GPU上並行化矩陣運算和K​​NN，您可以參考ParallelR