原創： 牛超   2017-02-21   Mail：[email protected]

引言：前一篇文章開始不再介紹簡單演算法，而是轉到資料探勘之旅。感謝CSDN將我前一篇機器學習C4.5決策樹演算法的博文推送到了首頁，也非常榮幸能夠得到雲和恩墨的蓋老師的肯定。個人會繼續以精品博文分享我的觀點，也希望結交更多的資料庫、資料分析領域的朋友。

從前面幾篇可以看到，SQL做為一種程式語言，能夠滿足各類資料處理的需要，關鍵就在於演算法與思維方式。個人經常調侃SQL思考問題比大部分流行的開發語言多一個維度，因為SQL主要是二維思考（集合）、區別於一維（資料結構）的思維方式。對於ORACLE，通過以SQL（相對巨集觀）為主體、PLSQL（微觀）為輔助，注入演算法（靈魂），貫徹效能優化（章程），資料的價值能夠充分有效發揮。

聚類問題，就是給定一個元素集合D，其中每個元素具有n個可觀察屬性，使用某種演算法將D劃分成k個子集，要求每個子集內部的元素之間相異度儘可能低，而不同子集的元素相異度儘可能高。其中每個子集叫做一個簇。本文將介紹聚類的經典演算法K均值聚類演算法，即K-MEANS，是一種觀察類學習，通過以元素間的相異度迭代地劃分簇並重新定位質心點重新聚類來達成的演算法，找了如下的圖以便加深理解。

演算法本身較容易理解，但為了定量說明，還是先祭出幾個數學公式

 用來定義兩個元素，各自具有n個度量屬性。

 用來標量X與Y的相異度（尤拉距離公式），本篇採用該公式。

 曼哈頓距離，即街區非直線段距離，很容易理解。也可以用來標量元素間相異度。

 標量規格化，為了平衡各個屬性因取值單位不同對距離的影響而按比例對映到相同的取值區間。通常將各個屬性均對映到[0,1]區間。

接著，我們來看看本次要用SQL實現的k-means演算法示例：以2016年的GDP統計資料給中國城市分級：

目前這份是經個人簡單加工的2016年真實國家統計資料，一共有100個城市的記錄，共5個度量（人口、面積、年GDP、人均GDP、單位面積GDP)。

我們總是說不相信眼淚的北上廣是一線城市，那麼該如何給城市分成四線是個很有趣的問題，看了一下國家統計局的劃分標準（此處省略1K字），但做為個人比較關心的是每個城市發達與市民富裕程度，於是簡單地拍腦袋定義瞭如下模型：

X ={人口，人口密度，年GDP，人均GDP，單位面積GDP}，共5個維度，權重都是1;

每個度量屬性a的標量規格化區間定義為[0,1];

聚合點相異度需要以質心點為基準計算，初始質心點取GDP名次的4等分點，簡化為MOD(名次,21)=0;

元素相異度直接使用尤拉距離公式，即

準備工作：

1.定義表象性的業務表，即城市GDP資料表

--城市GDP
CREATE TABLE CITY_GDP_T
(
CITY VARCHAR2(100),--城市名
POPULATION NUMBER ,--人口
AREA NUMBER,--面積
GDP_YEAR NUMBER,--年GDP
GDP_PER_CAPITA NUMBER ,--人均GDP
GDP_PER_AREA NUMBER --單位面積GDP
);

3.先預演一下質心點經過一次聚類後重新被選擇的演算法程式，其中第一代初始質心點根據GDP的分段城市的元素屬性，TA1，再根據TA1的聚類點用算術平均法計算得到第二代質心點，SQL如下：

WITH TA AS --整理度量值
 (SELECT T.CITY,
         AREA,
         ROW_NUMBER() OVER(ORDER BY GDP_YEAR DESC) RNUM,
         --人口
         POPULATION,
         MIN(POPULATION) OVER() MIN_POPULATION,
         MAX(POPULATION) OVER() MAX_POPULATION,
         --人口密度
         ROUND(POPULATION / AREA, 2) POLA_PER_AREA,
         MIN(ROUND(POPULATION / AREA, 2)) OVER() MIN_POLA_PER_AREA,
         MAX(ROUND(POPULATION / AREA, 2)) OVER() MAX_POLA_PER_AREA,
         --年GDP
         GDP_YEAR,
         MIN(GDP_YEAR) OVER() MIN_GDP_YEAR,
         MAX(GDP_YEAR) OVER() MAX_GDP_YEAR,
         --人均GDP
         GDP_PER_CAPITA,
         MIN(GDP_PER_CAPITA) OVER() MIN_GDP_PER_CAPITA,
         MAX(GDP_PER_CAPITA) OVER() MAX_GDP_PER_CAPITA,
         --單位面積GDP
         GDP_PER_AREA,
         MIN(GDP_PER_AREA) OVER() MIN_GDP_PER_AREA,
         MAX(GDP_PER_AREA) OVER() MAX_GDP_PER_AREA
    FROM CITY_GDP_T T),
