零.廣告

        本文所有程式碼實現均可以在 DML 找到，不介意的話請大家在github裡給我點個Star

一.引入 

K近鄰演算法作為資料探勘十大經典演算法之一，其演算法思想可謂是intuitive，就是從訓練集裡找離預測點最近的K個樣本來預測分類

        因為演算法思想簡單，你可以用很多方法實現它，這時效率就是我們需要慎重考慮的事情，最簡單的自然是求出測試樣本和訓練集所有點的距離然後排序選擇前K個，這個是O(nlogn)的，而其實從N個數據找前K個數據是一個很常見的演算法題，可以用最大堆(最小堆)實現，其效率是O（nlogk）的，而最廣泛的演算法是使用kd樹來減少掃描的點，這也就是這篇文章的主要內容，本文偏實現，詳細理論教程見

二.前提：堆的實現

        堆是一種二叉樹，用一個數組儲存，對於k號元素，k*2號是其左兒子，k*2+1號是其右兒子

        而大根堆就是跟比左兒子和右兒子都大，小根堆反之。

        要滿足這個條件我們需要通過up( index )操作和down( index )維護它的結構

        當然講這個的文章實在有些多了，隨便搜一篇大家看看：點選開啟連結

        大小根堆的作用是

              a) 優先佇列：因為第一個元素是最大或者最小的元素，所以可以實現優先佇列

              b) 前K個最大(最小)值

        因為事先KD-tree+BBF 要同時用到這兩個東西，所以把它們實現在了同一個類裡，感覺程式碼略漂亮，貼出來觀賞一下：

from __future__ import division
import numpy as np
import scipy as sp
def heap_judge(a,b):
        return a>b

class Heap:
        def __init__(self,K=None,compare=heap_judge):
                '''
                        'K'                 is the parameter to restrict the length of Heap
                                                !!! when K is confirmed,the Min heap contain Max K elements
                                                                  while Max heap contain Min K elements
                        'compare'         is the compare function which return a BOOL when pass two variable
                                                default is Max heap
                '''
                self.K=K
                self.compare=compare
                self.heap=['#']
                self.counter=0
        def insert(self,a):
                #print self.heap
                if self.K!=None:
                        print a.x,'==='
                if self.K==None:
                        self.heap.append(a)
                        self.counter+=1
                        self.up(self.counter)
                else:
                        if self.counter<self.K:
                                self.heap.append(a)
                                self.counter+=1
                                self.up(self.counter)
                        else:
                                if (not self.compare(a,self.heap[1])):
                                        self.heap[1]=a
                                        self.down(1)
                return
        def up(self,index):
                if (index==1):
                        return
                '''
                print index
                for t in range(index+1):
                        if t==0:
                                continue
                        print self.heap[t].x
                print 
                '''
                if self.compare(self.heap[index],self.heap[int(index/2)]):
                        #fit the condition
                        self.heap[index],self.heap[int(index/2)]=self.heap[int(index/2)],self.heap[index]
                        self.up(int(index/2))
                return
        def down(self,index):
                if 2*index>self.counter:
                        return
                tar_index=0
                if 2*index<self.counter:
                        if self.compare(self.heap[index*2],self.heap[index*2+1]):
                                tar_index=index*2
                        else:
                                tar_index=index*2+1
                else:
                        tar_index=index*2
                if not self.compare(self.heap[index],self.heap[tar_index]):
                        self.heap[index],self.heap[tar_index]=self.heap[tar_index],self.heap[index]
                        self.down(tar_index)
                return

        def delete(self,index):
                self.heap[index],self.heap[self.counter]=self.heap[self.counter],self.heap[index]
                self.heap.pop()
                self.counter-=1
                self.down(index)
                pass

        def delete_ele(self,a):
                try:
                        t=self.heap.index(a)
                except ValueError:
                        t=None
                if t!=None:
                        self.delete(t)
                return t
 

          compare引數是比較大小的，預設是“數”的大根堆，你可以往堆裡傳任何類，只要有相適應的compare引數，比如我們KD-tree傳的就是KD-Node

