DBSCAN聚類算法原理及其實現

jopen 9年前發布 | 40K 次閱讀算法

DBSCAN（Density-Based Spatial Clustering of Applications with Noise）聚類算法，它是一種基于高密度連通區域的、基于密度的聚類算法，能夠將具有足夠高密度的區域劃分為簇，并在具有噪聲的數據中發現任意形狀的簇。我們總結一下DBSCAN聚類算法原理的基本要點：

DBSCAN算法需要選擇一種距離度量，對于待聚類的數據集中，任意兩個點之間的距離，反映了點之間的密度，說明了點與點是否能夠聚到同一類中。由于DBSCAN算法對高維數據定義密度很困難，所以對于二維空間中的點，可以使用歐幾里德距離來進行度量。
DBSCAN算法需要用戶輸入2個參數：一個參數是半徑（Eps），表示以給定點P為中心的圓形鄰域的范圍；另一個參數是以點P為中心的鄰域內最少點的數量（MinPts）。如果滿足：以點P為中心、半徑為Eps的鄰域內的點的個數不少于MinPts，則稱點P為核心點。
DBSCAN聚類使用到一個k-距離的概念，k-距離是指：給定數據集P={p(i); i=0,1,…n}，對于任意點P(i)，計算點P(i)到集合D的子集S={p(1), p(2), …, p(i-1), p(i+1), …, p(n)}中所有點之間的距離，距離按照從小到大的順序排序，假設排序后的距離集合為D={d(1), d(2), …, d(k-1), d(k), d(k+1), …,d(n)}，則d(k)就被稱為k-距離。也就是說，k-距離是點p(i)到所有點（除了p(i)點）之間距離第k近的距離。對待聚類集合中每個點p(i)都計算k-距離，最后得到所有點的k-距離集合E={e(1), e(2), …, e(n)}。
根據經驗計算半徑Eps：根據得到的所有點的k-距離集合E，對集合E進行升序排序后得到k-距離集合E’，需要擬合一條排序后的E’集合中k-距離的變化曲線圖，然后繪出曲線，通過觀察，將急劇發生變化的位置所對應的k-距離的值，確定為半徑Eps的值。
根據經驗計算最少點的數量MinPts：確定MinPts的大小，實際上也是確定k-距離中k的值，DBSCAN算法取k=4，則MinPts=4。
另外，如果覺得經驗值聚類的結果不滿意，可以適當調整Eps和MinPts的值，經過多次迭代計算對比，選擇最合適的參數值。可以看出，如果MinPts不變，Eps取得值過大，會導致大多數點都聚到同一個簇中，Eps過小，會導致已一個簇的分裂；如果Eps不變，MinPts的值取得過大，會導致同一個簇中點被標記為噪聲點，MinPts過小，會導致發現大量的核心點。

我們需要知道的是，DBSCAN算法，需要輸入2個參數，這兩個參數的計算都來自經驗知識。半徑Eps的計算依賴于計算k-距離，DBSCAN取k=4，也就是設置MinPts=4，然后需要根據k-距離曲線，根據經驗觀察找到合適的半徑Eps的值，下面的算法實現過程中，我們會詳細說明。對于算法的實現，首先我們概要地描述一下實現的過程：

解析樣本數據文件
計算每個點與其他所有點之間的歐幾里德距離
計算每個點的k-距離值，并對所有點的k-距離集合進行升序排序，輸出的排序后的k-距離值
將所有點的k-距離值，在Excel中用散點圖顯示k-距離變化趨勢
根據散點圖確定半徑Eps的值
根據給定MinPts=4，以及半徑Eps的值，計算所有核心點，并建立核心點與到核心點距離小于半徑Eps的點的映射
根據得到的核心點集合，以及半徑Eps的值，計算能夠連通的核心點，得到噪聲點
將能夠連通的每一組核心點，以及到核心點距離小于半徑Eps的點，都放到一起，形成一個簇
選擇不同的半徑Eps，使用DBSCAN算法聚類得到的一組簇及其噪聲點，使用散點圖對比聚類效果