TB AS --規格化，以消除屬性值單位不同造成的影響
 (SELECT CITY,
         RNUM,
         (POPULATION - MIN_POPULATION) / (MAX_POPULATION - MIN_POPULATION) POLA_VALUE, --人口 ，規格化到[0,1]區間
         (POLA_PER_AREA - MIN_POLA_PER_AREA) /
         (MAX_POLA_PER_AREA - MIN_POLA_PER_AREA) PPA_VALUE, --人口密碼
         (GDP_YEAR - MIN_GDP_YEAR) / (MAX_GDP_YEAR - MIN_GDP_YEAR) GDP_VALUE, --GDP
         (GDP_PER_CAPITA - MIN_GDP_PER_CAPITA) /
         (MAX_GDP_PER_CAPITA - MIN_GDP_PER_CAPITA) GPC_VALUE, --人均GDP 
         (GDP_PER_AREA - MIN_GDP_PER_AREA) /
         (MAX_GDP_PER_AREA - MIN_GDP_PER_AREA) GPA_VALUE --單位面積GDP
    FROM TA
  --where mod(rnum,10) = 5
  ),
TA1 AS --第一代質心點選擇，根據GDP
 (SELECT RNUM       CLUSTERID,
         POLA_VALUE POLA_CENTROID,
         PPA_VALUE  PPA_CENTRIOID,
         GDP_VALUE  GDP_CENTROID,
         GPC_VALUE  GPC_CENTROID,
         GPA_VALUE  GPA_CENTROID
    FROM TB
   WHERE TB.RNUM > 0
     AND MOD(TB.RNUM, 21) = 0 --選質心點，利用GDP值
  ),
TC AS
 (SELECT A.CLUSTERID,
         B.CITY CTB,
         B.RNUM CTB_ROW,
         B.POLA_VALUE,
         B.GDP_VALUE,
         B.GPC_VALUE,
         B.GPA_VALUE,
         SQRT(POWER(B.POLA_VALUE - A.POLA_CENTROID, 2) +
              POWER(B.GDP_VALUE - A.GDP_CENTROID, 2) +
              POWER(B.GPC_VALUE - A.GPC_CENTROID, 2) +
              POWER(B.GPA_VALUE - A.GPA_CENTROID, 2)) DVALUE --尤拉距離
    FROM TA1 A
   CROSS JOIN TB B),
TD AS
 (SELECT TC.*,
         ROW_NUMBER() OVER(PARTITION BY CTB_ROW ORDER BY DVALUE) SUB_RNUM
    FROM TC),
TE AS --聚類選擇，各元素取相異度最低的質心點
 (SELECT * FROM TD WHERE SUB_RNUM = 1),
TA2 AS --第二代質心點選擇，根據各屬性算術平均值
 (SELECT TE.CLUSTERID CLUSTERID,
         AVG(POLA_VALUE) POLA_CENTROID,
         AVG(GDP_VALUE) GDP_CENTROID,
         AVG(GPC_VALUE) GPC_CENTROID,
         AVG(GPA_VALUE) GPA_CENTROID
    FROM TE
   GROUP BY CLUSTERID)
SELECT * FROM TA2

可以在集合TA1後面做一個SELECT看一下第一代的質心點，如下圖：

執行SQL後看一下第二代的質心點，發現維值都發生了變化，說明質心點還是不穩定的，需要迭代地尋找下去。

通過上面的預演，腦子也做了預熱。找到規律之後，霍然思路全部連通，K-MEANS聚類問題的關鍵就在於遞迴地尋找最穩定的質心點集合。為了保持演算法的通用性，抽象出了如下8個維度的聚類訓練集，同時定義了批次ID與初始質心標識（0，1）：

--K-MEANS聚類訓練集
CREATE TABLE DM_KMEANS_LEANING_T
(
EID NUMBER ,--元素ID
ELABEL VARCHAR2(100),
BATCH_ID NUMBER ,
C1 NUMBER ,--維度1質心座標 規格化後
C2 NUMBER ,
C3 NUMBER ,
C4 NUMBER ,
C5 NUMBER ,
C6 NUMBER ,
C7 NUMBER ,
C8 NUMBER ,
CENTRIOD_FLAG NUMBER(1)--初始質心 1/0
);
-- Create/Recreate indexes 
CREATE INDEX DM_KMEANS_LEANING_N1 ON DM_KMEANS_LEANING_T (BATCH_ID);