三.KD-BFF的原理：

          首先從KD-Tree的建立說起：（直接貼《統計學習方法》的內容了）

          事實上從選擇哪一個feature開始切割，還可以選擇方差最大的那個引數，但是考慮到簡便，以及我們可以選擇更多的相似性度量方法，還是用《統計學習方法》裡面的選擇方式了。

          然後是KD-tree搜尋的方法：（來自《統計學習方法》，但注意這裡是最近鄰，也就是k=1的時候）

          那麼我們要K近鄰要怎麼做呢？就是用堆的第二個應用，用大根堆保持K個最小的距離，然後用根的距離(也就是其中最大的一個)來作為判斷的依據是否有更近的點不在結果中，這一點很重要！

          同時摘錄july部落格的一段讀者留言講得非常好的：

              在某一層，分割面是第ki維，分割值是kv，那麼 abs(q[ki]-kv) 就是沒有選擇的那個分支的優先順序，也就是計算的是那一維上的距離； 同時，從優先佇列裡面取節點只在某次搜尋到葉節點後才發生，計算過距離的節點不會出現在佇列的，比如1~10這10個節點，你第一次搜尋到葉節點的路徑是1-5-7，那麼1，5，7是不會出現在優先佇列的。換句話說，優先佇列裡面存的都是查詢路徑上節點對應的相反子節點，比如：搜尋左子樹，就把對應這一層的右節點存進佇列。

          大致這就是我們實現的基本思路了

四.KD-BFF的實現：

        知道原理了，並且有了堆這個工具之後我們就可以著手實現這個演算法了：（終於要貼程式碼了）

       程式碼~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~此程式碼是 dml / KNN / kd.py

from __future__ import division
import numpy as np
import scipy as sp
from operator import itemgetter
from scipy.spatial.distance import euclidean
from dml.tool import Heap
class KDNode:
	def __init__(self,x,y,l):
		self.x=x
		self.y=y
		self.l=l
		self.F=None
		self.Lc=None
		self.Rc=None
		self.distsToNode=None

class KDTree:
	def __init__(self,X,y=None,dist=euclidean):
		self.X=X
		self.k=X.shape[0] #N
		self.y=y
		self.dist=dist
		self.P=self.maketree(X,y,0)
		self.P.F=None
	def maketree(self,data,y,deep):
		if data.size==0:
			return None
		lenght = data.shape[0]
		case = data.shape[1]
		p=int((case)/2)
		l = (deep%self.k)
		#print data
		data=np.vstack((data,y))
		data=np.array(sorted(data.transpose(),key=itemgetter(l))).transpose()
		#print data
		y=data[lenght,:]
		data=data[:lenght,:]
		
		v=data[l,p]
		rP=KDNode(data[:,p],y[p],l)
		#print data[:,p],y[p],l
		if case>1:
			ldata=data[:,data[l,:]<v]
			ly=y[data[l,:]<v]
			data[l,p]=v-1
			rdata=data[:,data[l,:]>=v]
			ry=y[data[l,:]>=v]
			data[l,p]=v
			rP.Lc=self.maketree(ldata,ly,deep+1)
			if rP.Lc!=None:
				rP.Lc.F=rP
			rP.Rc=self.maketree(rdata,ry,deep+1)
			if rP.Rc!=None:
				rP.Rc.F=rP
		return rP

	def search_knn(self,P,x,k,maxiter=200):
		def pf_compare(a,b):
			return self.dist(x,a.x)<self.dist(x,b.x)
		def ans_compare(a,b):
			return self.dist(x,a.x)>self.dist(x,b.x)
		pf_seq=Heap(compare=pf_compare)
		pf_seq.insert(P)    #prior sequence
		ans=Heap(k,compare=ans_compare)  #ans sequence
		while pf_seq.counter>0:
			t=pf_seq.heap[1]
			pf_seq.delete(1)
			flag=True
			if ans.counter==k:
				now=t.F
				#print ans.heap[1].x,'========'
				if now != None:
					q=x.copy()
					q[now.l]=now.x[now.l]
					length=self.dist(q,x)
					if length>self.dist(ans.heap[1].x,x):
						flag=False
					else:
						flag=True
				else:
					flag=True
			if flag:
				tp,pf_seq,ans=self.to_leaf(t,x,pf_seq,ans)
			#print "============="
			#ans.insert(tp)
		return ans


	def to_leaf(self,P,x,pf_seq,ans):
		tp=P
		if tp!=None:
			ans.insert(tp)
			if tp.x[tp.l]>x[tp.l]:
				if tp.Rc!=None:
					pf_seq.insert(tp.Rc)
				if tp.Lc==None:
					return tp,pf_seq,ans
				else:
					return self.to_leaf(tp.Lc,x,pf_seq,ans)
			if tp.Lc!=None:
				pf_seq.insert(tp.Lc)
			if tp.Rc==None:
					return tp,pf_seq,ans
			else:
					return self.to_leaf(tp.Rc,x,pf_seq,ans)

        然後KNN就是對上面這個類的一個包裝：

        程式碼~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~此程式碼是 dml / KNN / knn.py