然后，再詳細描述聚類過程的具體實現。

計算歐幾里德距離

我們使用的樣本數據，來自一組經緯度坐標數據，數據文件格式類似如下所示：

116.389463     39.87194
116.389463     39.874577
116.312984     39.887419
116.382798     39.853576
116.496648     39.872999
116.436246     39.911165
116.622074     40.061037
116.599267     40.062485
116.441824     39.940168
116.599267     40.062485
116.402096     39.942057
116.37319     39.93428
116.327812     39.899396
116.374739     39.898751
116.287195     39.959335
116.513574     39.878222
116.474355     39.962825
116.400651     40.008559
... ...

我們需要做的首先就是，解析樣本數據文件，將其轉換成我們需要的表示形式，我們定義了Point2D類，代碼如下所示：

package org.shirdrn.dm.clustering.common;
public class Point2D {
    protected final Double x;
    protected final Double y;
    public Point2D(Double x, Double y) {
        super();
        this.x = x;
        this.y = y;
    }
    @Override
    public int hashCode() {
        return 31 * x.hashCode() + 31 * y.hashCode();
    }
    @Override
    public boolean equals(Object obj) {
        Point2D other = (Point2D) obj;
        return this.x.doubleValue() == other.x.doubleValue() && this.y.doubleValue() == other.y.doubleValue();
    }
    public Double getX() {
        return x;
    }
    public Double getY() {
        return y;
    }
    @Override
    public String toString() {
        return "(" + x + ", " + y + ")";
    }
}

我們可以將解析后的點的對象放到一個List<Point2D> allPoints集合里面，以便后續使用時迭代集合。在計算兩點之間的歐幾里德距離時，需要迭代前面生成的Point2D的集合，計算歐幾里德距離，實現方法如下所示：

public static double euclideanDistance(Point2D p1, Point2D p2) {
    double sum = 0.0;
    double diffX = p1.getX() - p2.getX();
    double diffY = p1.getY() - p2.getY();
    sum += diffX * diffX + diffY * diffY;
    return Math.sqrt(sum);
}

如果需要，可以將計算的所有點之間的距離緩存，因為計算k-距離需要多次訪問點的歐幾里德距離，比如，可以使用Guava庫中的Cache工具：

Cache<Set<Point2D>, Double> distanceCache =
     CacheBuilder.newBuilder().maximumSize(Integer.MAX_VALUE).build();

上面代碼中，設置緩存容納足夠多（Integer.MAX_VALUE）的對象，將計算出的全部的歐幾里德距離放在內存中，便于后續迭代時重用。

計算k-距離

每個點都要計算k-距離，在計算一個點的k-距離的時候，首先要計算該點到其他所有點的歐幾里德距離，按照距離升序排序后，選擇第k小的距離作為k-距離的值，實現代碼如下所示：

Task task = q.poll(); // 從隊列q中取出一個Task，就是計算一個點的k-距離的任務
KPoint2D p1 = (KPoint2D) task.p;
final TreeSet<Double> sortedDistances = Sets.newTreeSet(new Comparator<Double>() { // 創建一個降序排序TreeSet
    @Override
    public int compare(Double o1, Double o2) {
        double diff = o1 - o2;
        if(diff > 0) {
            return -1;
        }
        if(diff < 0) {
            return 1;
        }
        return 0;
    }
});
for (int i = 0; i < allPoints.size(); i++) { // 計算點p1與allPoints集合中每個點的k-距離
    if(task.pos != i) { // 點p1與它自己的歐幾里德距離沒必要計算
        KPoint2D p2 = (KPoint2D) allPoints.get(i);
        Set<Point2D> set = Sets.newHashSet((Point2D) p1, (Point2D) p2);
        Double distance = distanceCache.getIfPresent(set); // 從緩存中取出歐幾里德距離（可能不存在）
        if(distance == null) {
            distance = MetricUtils.euclideanDistance(p1, p2);
            distanceCache.put(set, distance); // 不存在則加入緩存中
        }
        if(!sortedDistances.contains(distance)) {
            sortedDistances.add(distance);
        }
        if(sortedDistances.size() > k) { // TreeSet中只最多保留k個歐幾里德距離值
            Iterator<Double> iter = sortedDistances.iterator();
            iter.next();
            // remove (k+1)th minimum distance
            iter.remove(); // 將k+1個距離值中最大的刪除，剩余k個是最小的
        }
    }
}
// collect k-distance
p1.kDistance = sortedDistances.iterator().next();  // 此時，TreeSet中最大的，就是第k最小的距離