CREATE OR REPLACE TYPE DM_KMEANS_CENTROID_OBJ IS OBJECT
(
  CLUSTER_ID NUMBER, --簇ID
  C1         NUMBER, --維度1質心座標 規格化後
  C2         NUMBER,
  C3         NUMBER,
  C4         NUMBER,
  C5         NUMBER,
  C6         NUMBER,
  C7         NUMBER,
  C8         NUMBER,
  GENERATION NUMBER  --子代數
); 
CREATE OR REPLACE TYPE DM_KMEANS_CENTROID_TAB IS TABLE OF DM_KMEANS_CENTROID_OBJ;

最重要的演算法實現，考慮到ORACLE自定義函式本身是可遞迴的，我們便來實現這麼一個質心點選擇的遞迴函式，相異度選用尤拉距離公式，為了防止迭代次數過多影響效能，可以指定引數P_MAX_GENERATIOIN NUMBER，即最大遞迴子代數，具體指令碼如下：

CREATE OR REPLACE FUNCTION FUN_DM_KMEANS_CENTROID(P_BATCH_ID NUMBER, --批次ID
                                                  P_MAX_GENERATIOIN NUMBER DEFAULT NULL ,--最大遞迴子代數
                                                  P_CENTRIOD_LIST DM_KMEANS_CENTROID_TAB DEFAULT NULL --質心點
                                                 )
  RETURN DM_KMEANS_CENTROID_TAB AS
  V_PARENT_LIST  DM_KMEANS_CENTROID_TAB; --本次迭代質心點 D
  V_CURRENT_LIST DM_KMEANS_CENTROID_TAB; --質心點計算結果 D'
  V_CNT          NUMBER;-- D'- D元素數
BEGIN
  --取質心點 D
  IF P_CENTRIOD_LIST IS NOT NULL THEN
    V_PARENT_LIST := P_CENTRIOD_LIST;
  ELSE
    SELECT DM_KMEANS_CENTROID_OBJ(EID, C1, C2, C3, C4, C5, C6, C7, C8, 1) BULK COLLECT
      INTO V_PARENT_LIST
      FROM DM_KMEANS_LEANING_T
     WHERE BATCH_ID = P_BATCH_ID
       AND CENTRIOD_FLAG = 1;
    IF P_MAX_GENERATIOIN = 1 THEN
      RETURN V_PARENT_LIST ;
    END IF ;
  END IF;
  --計算修正後的質心點 D'
  WITH TC AS --計算各點相對質心點相異度（笛卡爾集）
   (SELECT 
     A.CLUSTER_ID,
     A.GENERATION,
     B.EID CTB_ROW,
     B.C1,
     B.C2,
     B.C3,
     B.C4,
     B.C5,
     B.C6,
     B.C7,
     B.C8,
     SQRT(NVL(POWER(B.C1 - A.C1, 2), 0) + NVL(POWER(B.C2 - A.C2, 2), 0) +
          NVL(POWER(B.C3 - A.C3, 2), 0) + NVL(POWER(B.C4 - A.C4, 2), 0) +
          NVL(POWER(B.C5 - A.C5, 2), 0) + NVL(POWER(B.C6 - A.C6, 2), 0) +
          NVL(POWER(B.C7 - A.C7, 2), 0) + NVL(POWER(B.C8 - A.C8, 2), 0)) DVALUE --歐氏距離
      FROM TABLE(V_PARENT_LIST) A
      JOIN DM_KMEANS_LEANING_T B
        ON B.BATCH_ID = P_BATCH_ID),
  TD AS--相異度排名
   (SELECT TC.*,
           ROW_NUMBER() OVER(PARTITION BY CTB_ROW ORDER BY DVALUE) SUB_RNUM
      FROM TC),
  TE AS--聚類選擇，各元素取相異度最低的質心點
   (SELECT * FROM TD WHERE SUB_RNUM = 1)
  SELECT DM_KMEANS_CENTROID_OBJ(TE.CLUSTER_ID,
                                AVG(C1), --平均值求新質心 C1維
                                AVG(C2),
                                AVG(C3),
                                AVG(C4),
                                AVG(C5),
                                AVG(C6),
                                AVG(C7),
                                AVG(C8),
                                GENERATION + 1 --子代+1
                                ) BULK COLLECT
    INTO V_CURRENT_LIST --生成D'
    FROM TE
   GROUP BY CLUSTER_ID, GENERATION;
   
   --判斷最大子代限制條件
   IF P_MAX_GENERATIOIN IS NOT NULL AND V_CURRENT_LIST(1).GENERATION = P_MAX_GENERATIOIN THEN
     RETURN V_CURRENT_LIST;
   ELSE
    SELECT COUNT(1)
      INTO V_CNT
      FROM (SELECT C1, C2, C3, C4, C5, C6, C7, C8
             FROM TABLE(V_CURRENT_LIST)
            MINUS
            SELECT C1, C2, C3, C4, C5, C6, C7, C8
             FROM TABLE(V_PARENT_LIST) );
    IF V_CNT = 0 THEN
      RETURN V_CURRENT_LIST; --質心點集穩定後返回當前子代
    ELSE
      RETURN FUN_DM_KMEANS_CENTROID(P_BATCH_ID,P_MAX_GENERATIOIN, V_CURRENT_LIST); --遞迴計算質心點集
    END IF;
   END IF ;
