大型矢量距离计算[性能]

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