上述代碼中，KPoint2D類是Point2D的子類，不過比基類Point2D多了一個k-距離的屬性，代碼如下所示：

private class KPoint2D extends Point2D {

     private Double kDistance = 0.0;

     public KPoint2D(Double x, Double y) {
          super(x, y);
     }

     @Override
     public int hashCode() {
          return super.hashCode();
     }

     @Override
     public boolean equals(Object obj) {
          return super.equals(obj);
     }

}

代碼比較容易，可以查看相關注釋信息。

繪制k-距離曲線，尋找半徑Eps

x軸坐標點我們直接使用遞增的自然數序列，每個點對應一個自然數，y軸就是所有點的k-距離的大小，我們在代碼中可以進行處理，實現如下：

for(int i=0; i<allPoints.size(); i++) {
                    KPoint2D kp = (KPoint2D) allPoints.get(i);
                    System.out.println(i + "\t" + kp.kDistance);
               }

最終生成的數據，輸出格式類似如下：

0     6.795895820371257E-4
1     8.305064719800753E-4
2     8.692537028985709E-4
3     8.81717074805001E-4
4     9.38043175973106E-4
5     0.0010181414440047032
6     0.0011109837982601679
7     0.0011109837982601679
8     0.0011414013316968501
9     0.0011533646431092647
10     0.0011540277293025107
11     0.0011712783614491256
12     0.001171973122556046
13     0.001171973122556046
14     0.0012320292204251713
15     0.0012371273176228466
... ...

我們把輸出數據復制到Excel表格中，使用上述數據生成散點圖，基于x坐標取了4個不同的范圍，觀察曲線的變化情況，0~2600、0~2000、2001~2630、0~2500各個x坐標范圍內的點，對應的散點圖分別如下所示：

通過上圖可以看出：

左上圖（0~2600）：由于x=2500之后店的k-距離變化太快（可以忽略），導致前面點的k-距離的變化趨勢無法觀察出來。
右上圖（0~2000）：去掉尾部的一些點的k-距離，可以看出曲線的變化趨勢。
左下圖（2001~2630）：x坐標軸后半部分的距離的變化趨勢。
右下圖（0~2500）：去掉尾部一些k-距離點，展示大部分k-距離點的變化趨勢。

綜合上面4個圖，可以選擇得到半徑Eps的范圍大致在0.002~0.006之間，已知MinPts=4，具體我們可以選擇下面3個k-距離的值作為半徑Eps：

0.0025094814205335555
0.004417483559674606
0.006147849217403014

通過下一步進行聚類，看一下使用選擇的Eps的這幾個值進行聚類的效果對比。另外，對半徑Eps=8也進行聚類，主要是為了看一下半徑變化對聚類效果的影響。

計算核心點

聚類過程需要知道半徑Eps，半徑已知，首先需要計算有哪些核心點，給定點，如果以該點為中心半徑為Eps的鄰域內的其他點的數量大于等于MinPts，則該點就為核心點。計算核心點的實現代碼如下所示：