END;

首先我們要把業務資料轉換載入到訓練集中，這是個簡單的ETL過程，將城市GDP表資料經過抽取、維值[0,1]規格化轉換、分配批次號3後最終載入到目標K-MEAN訓練集：

INSERT INTO DM_KMEANS_LEANING_T
  (EID, ELABEL, BATCH_ID, C1, C2,C3,C4,C5,CENTRIOD_FLAG)
  WITH TA AS --整理度量值
   (SELECT T.CITY,
           AREA,
           ROW_NUMBER() OVER(ORDER BY GDP_YEAR DESC) RNUM,
           --人口
           POPULATION,
           MIN(POPULATION) OVER () MIN_POPULATION,
           MAX(POPULATION) OVER () MAX_POPULATION,
           --人口密度
           ROUND(POPULATION / AREA, 2) POLA_PER_AREA,
           MIN(ROUND(POPULATION / AREA, 2)) OVER() MIN_POLA_PER_AREA,
           MAX(ROUND(POPULATION / AREA, 2)) OVER() MAX_POLA_PER_AREA,
           --年GDP
           GDP_YEAR,
           MIN(GDP_YEAR) OVER() MIN_GDP_YEAR,
           MAX(GDP_YEAR) OVER() MAX_GDP_YEAR,
           --人均GDP
           GDP_PER_CAPITA,
           MIN(GDP_PER_CAPITA) OVER() MIN_GDP_PER_CAPITA,
           MAX(GDP_PER_CAPITA) OVER() MAX_GDP_PER_CAPITA,
           --單位面積GDP
           GDP_PER_AREA,
           MIN(GDP_PER_AREA) OVER() MIN_GDP_PER_AREA,
           MAX(GDP_PER_AREA) OVER() MAX_GDP_PER_AREA
      FROM CITY_GDP_T T),
  TB AS --規格化，以消除屬性值單位不同造成的影響
   (SELECT CITY,
           RNUM,
           (POPULATION - MIN_POPULATION ) / (MAX_POPULATION - MIN_POPULATION) POLA_VALUE ,--人口 ，規格化到[0,1]區間
           (POLA_PER_AREA - MIN_POLA_PER_AREA ) / (MAX_POLA_PER_AREA - MIN_POLA_PER_AREA) PPA_VALUE ,--人口密碼
           (GDP_YEAR - MIN_GDP_YEAR ) / (MAX_GDP_YEAR - MIN_GDP_YEAR) GDP_VALUE ,--GDP
           (GDP_PER_CAPITA - MIN_GDP_PER_CAPITA) / (MAX_GDP_PER_CAPITA - MIN_GDP_PER_CAPITA) GPC_VALUE,--人均GDP 
           (GDP_PER_AREA - MIN_GDP_PER_AREA) / (MAX_GDP_PER_AREA - MIN_GDP_PER_AREA) GPA_VALUE --單位面積GDP
      FROM TA
    ), -- ,
  TR1 AS
   (SELECT RNUM CLUSTERID,
           CITY,
           3 BATCH_ID,
           POLA_VALUE POLA_CENTRIOD,
           PPA_VALUE PPA_CENTRIOD,
           GDP_VALUE GDP_CENTRIOD,
           GPC_VALUE GPC_CENTROID,           
           GPA_VALUE GPA_CENTROID,
           CASE WHEN TB.RNUM > 0 AND MOD(TB.RNUM, 21) = 0 THEN 1 ELSE 0 END CENTRIOD_FLAG --選質心點，利用GDP值
      FROM TB)
  SELECT * FROM TR1 ;
  COMMIT ;

