二分k-均值演算法步驟：

首先將所有點作為一個屬，然後將該簇-分為二，之 後選擇其中-個簇進續進行劃分，選擇哪一個簇進行劃取決於對其劃分是否可以最大程度降低SSE的值，上述基於SSE的別分過程不斷重複，直到得到使用者指定的屬數目為止，
將所有點看成一個簇
當簇數目小於k時
      對於每一個簇:
          計算總誤差
           在給定的簇上面進行K-均值聚類(k=2)計算將該簇一分為二後的總誤差
      選擇使得誤差最小的那個簇進行劃分操作
另一種做法選擇S正最大的簇進行劃分，直到簇數目達到使用者指定的數目為止。

Python實現

# -*- coding: utf-8 -*-
"""
Created on Thu Aug 16 11:09:04 2018

@author: wjw
"""

from numpy import  *

def loadDataSet(fileName):
    dataMat=[]
    fr=open(fileName)
    for line in fr.readlines():
        arr=line.strip().split('\t')
        #map (回撥函式，列表)-》迴圈列表中的每個值，呼叫迴歸函式得到結果，存到map -》 list
        l=list(map(float,arr))#float 強制型別轉換
        dataMat.append(l)
    return dataMat
    
dataMat=loadDataSet('dataset/testSet.txt')
#print(dataMat[0:5])
    
#距離 度量方式：
#1.歐氏距離
def disEuclid(vecA,vecB):
    return sqrt(sum(power((vecA-vecB),2)))

#初始質心選取
def randCenter(dataSet,k):
    n=shape(dataSet)[1]#n=..列
    centers=mat(zeros((k,n)))
    #print(centers)
    for j in range(n): #每迴圈一次，產生的一個特徵值
        #當前第j列的最小值，與最大值，求範圍
        minJ=min(dataSet[:,j])
        #print('minJ:',minJ)
        maxJ=max(dataSet[:,j])
        rangeJ=float(maxJ-minJ)
        #生成隨機數 ，k行（代表k箇中心點）#0-1
        centers[:,j]=minJ+rangeJ*random.rand(k,1)#random.rand(k,1)->k行一列的資料 randej標量
        print(centers[:,j])
    return centers
    
def myKMeans(dataSet,k,disMea=disEuclid,initCenter=randCenter):
    '''
    
    '''
    m=shape(dataSet)[0]
    #這個zeros是一個m行2列的資料（記錄這個點所屬的簇的索引，，記錄這個點到其質心的距離）
    clusterAssment=mat(zeros((m,2)))
    
    #生成隨機質心
    centers=initCenter(dataSet,k)
    
    clusterChanged=True
    while clusterChanged:
        clusterChanged=False
        #迴圈每個點，計算他與每個質心的位置
        for i in range (m):
            #這個點到某質心的最小距離及質心所在的索引
            minDist=inf #某個點距離質心點的最小距離 inf 無窮大
            minIndex=-1#有最小距離的質心的索引——知道這個簇了
            for j in range(k):
                distJI=disMea(centers[j,:],dataSet[i,:])#算距離
                #print('===:',centers[j,:])
                #print(']][][[][]]',dataSet[i,:])
                #print('================================================================')
                if distJI < minDist:
                    minDist=distJI
                    minIndex=j
                #for I in   range(m)迴圈完，表明已經找到了第i個數據點所屬的簇，且計算出了距離         
            if clusterAssment[i,0]!=minIndex:
                clusterChanged=True
                #更新這個點到質心的索引及誤差
            clusterAssment[i,:]=minIndex,minDist**2
            
            #============以上的迴圈個更新每個點的簇======================================================================
            
        #遍歷所有的簇，重新找質心
        for cent in range(k):
            flag=clusterAssment[:,0].A==cent #查詢這個 cent 簇所有的點
            #print('flag:',flag)
            pointsInCluster=dataSet[nonzero(flag)[0]]#第cent個簇所有點
            #print('pointsInCluster:',pointsInCluster)
            centers[cent,:]=mean(pointsInCluster,axis=0)#對於這個簇中每個點的列取均值，更新中心點centers[cent]
    return centers,clusterAssment
    
    
dataMat=mat(loadDataSet('dataset/testSet.txt'))
centers,clusterAssment=myKMeans(dataMat,4)
print(centers)
print(clusterAssment)#第一個列為簇的編號，第二列是當前點到這個簇的質心的距離
    
def biKmeans(dataSet,k,distMea=disEuclid):
    #1.將所有點看成一個簇
    #2.取出這個簇的中心點
    #3.計算這個簇中每個點到中心點的距離
    #4.這個zeros是一個m行2列的資料（記錄這個點所屬的簇的索引，，記錄這個點到其質心的距離）
    #clusterAssment=mat(zeros((m,2)))
    #存所有的中心點
    
    m=shape(dataSet)[0]
    #這個zeros是一個m行2列的資料（記錄這個點所屬的簇的索引，，記錄這個點到其質心的距離）
    clusterAssment=mat(zeros((m,2)))
    center0=mean(dataSet,axis=0).tolist()[0]
    print('第一個中心點：',center0)
    centList=[center0]
    for j in range(m):
        clusterAssment[j,1]=distMea(mat(center0),dataSet[j,:])**2
    #迴圈來產生質心
    
    while(len(centList)<k):#當簇數目小於k時
        lowestSSE=inf #初始化最小 誤差平方和
        #迴圈每個簇
        for i in range (len(centList)):
            #到dataset中篩選出屬於第i個簇的資料樣本
            pointsInCluster=dataSet[nonzero(clusterAssment[:,0].A==i)[0],:]
            #對這個簇中的樣本進行一次k=2的聚類
            centroidMat,splitClustAss=myKMeans(pointsInCluster,2,distMea)
            sseSplit=sum(splitClustAss[:,1])
            #剩餘的資料集的誤差
            #sseSplit=sum(splitClustAss[:,1])-sseSplit
            sseNotSplit=sum(clusterAssment[nonzero(clusterAssment[:,0].A!=i)[0],1])
            #總誤差
            totalSplit=sseSplit+sseNotSplit
            if totalSplit<lowestSSE:
                bestCentToSplit=i#最好的質心索引
                bestNewCent=centroidMat
                bestClustAss=splitClustAss.copy()
                lowestSSE=totalSplit
        #2分k聚類返回係數0或1，需要把1換成當前簇數目，以免造成重複
        bestClustAss[ nonzero(bestClustAss[:,0].A==1)[0],0 ]=len(centList)
        #返回的是第幾行第幾行
        #把0換成別切分的簇，或者與上面的交換賦值也可以
        bestClustAss[ nonzero(bestClustAss[:,0].A==0)[0],0]=bestCentToSplit
        #print('===:',bestCentToSplit)
        #將centlist指定位置上的質心換成分割後的質心
        centList[bestCentToSplit]=bestNewCent[0,:].tolist()[0]
        #將另一個質心新增上去
        centList.append(bestNewCent[1,:].tolist()[0])
        #將劃分後的新質點及點分佈賦值給結果矩陣
        clusterAssment[nonzero(clusterAssment[:,0].A==bestCentToSplit)[0],:]=bestClustAss
    return mat(centList),clusterAssment

datMat=mat(loadDataSet('dataset/testSet.txt'))
centList,clusterAssment=biKmeans(dataMat,4)
print(centList)
        
 

注：biKmeans 函式就是二分k-均值演算法 

可藉助spyder的debug模式進行除錯，幫助理解程式碼，如果spyder版本太低可以如下命令進行升級

conda update spyder