Point2D p1 = taskQueue.poll();
if(p1 != null) {
    ++processedPoints;
    Set<Point2D> set = Sets.newHashSet();
    Iterator<Point2D> iter = epsEstimator.allPointIterator();
    while(iter.hasNext()) {
        Point2D p2 = iter.next();
        if(!p2.equals(p1)) {
            double distance = epsEstimator.getDistance(Sets.newHashSet(p1, p2));
            // collect a point belonging to the point p1
            if(distance <= eps) { // 收集每個點與其鄰域內的點之間距離小于等于Eps的點
                set.add(p2);
            }
        }
    }
    // decide whether p1 is core point
    if(set.size() >= minPts) { // 計算收集到的鄰域內的點的個數，小于等于MinPts，則加入到映射表Map<Point2D, Set<Point2D>> corePointScopeSet中
        corePointScopeSet.put(p1, set);
        LOG.debug("Decide core point: point" + p1 + ", set=" + set);
    } else {
        // here, perhaps a point was wrongly put into noise point set
        // afterwards we should remedy noise point set
        if(!noisePoints.contains(p1)) { // 這里，會把一些點錯誤地加入到噪聲點集合noisePoints中，后面會進行修正
            noisePoints.add(p1);
        }
    }
}

那些被錯誤放入噪聲點集合的點，需要在計算核心點完成之后，遍歷噪聲點集合，與核心點集合（及其對應的映射點集合）進行比對，代碼如下所示：

// process noise points
Iterator<Point2D> iter = noisePoints.iterator();
while(iter.hasNext()) {
    Point2D np = iter.next();
    if(corePointScopeSet.containsKey(np)) {
        iter.remove();
    } else {
        for(Set<Point2D> set : corePointScopeSet.values()) {
            if(set.contains(np)) {
                iter.remove();
                break;
            }
        }
    }
}

這樣，有些非噪聲點就從噪聲點集合中被排除了。

連通核心點生成簇

核心點能夠連通（有些書籍中稱為：“密度可達”），它們構成的以Eps長度為半徑的圓形鄰域相互連接或重疊，這些連通的核心點及其所處的鄰域內的全部點構成一個簇。假設MinPts=4，則連通的核心點示例，如下圖所示：

計算連通的核心點的思路是，基于廣度遍歷與深度遍歷集合的方式：從核心點集合S中取出一個點p，計算點p與S集合中每個點（除了p點）是否連通，可能會得到一個連通核心點的集合C1，然后從集合S中刪除點p和C1集合中的點，得到核心點集合S1；再從S1中取出一個點p1，計算p1與核心點集合S1集中每個點（除了p1點）是否連通，可能得到一個連通核心點集合C2，再從集合S1中刪除點p1和C2集合中所有點，得到核心點集合S2，……最后得到p、p1、p2、……，以及C1、C2、……就構成一個簇的核心點。最終將核心點集合S中的點都遍歷完成，得到所有的簇。具體代碼實現，如下所示：

// join connected core points
LOG.info("Joining connected points ...");
Set<Point2D> corePoints = Sets.newHashSet(corePointScopeSet.keySet());
while(true) {
    Set<Point2D> set = Sets.newHashSet();
    Iterator<Point2D> iter = corePoints.iterator();
    if(iter.hasNext()) {
        Point2D p = iter.next();
        iter.remove();
        Set<Point2D> connectedPoints = joinConnectedCorePoints(p, corePoints);
        set.addAll(connectedPoints);
        while(!connectedPoints.isEmpty()) {
            connectedPoints = joinConnectedCorePoints(connectedPoints, corePoints);
            set.addAll(connectedPoints);
        }
        clusteredPoints.put(p, set);
    } else {
        break;
    }
}

上面調用了重載的兩個方法joinConnectedCorePoints，分別根據參數不同計算連通的核心點集合：計算核心點集合中一個點，與該核心點集合中其它的所有核心點是否連通，調用如下方法：

private Set<Point2D> joinConnectedCorePoints(Point2D p1, Set<Point2D> leftCorePoints) {
    Set<Point2D> set = Sets.newHashSet();
    for(Point2D p2 : leftCorePoints) {
        double distance = epsEstimator.getDistance(Sets.newHashSet(p1, p2));
        if(distance <= eps) {
            // join 2 core points to the same cluster
            set.add(p2);
        }
    }
    // remove connected points
    leftCorePoints.removeAll(set); // 刪除已經確定為與p1連通的核心點
    return set;
}

還有一個方法是，上面第一個參數變為一個集合，它調用上面的方式處理每一個點，方法代碼實現如下所示：