#coding:utf-8 
import numpy as np
import scipy as sp
from scipy.spatial.distance import cdist
from scipy.spatial.distance import euclidean
from dml.KNN.kd import KDTree

#import pylab as py
class KNNC:
	"""docstring for KNNC"""
	def __init__(self,X,K,labels=None,dist=euclidean):
		'''
			X is a N*M matrix where M is the case 
			labels is prepare for the predict.
			dist is the similarity measurement way,

			The distance function can be ‘braycurtis’, ‘canberra’, 
			‘chebyshev’, ‘cityblock’, ‘correlation’, ‘cosine’, 
			‘dice’, ‘euclidean’, ‘hamming’, ‘jaccard’, ‘kulsinski’, 
			‘mahalanobis’, 

		'''
		self.X = np.array(X)
		if labels==None:
			np.zeros((1,self.X.shape[1]))
		self.labels = np.array(labels)
		self.K = K
		self.dist = dist
		self.KDTrees=KDTree(X,labels,self.dist)

	def predict(self,x,k):
		ans=self.KDTrees.search_knn(self.KDTrees.P,x,k)
		dc={}
		maxx=0
		y=0
		for i in range(ans.counter+1):
			if i==0:
				continue
			dc.setdefault(ans.heap[i].y,0)
			dc[ans.heap[i].y]+=1
			if dc[ans.heap[i].y]>maxx:
				maxx=dc[ans.heap[i].y]
				y=ans.heap[i].y
		return y
	def pred(self,test_x,k=None):
		'''
			test_x is a N*TM matrix,and indicate TM test case
			you can redecide the k
		'''
		if k==None:
			k=self.K
		test_case=np.array(test_x)
		y=[]
		for i in range(test_case.shape[1]):
			y.append(self.predict(test_case[:,i].transpose(),k))
		return y

       因為KNN畢竟是一個分類演算法，所以我在predict是加上了分類的程式碼，如果只想檢驗Kd-tree的話，你可以直接用for_point()找最近k個點

五.測試+後記

       測試：

       我們選取《統計學習方法》上面的例子：

       使用程式碼：

X=np.array([[2,5,9,4,8,7],[3,4,6,7,1,2]])
y=np.array([2,5,9,4,8,7])
knn=KNNC(X,1,y)
print knn.for_point([[6.5],[7]],1)

       輸出中後面帶了“===”的是掃描過的點，最後的是搜尋的結果：

       我們可以看到的確避免掃描了(2,3)，Bingo！！

       我們再knn.for_point([[2],[2]])：可以看到避免掃了很多點！！！

       後記：

       從實現寫此文前後耗時兩天，昨天寫程式碼寫到熄燈且剛好測試通過，怎一個爽字了得！！最後，再在github上求個Star

reference：

【2】《統計學習方法》 李航