SELECT *
FROM TABLE( FUN_DM_KMEANS_CENTROID(3)) A  ;

通過計算結果，我們可以看到質心點在第10代(GENERATION)才得以穩定下來，如果想看看第5代長什麼樣怎麼辦？簡單，如下：

SELECT *
FROM TABLE( FUN_DM_KMEANS_CENTROID(3,5)) A  ;

接下來我們回到問題的初衷，如何把城市分成四線，具體哪些城市是一線，只需要下面這一個SQL，和質心點選擇函式中功能大同小異：

WITH PARAMS AS (
   SELECT 3 BATCH_ID , --批次ID
          NULL MAX_GENERATION--最大子代數 
   FROM DUAL 
   ) , TC AS
   (SELECT /*+ORDERED USE_NL(PM,A,B)*/ --HINT，引導選擇正確執行計劃
     A.CLUSTER_ID,
     A.GENERATION,
     SQRT(NVL(POWER(A.C1,2), 0) + NVL(POWER(A.C2,2), 0) + NVL(POWER(A.C3,2), 0) + NVL(POWER(A.C4,2), 0) 
     + NVL(POWER(A.C5,2), 0) + NVL(POWER(A.C6,2), 0) + NVL(POWER(A.C7,2), 0) + NVL(POWER(A.C8,2), 0) ) D_CLUSTER,--質點到原點尤拉距離,
     B.EID CTB_ROW,
     B.ELABEL ,
     B.C1,--人口
     B.C2,--人口密度
     B.C3,--年GDP
     B.C4,--人均GDP
     B.C5,--單位面積GDP
     B.C6,
     B.C7,
     B.C8,
     SQRT(NVL(POWER(B.C1,2), 0) + NVL(POWER(B.C2,2), 0) + NVL(POWER(B.C3,2), 0) + NVL(POWER(B.C4,2), 0) 
     + NVL(POWER(B.C5,2), 0) + NVL(POWER(B.C6,2), 0) + NVL(POWER(B.C7,2), 0) + NVL(POWER(B.C8,2), 0) ) D_EID,--元素點到原點尤拉距離,     
     SQRT(NVL(POWER(B.C1 - A.C1, 2), 0) + NVL(POWER(B.C2 - A.C2, 2), 0) +
          NVL(POWER(B.C3 - A.C3, 2), 0) + NVL(POWER(B.C4 - A.C4, 2), 0) +
          NVL(POWER(B.C5 - A.C5, 2), 0) + NVL(POWER(B.C6 - A.C6, 2), 0) +
          NVL(POWER(B.C7 - A.C7, 2), 0) + NVL(POWER(B.C8 - A.C8, 2), 0)) DVALUE --尤拉距離 質點到元素點尤拉距離
      FROM PARAMS PM
      CROSS JOIN TABLE( FUN_DM_KMEANS_CENTROID(PM.BATCH_ID , PM.MAX_GENERATION)) A
      JOIN DM_KMEANS_LEANING_T B
        ON B.BATCH_ID = PM.BATCH_ID),
  TD AS
   (SELECT TC.*,
           ROW_NUMBER() OVER(PARTITION BY CTB_ROW ORDER BY DVALUE) SUB_RNUM
      FROM TC)
  SELECT * FROM TD WHERE SUB_RNUM = 1
  ORDER BY D_CLUSTER DESC, D_EID DESC;

是不是和我一樣迫不及待地想看結果了，我所關心的城市到底被分到了哪一級，輸出結果:

如此便計算出了我心目中的四線城市。對結果不喜的，莫爭議，這就是個一個數字遊戲，畢竟只是堆疊出來的度量模型沒什麼權威。簡單分析一下，CUSTER_ID值的大小不能說明什麼，只是用來給簇編號確定分類的。根據CLUSTER_ID分類，可以看到北上廣深以及其他的直轄市都在最繁榮的分類中，蘇州、成都能夠擠進去說明很有實力。鄂爾多斯領跑二線。。。這個城市也很有趣。而我的家鄉煙臺只能搭上三線的邊，難免有些失落。

至此，SQL版本的K-MEANS聚類演算法已經介紹完，個人舉的例子可能沒有那麼貼切。因為對資料探勘來說，資料量太小，結果的偶然性會比較高。但麻雀雖小，卻較為完整地用SQL表述了K-MEANS聚類的思想。實現這麼個演算法，全篇沒有用到一個迴圈處理，還是那句話，資料處理，SQL為王。Thinking in SQL，處處是集合，或許，只有你想不到的。