private Set<Point2D> joinConnectedCorePoints(Set<Point2D> connectedPoints, Set<Point2D> leftCorePoints) {
          Set<Point2D> set = Sets.newHashSet();
          for(Point2D p1 : connectedPoints) {
               set.addAll(joinConnectedCorePoints(p1, leftCorePoints));
          }
          return set;
     }

最后，集合clusteredPoints存儲的就是聚類后的核心點的信息，再根據核心點到其鄰域內半徑小于等于Eps的點的集合的映射，就能將所有的點生成聚類了。

聚類效果比較

選擇不同的半徑Eps進行聚類分析，聚類的結果也不相同，有些情況下差異很大，有些情況差異較小。比如，在MinPts=4的情況下，如何繪制k-距離曲線，已經在前面詳細說明了處理過程，我們根據k-距離趨勢增長曲線，選擇了一組半徑Eps的值，執行DBSCAN聚類算法，下面我們比較在MinPts=4和MinPts=8的情況下，再分別選擇不同的半徑Eps，執行DBSCAN聚類算法生成簇，對分布情況進行對比。

參數：MinPts=4

選擇半徑Eps的值分別為如下3個觀察值：

0.0025094814205335555
0.004417483559674606
0.006147849217403014

最終得到的聚類效果圖在下面可以看到，其中，左側為沒有噪聲點的情況各個簇的分布情況，右側是將噪聲點全部加入到圖上顯示，圖中圖例中的9999表示噪聲點，其他的都是聚類生成的簇，如下圖所示：

通過上圖可以看出，半徑Eps選擇的較小時，會產生大量的噪聲點，其實我們想一下，半徑小了，自然落在某個點的鄰域內的點減少的可能性就增加了，導致很多點可能就無法成為核心點，自然也就無法成為某個簇的點，而且很生成的簇包含的點的數量都比較少，某些本來很接近的點應該可以成為同一個簇，但是被分裂了。

當半徑比較大一些時，生成的一個簇包含了比較多的點，而且這個簇的形狀中間可能出現一些“空洞”，因為點之間的距離大一些也能滿足成為核心點的要求，所以這些點聚到一簇中可能確實比較合理。

參數：MinPts=8

選擇半徑Eps的值分別為如下3個觀察值：

0.004900098978598581
0.009566439044911
0.013621050253196359

使用我們實現的DBSCAN聚類算法進行分析處理，得到的聚類結果，如下圖所示：

總結

因為DBSCAN聚類算法，是基于密度的聚類算法，所以對于密度分別不均，各個簇的密度變化較大時，可能會導致一些問題：比如半徑Eps較大時，本來不屬于同一個的簇的點被聚到一個簇中；比如半徑Eps比較小時，會出現大量比較小的簇（即每個簇中含有的點不較少，但是這些點組成的小簇密度確實很大），同時出現了大量的點不滿足成為核心點的要求，MinPts越大越容易出現這種情況。DBSCAN聚類算法的思想非常明確易懂，雖然需要用戶輸入2個參數，但是正是輸入參數的靈活性，我們可以根據自己實際應用的需要，適當調整參數值，在特定應用場景下進行聚類分析，得到合適的簇劃分。通過上面選擇不同參數進行聚類的結果對比，如果希望局部比較密集的點都能夠生成簇，那么在固定MinPts的值的條件下，半徑Eps通過調整變小，可能達到這一需要，產生的簇的數量比較多，但是同時也產生了大量分散的噪聲點，實際上應該可以進行二次聚類，將噪聲點也通過合適的方式歸到指定的簇中；如果希望生成全局比較大的簇，可以適當調整半徑Eps變大，生成的簇的數量較少，噪聲點的數量也相對較少。

原文 http://shiyanjun.cn/archives/1288.html

本文由用戶 jopen 自行上傳分享，僅供網友學習交流。所有權歸原作者，若您的權利被侵害，請聯系管理員。

轉載本站原創文章，請注明出處，并保留原始鏈接、圖片水印。

本站是一個以用戶分享為主的開源技術平臺，歡迎各類分享！

本文地址：http://www.baiduhome.net/lib/view/open1447581496902.html

算法