wiki百科:http://zh.wikipedia.org/wiki/%E5%86%B3%E7%AD%96%E6%A0%91%E5%AD%A6%E4%B9%A0
html
opencv學習筆記--二杈決策樹:http://blog.csdn.net/homechao/article/details/9061921node
前言
前兩日,在微博上說:「到今天爲止,我至少虧欠了3篇文章待寫:一、KD樹;二、神經網絡;三、編程藝術第28章。你看到,blog內的文章與你於別處所見的任何都不一樣。因而,等啊等,等一臺電腦,只好等待..」。得益於田,借了我一臺電腦(借他電腦的時候,我連表示感謝,他說「能找到工做全靠你的博客,這點兒小忙還說,不地道」,有的時候,稍許感覺到受人信任也是一種壓力,願我不辜負你們對個人信任),因而今天開始Top 10 Algorithms in Data Mining系列第三篇文章,即本文「從K近鄰算法談到KD樹、SIFT+BBF算法」的創做。python
一我的堅持本身的興趣是比較難的,由於太多的人太容易爲外界所動了,而尤爲當你沒法從中獲得多少實際性的回報時,所幸,我能一直堅持下來。畢達哥拉斯學派有句名言:「萬物皆數」,最近讀完「微積分概念發展史」後也感覺到了這一點。同時,從算法到數據挖掘、機器學習,再到數學,其中每個領域任何一個細節都值得探索終生,或許,這就是「終生爲學」的意思。linux
本文各部份內容分佈以下:程序員
- 第一部分講K近鄰算法,其中重點闡述了相關的距離度量表示法,
- 第二部分着重講K近鄰算法的實現--KD樹,和KD樹的插入,刪除,最近鄰查找等操做,及KD樹的一系列相關改進(包括BBF,M樹等);
- 第三部分講KD樹的應用:SIFT+kd_BBF搜索算法。
同時,你將看到,K近鄰算法同本系列的前兩篇文章所講的決策樹分類貝葉斯分類,及支持向量機SVM同樣,也是用於解決分類問題的算法,web
而本數據挖掘十大算法系列也會按照分類,聚類,關聯分析,預測迴歸等問題依次展開闡述。面試
OK,行文倉促,本文如有任何漏洞,問題或者錯誤,歡迎朋友們隨時不吝指正,各位的批評也是我繼續寫下去的動力之一。感謝。算法
第一部分、K近鄰算法
1.一、什麼是K近鄰算法
何謂K近鄰算法,即K-Nearest Neighbor algorithm,簡稱KNN算法,單從名字來猜測,能夠簡單粗暴的認爲是:K個最近的鄰居,當K=1時,算法便成了最近鄰算法,即尋找最近的那個鄰居。爲什麼要找鄰居?打個比方來講,假設你來到一個陌生的村莊,如今你要找到與你有着類似特徵的人羣融入他們,所謂入夥。spring
用官方的話來講,所謂K近鄰算法,便是給定一個訓練數據集,對新的輸入實例,在訓練數據集中找到與該實例最鄰近的K個實例(也就是上面所說的K個鄰居),這K個實例的多數屬於某個類,就把該輸入實例分類到這個類中。根據這個說法,我們來看下引自維基百科上的一幅圖:shell
如上圖所示,有兩類不一樣的樣本數據,分別用藍色的小正方形和紅色的小三角形表示,而圖正中間的那個綠色的圓所標示的數據則是待分類的數據。也就是說,如今,咱們不知道中間那個綠色的數據是從屬於哪一類(藍色小正方形or紅色小三角形),下面,咱們就要解決這個問題:給這個綠色的圓分類。
咱們常說,物以類聚,人以羣分,判別一我的是一個什麼樣品質特徵的人,經常能夠從他/她身邊的朋友入手,所謂觀其友,而識其人。咱們不是要判別上圖中那個綠色的圓是屬於哪一類數據麼,好說,從它的鄰居下手。但一次性看多少個鄰居呢?從上圖中,你還能看到:
- 若是K=3,綠色圓點的最近的3個鄰居是2個紅色小三角形和1個藍色小正方形,少數從屬於多數,基於統計的方法,斷定綠色的這個待分類點屬於紅色的三角形一類。
- 若是K=5,綠色圓點的最近的5個鄰居是2個紅色三角形和3個藍色的正方形,仍是少數從屬於多數,基於統計的方法,斷定綠色的這個待分類點屬於藍色的正方形一類。
於此咱們看到,當沒法斷定當前待分類點是從屬於已知分類中的哪一類時,咱們能夠依據統計學的理論看它所處的位置特徵,衡量它周圍鄰居的權重,而把它歸爲(或分配)到權重更大的那一類。這就是K近鄰算法的核心思想。
1.二、近鄰的距離度量表示法
上文第一節,咱們看到,K近鄰算法的核心在於找到實例點的鄰居,這個時候,問題就接踵而至了,如何找到鄰居,鄰居的斷定標準是什麼,用什麼來度量。這一系列問題即是下面要講的距離度量表示法。但有的讀者可能就有疑問了,我是要找鄰居,找類似性,怎麼又跟距離扯上關係了?
這是由於特徵空間中兩個實例點的距離和反應出兩個實例點之間的類似性程度。K近鄰模型的特徵空間通常是n維實數向量空間,使用的距離能夠使歐式距離,也是能夠是其它距離,既然扯到了距離,下面就來具體闡述下都有哪些距離度量的表示法,權當擴展。
- 1. 歐氏距離,最多見的兩點之間或多點之間的距離表示法,又稱之爲歐幾里得度量,它定義於歐幾里得空間中,如點 x = (x1,...,xn) 和 y = (y1,...,yn) 之間的距離爲:
(1)二維平面上兩點a(x1,y1)與b(x2,y2)間的歐氏距離:
(2)三維空間兩點a(x1,y1,z1)與b(x2,y2,z2)間的歐氏距離:
(3)兩個n維向量a(x11,x12,…,x1n)與 b(x21,x22,…,x2n)間的歐氏距離:
也能夠用表示成向量運算的形式:
其上,二維平面上兩點歐式距離,代碼能夠以下編寫:
-
- double euclideanDistance(const vector<double>& v1, const vector<double>& v2)
- {
- assert(v1.size() == v2.size());
- double ret = 0.0;
- for (vector<double>::size_type i = 0; i != v1.size(); ++i)
- {
- ret += (v1[i] - v2[i]) * (v1[i] - v2[i]);
- }
- return sqrt(ret);
- }
- 2. 曼哈頓距離,咱們能夠定義曼哈頓距離的正式意義爲L1-距離或城市區塊距離,也就是在歐幾里得空間的固定直角座標系上兩點所造成的線段對軸產生的投影的距離總和。例如在平面上,座標(x1, y1)的點P1與座標(x2, y2)的點P2的曼哈頓距離爲:,要注意的是,曼哈頓距離依賴座標系統的轉度,而非系統在座標軸上的平移或映射。
通俗來說,想象你在曼哈頓要從一個十字路口開車到另一個十字路口,駕駛距離是兩點間的直線距離嗎?顯然不是,除非你能穿越大樓。而實際駕駛距離就是這個「曼哈頓距離」,此即曼哈頓距離名稱的來源, 同時,曼哈頓距離也稱爲城市街區距離(City Block distance)。
(1)二維平面兩點a(x1,y1)與b(x2,y2)間的曼哈頓距離
(2)兩個n維向量a(x11,x12,…,x1n)與 b(x21,x22,…,x2n)間的曼哈頓距離
- 3. 切比雪夫距離,若二個向量或二個點p 、and q,其座標分別爲及,則二者之間的切比雪夫距離定義以下:,
這也等於如下Lp度量的極值:
,所以切比雪夫距離也稱爲L∞度量。
以數學的觀點來看,切比雪夫距離是由一致範數(uniform norm)(或稱爲上確界範數)所衍生的度量,也是超凸度量(injective metric space)的一種。
在平面幾何中,若二點p及q的直角座標系座標爲
及
,則切比雪夫距離爲:
。
玩過國際象棋的朋友或許知道,國王走一步可以移動到相鄰的8個方格中的任意一個。那麼國王從格子(x1,y1)走到格子(x2,y2)最少須要多少步?。你會發現最少步數老是max( | x2-x1 | , | y2-y1 | ) 步 。有一種相似的一種距離度量方法叫切比雪夫距離。
(1)二維平面兩點a(x1,y1)與b(x2,y2)間的切比雪夫距離
(2)兩個n維向量a(x11,x12,…,x1n)與 b(x21,x22,…,x2n)間的切比雪夫距離
這個公式的另外一種等價形式是
- 4. 閔可夫斯基距離(Minkowski Distance),閔氏距離不是一種距離,而是一組距離的定義。
(1) 閔氏距離的定義
兩個n維變量a(x11,x12,…,x1n)與 b(x21,x22,…,x2n)間的閔可夫斯基距離定義爲:
-
5. 標準化歐氏距離 (Standardized Euclidean distance ),標準化歐氏距離是針對簡單歐氏距離的缺點而做的一種改進方案。標準歐氏距離的思路:既然數據各維份量的分佈不同,那先將各個份量都「標準化」到均值、方差相等。至於均值和方差標準化到多少,先複習點統計學知識。
假設樣本集X的數學指望或均值(mean)爲m,標準差(standard deviation,方差開根)爲s,那麼X的「標準化變量」X*表示爲:(X-m)/s,並且標準化變量的數學指望爲0,方差爲1。
即,樣本集的標準化過程(standardization)用公式描述就是:
標準化後的值 = ( 標準化前的值 - 份量的均值 ) /份量的標準差
通過簡單的推導就能夠獲得兩個n維向量a(x11,x12,…,x1n)與 b(x21,x22,…,x2n)間的標準化歐氏距離的公式:
若是將方差的倒數當作是一個權重,這個公式能夠當作是一種加權歐氏距離(Weighted Euclidean distance)。
-
6. 馬氏距離(Mahalanobis Distance)
(1)馬氏距離定義
有M個樣本向量X1~Xm,
協方差矩陣記爲S,均值記爲向量μ,則其中樣本向量X到u的馬氏距離表示爲:
(
協方差矩陣中每一個元素是各個矢量元素之間的協方差Cov(X,Y),Cov(X,Y) = E{ [X-E(X)] [Y-E(Y)]},其中E爲數學指望)
而其中向量Xi與Xj之間的馬氏距離定義爲:
若協方差矩陣是單位矩陣(各個樣本向量之間獨立同分布),則公式就成了:
也就是歐氏距離了。
若協方差矩陣是對角矩陣,公式變成了標準化歐氏距離。
(2)馬氏距離的優缺點:量綱無關,排除變量之間的相關性的干擾。
「
微博上的seafood高清版點評道:原來馬氏距離是根據協方差矩陣演變,一直被老師誤導了,怪不得看Killian在05年NIPS發表的LMNN論文時候總是看到協方差矩陣和半正定,原來是這回事」
- 七、巴氏距離(Bhattacharyya Distance),在統計中,Bhattacharyya距離測量兩個離散或連續機率分佈的類似性。它與衡量兩個統計樣品或種羣之間的重疊量的Bhattacharyya係數密切相關。Bhattacharyya距離和Bhattacharyya係數以20世紀30年代曾在印度統計研究所工做的一個統計學家A. Bhattacharya命名。同時,Bhattacharyya係數能夠被用來肯定兩個樣本被認爲相對接近的,它是用來測量中的類分類的可分離性。
對於離散機率分佈 p和q在同一域 X,它被定義爲:
對於連續機率分佈,Bhattacharyya係數被定義爲:
在
這兩種狀況下,巴氏距離
並無服從三角不等式.(值得一提的是,Hellinger距離不服從三角不等式
)。
對於多變量的高斯分佈
,
,
和是手段和協方差的分佈
。
須要注意的是,在這種狀況下,第一項中的Bhattacharyya距離與馬氏距離有關聯。
Bhattacharyya係數是兩個統計樣本之間的重疊量的近似測量,能夠被用於肯定被考慮的兩個樣本的相對接近。
計算Bhattacharyya係數涉及集成的基本形式的兩個樣本的重疊的時間間隔的值的兩個樣本被分裂成一個選定的分區數,而且在每一個分區中的每一個樣品的成員的數量,在下面的公式中使用
- 8. 漢明距離(Hamming distance), 兩個等長字符串s1與s2之間的漢明距離定義爲將其中一個變爲另一個所須要做的最小替換次數。例如字符串「1111」與「1001」之間的漢明距離爲2。應用:信息編碼(爲了加強容錯性,應使得編碼間的最小漢明距離儘量大)。
或許,你還沒明白我再說什麼,不急,看下
上篇blog中第78題的第3小題整理的一道面試題目,便一目瞭然了。以下圖所示:
-
-
-
- f[i,j] = min { f[i-1,j]+1, f[i,j-1]+1, f[i-1,j-1]+(s[i]==t[j]?0:1) }
-
-
(上篇blog中第78題的第3小題給出了多種方法,讀者能夠參看之。同時,程序員編程藝術系列第二十八章將詳細闡述這個問題)
- 9. 夾角餘弦(Cosine) ,幾何中夾角餘弦可用來衡量兩個向量方向的差別,機器學習中借用這一律念來衡量樣本向量之間的差別。
(1)在二維空間中向量A(x1,y1)與向量B(x2,y2)的夾角餘弦公式:
(2) 兩個n維樣本點a(x11,x12,…,x1n)和b(x21,x22,…,x2n)的夾角餘弦
相似的,對於兩個n維樣本點a(x11,x12,…,x1n)和b(x21,x22,…,x2n),能夠使用相似於夾角餘弦的概念來衡量它們間的類似程度,即:
夾角餘弦取值範圍爲[-1,1]。夾角餘弦越大表示兩個向量的夾角越小,夾角餘弦越小表示兩向量的夾角越大。當兩個向量的方向重合時夾角餘弦取最大值1,當兩個向量的方向徹底相反夾角餘弦取最小值-1。
- 10. 傑卡德類似係數(Jaccard similarity coefficient)
兩個集合A和B的交集元素在A,B的並集中所佔的比例,稱爲兩個集合的傑卡德類似係數,用符號J(A,B)表示。
與傑卡德類似係數相反的概念是傑卡德距離(Jaccard distance)。
傑卡德距離用兩個集合中不一樣元素佔全部元素的比例來衡量兩個集合的區分度。
可將傑卡德類似係數用在衡量樣本的類似度上。
舉例:樣本A與樣本B是兩個n維向量,並且全部維度的取值都是0或1,例如:A(0111)和B(1011)。咱們將樣本當作是一個集合,1表示集合包含該元素,0表示集合不包含該元素。
依據上文給的傑卡德類似係數及傑卡德距離的相關定義,樣本A與B的傑卡德類似係數J能夠表示爲:
這裏M
11+
M
01+
M
10可理解爲A與B的並集的元素個數,而M
11是A與B的交集的元素個數。而樣本A與B的傑卡德距離表示爲J':
- 11.皮爾遜係數(Pearson Correlation Coefficient)
在具體闡述皮爾遜相關係數以前,有必要解釋下什麼是相關係數 ( Correlation coefficient )與相關距離(Correlation distance)。
相關係數 ( Correlation coefficient )的定義是:
(
其中,E爲數學指望或均值,D爲方差,D開根號爲標準差,E{ [X-E(X)] [Y-E(Y)]}稱爲隨機變量X與Y的協方差,記爲Cov(X,Y),即Cov(X,Y) = E{ [X-E(X)] [Y-E(Y)]},而兩個變量之間的協方差和標準差的商則稱爲隨機變量X與Y的相關係數,記爲)
相關係數衡量隨機變量X與Y相關程度的一種方法,相關係數的取值範圍是[-1,1]。相關係數的絕對值越大,則代表X與Y相關度越高。當X與Y線性相關時,相關係數取值爲1(正線性相關)或-1(負線性相關)。
具體的,若是有兩個變量:X、Y,最終計算出的相關係數的含義能夠有以下理解:
- 當相關係數爲0時,X和Y兩變量無關係。
- 當X的值增大(減少),Y值增大(減少),兩個變量爲正相關,相關係數在0.00與1.00之間。
- 當X的值增大(減少),Y值減少(增大),兩個變量爲負相關,相關係數在-1.00與0.00之間。
相關距離的定義是:
OK,接下來,我們來重點了解下皮爾遜相關係數。
在統計學中,皮爾遜積矩相關係數(英語:Pearson product-moment correlation coefficient,又稱做 PPMCC或PCCs, 用r表示)用於度量兩個變量X和Y之間的相關(線性相關),其值介於-1與1之間。
一般狀況下經過如下取值範圍判斷變量的相關強度:
相關係數 0.8-1.0 極強相關
0.6-0.8 強相關
0.4-0.6 中等程度相關
0.2-0.4 弱相關
0.0-0.2 極弱相關或無相關
在天然科學領域中,該係數普遍用於度量兩個變量之間的相關程度。它是由卡爾·皮爾遜從弗朗西斯·高爾頓在19世紀80年代提出的一個類似卻又稍有不一樣的想法演變而來的。這個相關係數也稱做「皮爾森相關係數r」。
(1)皮爾遜係數的定義:
兩個變量之間的皮爾遜相關係數定義爲兩個變量之間的協方差和標準差的商:
以上方程定義了整體相關係數, 通常表示成希臘字母ρ(rho)。基於樣本對協方差和方差進行估計,能夠獲得樣本標準差, 通常表示成r:
一種等價表達式的是表示成標準分的均值。基於(Xi, Yi)的樣本點,樣本皮爾遜係數是
其中
、
及
,分別是標準分、樣本平均值和樣本標準差。
或許上面的講解令你頭腦混亂不堪,不要緊,我換一種方式講解,以下:
假設有兩個變量X、Y,那麼兩變量間的皮爾遜相關係數可經過如下公式計算:
- 公式一:
注:勿忘了上面說過,「皮爾遜相關係數定義爲兩個變量之間的協方差和標準差的商」,其中標準差的計算公式爲:
以上列出的四個公式等價,其中E是數學指望,cov表示協方差,N表示變量取值的個數。
(2)皮爾遜相關係數的適用範圍
當兩個變量的標準差都不爲零時,相關係數纔有定義,皮爾遜相關係數適用於:
- 兩個變量之間是線性關係,都是連續數據。
- 兩個變量的整體是正態分佈,或接近正態的單峯分佈。
- 兩個變量的觀測值是成對的,每對觀測值之間相互獨立。
(3)如何理解皮爾遜相關係數
rubyist:皮爾遜相關係數理解有兩個角度
其一, 按照高中數學水平來理解, 它很簡單, 能夠看作將兩組數據首先作Z分數處理以後, 而後兩組數據的乘積和除以樣本數,Z分數通常表明正態分佈中, 數據偏離中心點的距離.等於變量減掉平均數再除以標準差.(就是高考的標準分相似的處理)
樣本標準差則等於變量減掉平均數的平方和,再除以樣本數,最後再開方,也就是說,方差開方即爲標準差,樣本標準差計算公式爲:
因此, 根據這個最樸素的理解,咱們能夠將公式依次精簡爲:
其二, 按照大學的線性數學水平來理解, 它比較複雜一點,能夠看作是兩組數據的向量夾角的餘弦。下面是關於此皮爾遜係數的幾何學的解釋,先來看一幅圖,以下所示:
迴歸直線: y=gx(x) [紅色] 和 x=gy(y) [藍色]
如上圖,對於沒有中心化的數據, 相關係數與兩條可能的迴歸線y=gx(x) 和 x=gy(y) 夾角的餘弦值一致。
對於沒有中心化的數據 (也就是說, 數據移動一個樣本平均值以使其均值爲0), 相關係數也能夠被視做由兩個隨機變量 向量 夾角 的 餘弦值(見下方)。
舉個例子,例如,有5個國家的國民生產總值分別爲 10, 20, 30, 50 和 80 億美圓。 假設這5個國家 (順序相同) 的貧困百分比分別爲 11%, 12%, 13%, 15%, and 18% 。 令 x 和 y 分別爲包含上述5個數據的向量: x = (1, 2, 3, 5, 8) 和 y = (0.11, 0.12, 0.13, 0.15, 0.18)。
利用一般的方法計算兩個向量之間的夾角 (參見 數量積), 未中心化 的相關係數是:
咱們發現以上的數據特地選定爲徹底相關: y = 0.10 + 0.01 x。 因而,皮爾遜相關係數應該等於1。將數據中心化 (經過E(x) = 3.8移動 x 和經過 E(y) = 0.138 移動 y ) 獲得 x = (−2.8, −1.8, −0.8, 1.2, 4.2) 和 y = (−0.028, −0.018, −0.008, 0.012, 0.042), 從中
(4)皮爾遜相關的約束條件
從以上解釋, 也能夠理解皮爾遜相關的約束條件:
- 1 兩個變量間有線性關係
- 2 變量是連續變量
- 3 變量均符合正態分佈,且二元分佈也符合正態分佈
- 4 兩變量獨立
在實踐統計中,通常只輸出兩個係數,一個是相關係數,也就是計算出來的相關係數大小,在-1到1之間;另外一個是獨立樣本檢驗係數,用來檢驗樣本一致性。
簡單說來,各類「距離」的應用場景簡單歸納爲,空間:歐氏距離,路徑:曼哈頓距離,國際象棋國王:切比雪夫距離,以上三種的統一形式:閔可夫斯基距離,加權:標準化歐氏距離,排除量綱和依存:馬氏距離,向量差距:夾角餘弦,編碼差異:漢明距離,集合近似度:傑卡德相似係數與距離,相關:相關係數與相關距離。
1.三、K值的選擇
除了上述1.2節如何定義鄰居的問題以外,還有一個選擇多少個鄰居,即K值定義爲多大的問題。不要小看了這個K值選擇問題,由於它對K近鄰算法的結果會產生重大影響。如李航博士的一書「統計學習方法」上所說:
- 若是選擇較小的K值,就至關於用較小的領域中的訓練實例進行預測,「學習」近似偏差會減少,只有與輸入實例較近或類似的訓練實例纔會對預測結果起做用,與此同時帶來的問題是「學習」的估計偏差會增大,換句話說,K值的減少就意味着總體模型變得複雜,容易發生過擬合;
- 若是選擇較大的K值,就至關於用較大領域中的訓練實例進行預測,其優勢是能夠減小學習的估計偏差,但缺點是學習的近似偏差會增大。這時候,與輸入實例較遠(不類似的)訓練實例也會對預測器做用,使預測發生錯誤,且K值的增大就意味着總體的模型變得簡單。
- K=N,則徹底不足取,由於此時不管輸入實例是什麼,都只是簡單的預測它屬於在訓練實例中最多的累,模型過於簡單,忽略了訓練實例中大量有用信息。
在實際應用中,K值通常取一個比較小的數值,例如採用
交叉驗證法(簡單來講,就是一部分樣本作訓練集,一部分作測試集)來選擇最優的K值。
第二部分、K近鄰算法的實現:KD樹
2.0、背景
以前blog內曾經介紹過SIFT特徵匹配算法,特徵點匹配和數據庫查、圖像檢索本質上是同一個問題,均可以歸結爲一個經過距離函數在高維矢量之間進行類似性檢索的問題,如何快速而準確地找到查詢點的近鄰,很多人提出了不少高維空間索引結構和近似查詢的算法。
通常說來,索引結構中類似性查詢有兩種基本的方式:
- 一種是範圍查詢,範圍查詢時給定查詢點和查詢距離閾值,從數據集中查找全部與查詢點距離小於閾值的數據
- 另外一種是K近鄰查詢,就是給定查詢點及正整數K,從數據集中找到距離查詢點最近的K個數據,當K=1時,它就是最近鄰查詢。
一樣,針對特徵點匹配也有兩種方法:
- 最容易的辦法就是線性掃描,也就是咱們常說的窮舉搜索,依次計算樣本集E中每一個樣本到輸入實例點的距離,而後抽取出計算出來的最小距離的點即爲最近鄰點。此種辦法簡單直白,但當樣本集或訓練集很大時,它的缺點就立馬暴露出來了,舉個例子,在物體識別的問題中,可能有數千個甚至數萬個SIFT特徵點,而去一一計算這成千上萬的特徵點與輸入實例點的距離,明顯是不足取的。
- 另一種,就是構建數據索引,由於實際數據通常都會呈現簇狀的聚類形態,所以咱們想到創建數據索引,而後再進行快速匹配。索引樹是一種樹結構索引方法,其基本思想是對搜索空間進行層次劃分。根據劃分的空間是否有混疊能夠分爲Clipping和Overlapping兩種。前者劃分空間沒有重疊,其表明就是k-d樹;後者劃分空間相互有交疊,其表明爲R樹。
而關於R樹本blog內以前已有介紹(同時,關於基於R樹的最近鄰查找,還能夠看下這篇文章:http://blog.sina.com.cn/s/blog_72e1c7550101dsc3.html),本文着重介紹k-d樹。
1975年,來自斯坦福大學的Jon Louis Bentley在ACM雜誌上發表的一篇論文:Multidimensional Binary Search Trees Used for Associative Searching 中正式提出和闡述的了以下圖形式的把空間劃分爲多個部分的k-d樹。
2.一、什麼是KD樹
Kd-樹是K-dimension tree的縮寫,是對數據點在k維空間(如二維(x,y),三維(x,y,z),k維(x1,y,z..))中劃分的一種數據結構,主要應用於多維空間關鍵數據的搜索(如:範圍搜索和最近鄰搜索)。本質上說,Kd-樹就是一種平衡二叉樹。
首先必須搞清楚的是,k-d樹是一種空間劃分樹,說白了,就是把整個空間劃分爲特定的幾個部分,而後在特定空間的部份內進行相關搜索操做。想像一個三維(多維有點爲難你的想象力了)空間,kd樹按照必定的劃分規則把這個三維空間劃分了多個空間,以下圖所示:
2.二、KD樹的構建
kd樹構建的僞代碼以下圖所示:
再舉一個簡單直觀的實例來介紹k-d樹構建算法。假設有6個二維數據點{(2,3),(5,4),(9,6),(4,7),(8,1),(7,2)},數據點位於二維空間內,以下圖所示。爲了能有效的找到最近鄰,k-d樹採用分而治之的思想,即將整個空間劃分爲幾個小部分,首先,粗黑線將空間一分爲二,而後在兩個子空間中,細黑直線又將整個空間劃分爲四部分,最後虛黑直線將這四部分進一步劃分。
6個二維數據點{(2,3),(5,4),(9,6),(4,7),(8,1),(7,2)}構建kd樹的具體步驟爲:
- 肯定:split域=x。具體是:6個數據點在x,y維度上的數據方差分別爲39,28.63,因此在x軸上方差更大,故split域值爲x;
- 肯定:Node-data = (7,2)。具體是:根據x維上的值將數據排序,6個數據的中值(所謂中值,即中間大小的值)爲7,因此Node-data域位數據點(7,2)。這樣,該節點的分割超平面就是經過(7,2)並垂直於:split=x軸的直線x=7;
- 肯定:左子空間和右子空間。具體是:分割超平面x=7將整個空間分爲兩部分:x<=7的部分爲左子空間,包含3個節點={(2,3),(5,4),(4,7)};另外一部分爲右子空間,包含2個節點={(9,6),(8,1)};
如上算法所述,kd樹的構建是一個遞歸過程,咱們對左子空間和右子空間內的數據重複根節點的過程就能夠獲得一級子節點(5,4)和(9,6),同時將空間和數據集進一步細分,如此往復直到空間中只包含一個數據點。
與此同時,通過對上面所示的空間劃分以後,咱們能夠看出,點(7,2)能夠爲根結點,從根結點出發的兩條紅粗斜線指向的(5,4)和(9,6)則爲根結點的左右子結點,而(2,3),(4,7)則爲(5,4)的左右孩子(經過兩條細紅斜線相連),最後,(8,1)爲(9,6)的左孩子(經過細紅斜線相連)。如此,便造成了下面這樣一棵k-d樹:
k-d樹的數據結構
針對上表給出的kd樹的數據結構,轉化成具體代碼以下所示(注,本文如下代碼分析基於Rob Hess維護的sift庫):
-
- struct kd_node
- {
- int ki; //關鍵點直方圖方差最大向量系列位置
- double kv; //直方圖方差最大向量系列中最中間模值
- int leaf;
- struct feature* features;
- int n;
- struct kd_node* kd_left;
- struct kd_node* kd_right;
- };
也就是說,如以前所述,kd樹中,kd表明k-dimension,每一個節點即爲一個k維的點。每一個非葉節點能夠想象爲一個分割超平面,用垂直於座標軸的超平面將空間分爲兩個部分,這樣遞歸的從根節點不停的劃分,直到沒有實例爲止。經典的構造k-d tree的規則以下:
- 隨着樹的深度增長,循環的選取座標軸,做爲分割超平面的法向量。對於3-d tree來講,根節點選取x軸,根節點的孩子選取y軸,根節點的孫子選取z軸,根節點的曾孫子選取x軸,這樣循環下去。
- 每次均爲全部對應實例的中位數的實例做爲切分點,切分點做爲父節點,左右兩側爲劃分的做爲左右兩子樹。
對於n個實例的k維數據來講,創建kd-tree的時間複雜度爲O(k*n*logn)。
如下是構建k-d樹的代碼:
- struct kd_node* kdtree_build( struct feature* features, int n )
- {
- struct kd_node* kd_root;
-
- if( ! features || n <= 0 )
- {
- fprintf( stderr, "Warning: kdtree_build(): no features, %s, line %d\n",
- __FILE__, __LINE__ );
- return NULL;
- }
-
-
- kd_root = kd_node_init( features, n );
- expand_kd_node_subtree( kd_root );
-
- return kd_root;
- }
上面的涉及初始化操做的兩個函數kd_node_init,及expand_kd_node_subtree代碼分別以下所示:
- static struct kd_node* kd_node_init( struct feature* features, int n )
- {
- struct kd_node* kd_node;
-
- kd_node = (struct kd_node*)(malloc( sizeof( struct kd_node ) ));
- memset( kd_node, 0, sizeof( struct kd_node ) );
- kd_node->ki = -1;
- kd_node->features = features;
- kd_node->n = n;
-
- return kd_node;
- }
- static void expand_kd_node_subtree( struct kd_node* kd_node )
- {
-
- if( kd_node->n == 1 || kd_node->n == 0 )
- {
- kd_node->leaf = 1;
- return;
- }
-
- assign_part_key( kd_node );
- partition_features( kd_node );
-
- if( kd_node->kd_left )
- expand_kd_node_subtree( kd_node->kd_left );
- if( kd_node->kd_right )
- expand_kd_node_subtree( kd_node->kd_right );
- }
構建完kd樹以後,現在進行最近鄰搜索呢?從下面的動態gif圖中,你是否能看出些許端倪呢?
k-d樹算法能夠分爲兩大部分,除了上部分有關k-d樹自己這種數據結構創建的算法,另外一部分是在創建的k-d樹上各類諸如插入,刪除,查找(最鄰近查找)等操做涉及的算法。下面,我們依次來看kd樹的插入、刪除、查找操做。
2.三、KD樹的插入
元素插入到一個K-D樹的方法和二叉檢索樹相似。本質上,在偶數層比較x座標值,而在奇數層比較y座標值。當咱們到達了樹的底部,(也就是當一個空指針出現),咱們也就找到告終點將要插入的位置。生成的K-D樹的形狀依賴於結點插入時的順序。給定N個點,其中一個結點插入和檢索的平均代價是O(log2N)。
下面4副圖(來源:中國地質大學電子課件)說明了插入順序爲(a) Chicago, (b) Mobile, (c) Toronto, and (d) Buffalo,創建空間K-D樹的示例:
應該清楚,這裏描述的插入過程當中,每一個結點將其所在的平面分割成兩部分。因比,Chicago 將平面上全部結點分紅兩部分,一部分全部的結點x座標值小於35,另外一部分結點的x座標值大於或等於35。一樣Mobile將全部x座標值大於35的結點以分紅兩部分,一部分結點的Y座標值是小於10,另外一部分結點的Y座標值大於或等於10。後面的Toronto、Buffalo也按照一分爲二的規則繼續劃分。
2.四、KD樹的刪除
KD樹的刪除能夠用遞歸程序來實現。咱們假設但願從K-D樹中刪除結點(a,b)。若是(a,b)的兩個子樹都爲空,則用空樹來代替(a,b)。不然,在(a,b)的子樹中尋找一個合適的結點來代替它,譬如(c,d),則遞歸地從K-D樹中刪除(c,d)。一旦(c,d)已經被刪除,則用(c,d)代替(a,b)。假設(a,b)是一個X識別器,那麼,
它得替代節點要麼是(a,b)左子樹中的X座標最大值的結點,要麼是(a,b)右子樹中x座標最小值的結點。
也就是說,跟普通二叉樹(
包括以下圖所示的紅黑樹)結點的刪除是一樣的思想:用被刪除節點A的左子樹的最右節點或者A的右子樹的最左節點做爲替代A的節點(
好比,下圖紅黑樹中,若要刪除根結點26,第一步即是用23或28取代根結點26)。
當(a,b)的右子樹爲空時,找到(a,b)左子樹中具備x座標最大的結點,譬如(c,d),將(a,b)的左子樹放到(c,d)的右子樹中,且在樹中從它的上一層遞歸地應用刪除過程(也就是(a,b)的左子樹) 。
下面來舉一個實際的例子(
來源:中國地質大學電子課件,原課件錯誤已經在下文中訂正),以下圖所示,原始圖像及對應的kd樹,如今要刪除圖中的A結點,請看一系列刪除步驟:
要刪除上圖中結點A,選擇結點A的右子樹中X座標值最小的結點,這裏是C,C成爲根,以下圖:
從C的右子樹中找出一個結點代替先前C的位置,
這裏是D,並將D的左子樹轉爲它的右子樹,D代替先前C的位置,以下圖:
在D的新右子樹中,找X座標最小的結點,這裏爲H,H代替D的位置,
在D的右子樹中找到一個Y座標最小的值,這裏是I,將I代替原先H的位置,從而A結點從圖中順利刪除,以下圖所示:
從一個K-D樹中刪除結點(a,b)的問題變成了在(a,b)的子樹中尋找x座標爲最小的結點。不幸的是尋找最小x座標值的結點比二叉檢索樹中解決相似的問題要複雜得多。特別是雖然最小x座標值的結點必定在x識別器的左子樹中,但它一樣可在y識別器的兩個子樹中。所以關係到檢索,且必須注意檢索座標,以使在每一個奇數層僅檢索2個子樹中的一個。
從K-D樹中刪除一個結點是代價很高的,很清楚刪除子樹的根受到子樹中結點個數的限制。用TPL(T)表示樹T總的路徑長度。可看出樹中子樹大小的總和爲TPL(T)+N。 以隨機方式插入N個點造成樹的TPL是O(N*log2N),這就意味着從一個隨機造成的K-D樹中刪除一個隨機選取的結點平均代價的上界是O(log2N) 。
2.五、KD樹的最近鄰搜索算法
現實生活中有許多問題須要在多維數據的快速分析和快速搜索,對於這個問題最經常使用的方法是所謂的kd樹。在k-d樹中進行數據的查找也是特徵匹配的重要環節,其目的是檢索在k-d樹中與查詢點距離最近的數據點。在一個N維的笛卡兒空間在兩個點之間的距離是由下述公式肯定:
2.5.一、k-d樹查詢算法的僞代碼
k-d樹查詢算法的僞代碼以下所示:
- 算法:k-d樹最鄰近查找
- 輸入:Kd,
- target
- 輸出:nearest,
- dist
-
- 1. If Kd爲NULL,則設dist爲infinite並返回
- 2.
- Kd_point = &Kd;
- nearest = Kd_point -> Node-data;
-
- while(Kd_point)
- push(Kd_point)到search_path中;
-
- If Dist(nearest,target) > Dist(Kd_point -> Node-data,target)
- nearest = Kd_point -> Node-data;
- Min_dist = Dist(Kd_point,target);
- s = Kd_point -> split;
-
- If target[s] <= Kd_point -> Node-data[s]
- Kd_point = Kd_point -> left;
- else
- Kd_point = Kd_point ->right;
- End while
-
- 3.
- while(search_path != NULL)
- back_point = 從search_path取出一個節點指針;
- s = back_point -> split;
-
- If Dist(target[s],back_point -> Node-data[s]) < Max_dist
- If target[s] <= back_point -> Node-data[s]
- Kd_point = back_point -> right;
- else
- Kd_point = back_point -> left;
- 將Kd_point壓入search_path堆棧;
-
- If Dist(nearest,target) > Dist(Kd_Point -> Node-data,target)
- nearest = Kd_point -> Node-data;
- Min_dist = Dist(Kd_point -> Node-data,target);
- End while
讀者來信點評@yhxyhxyhx,在「將Kd_point壓入search_path堆棧;」這行代碼後,應該是調到步驟2再往下走二分搜索的邏輯一直到葉結點,我寫了一個遞歸版本的二維kd tree的搜索函數你對比的看看:
- void innerGetClosest(NODE* pNode, PT point, PT& res, int& nMinDis)
- {
- if (NULL == pNode)
- return;
- int nCurDis = abs(point.x - pNode->pt.x) + abs(point.y - pNode->pt.y);
- if (nMinDis < 0 || nCurDis < nMinDis)
- {
- nMinDis = nCurDis;
- res = pNode->pt;
- }
- if (pNode->splitX && point.x <= pNode->pt.x || !pNode->splitX && point.y <= pNode->pt.y)
- innerGetClosest(pNode->pLft, point, res, nMinDis);
- else
- innerGetClosest(pNode->pRgt, point, res, nMinDis);
- int rang = pNode->splitX ? abs(point.x - pNode->pt.x) : abs(point.y - pNode->pt.y);
- if (rang > nMinDis)
- return;
- NODE* pGoInto = pNode->pLft;
- if (pNode->splitX && point.x > pNode->pt.x || !pNode->splitX && point.y > pNode->pt.y)
- pGoInto = pNode->pRgt;
- innerGetClosest(pGoInto, point, res, nMinDis);
- }
下面,以兩個簡單的實例(例子來自圖像局部不變特性特徵與描述一書)來描述最鄰近查找的基本思路。
2.5.二、舉例:查詢點(2.1,3.1)
星號表示要查詢的點(2.1,3.1)。經過二叉搜索,順着搜索路徑很快就能找到最鄰近的近似點,也就是葉子節點(2,3)。而找到的葉子節點並不必定就是最鄰近的,最鄰近確定距離查詢點更近,應該位於以查詢點爲圓心且經過葉子節點的圓域內。爲了找到真正的最近鄰,還須要進行相關的‘回溯'操做。也就是說,算法首先沿搜索路徑反向查找是否有距離查詢點更近的數據點。
以查詢(2.1,3.1)爲例:
- 二叉樹搜索:先從(7,2)點開始進行二叉查找,而後到達(5,4),最後到達(2,3),此時搜索路徑中的節點爲<(7,2),(5,4),(2,3)>,首先以(2,3)做爲當前最近鄰點,計算其到查詢點(2.1,3.1)的距離爲0.1414,
- 回溯查找:在獲得(2,3)爲查詢點的最近點以後,回溯到其父節點(5,4),並判斷在該父節點的其餘子節點空間中是否有距離查詢點更近的數據點。以(2.1,3.1)爲圓心,以0.1414爲半徑畫圓,以下圖所示。發現該圓並不和超平面y = 4交割,所以不用進入(5,4)節點右子空間中(圖中灰色區域)去搜索;
- 最後,再回溯到(7,2),以(2.1,3.1)爲圓心,以0.1414爲半徑的圓更不會與x = 7超平面交割,所以不用進入(7,2)右子空間進行查找。至此,搜索路徑中的節點已經所有回溯完,結束整個搜索,返回最近鄰點(2,3),最近距離爲0.1414。
2.5.三、舉例:查詢點(2,4.5)
一個複雜點了例子如查找點爲(2,4.5),具體步驟依次以下:
- 一樣先進行二叉查找,先從(7,2)查找到(5,4)節點,在進行查找時是由y = 4爲分割超平面的,因爲查找點爲y值爲4.5,所以進入右子空間查找到(4,7),造成搜索路徑<(7,2),(5,4),(4,7)>,但(4,7)與目標查找點的距離爲3.202,而(5,4)與查找點之間的距離爲3.041,因此(5,4)爲查詢點的最近點;
- 以(2,4.5)爲圓心,以3.041爲半徑做圓,以下圖所示。可見該圓和y = 4超平面交割,因此須要進入(5,4)左子空間進行查找,也就是將(2,3)節點加入搜索路徑中得<(7,2),(2,3)>;因而接着搜索至(2,3)葉子節點,(2,3)距離(2,4.5)比(5,4)要近,因此最近鄰點更新爲(2,3),最近距離更新爲1.5;
- 回溯查找至(5,4),直到最後回溯到根結點(7,2)的時候,以(2,4.5)爲圓心1.5爲半徑做圓,並不和x = 7分割超平面交割,以下圖所示。至此,搜索路徑回溯完,返回最近鄰點(2,3),最近距離1.5。
上述兩次實例代表,當查詢點的鄰域與分割超平面兩側空間交割時,須要查找另外一側子空間,致使檢索過程複雜,效率降低。
通常來說,最臨近搜索只須要檢測幾個葉子結點便可,以下圖所示:
可是,若是當實例點的分佈比較糟糕時,幾乎要遍歷全部的結點,以下所示:
研究代表N個節點的K維k-d樹搜索過程時間複雜度爲:tworst=O(kN1-1/k)。
同時,以上爲了介紹方便,討論的是二維或三維情形。但在實際的應用中,如SIFT特徵矢量128維,SURF特徵矢量64維,維度都比較大,直接利用k-d樹快速檢索(維數不超過20)的性能急劇降低,幾乎接近貪婪線性掃描。假設數據集的維數爲D,通常來講要求數據的規模N知足N»2D,才能達到高效的搜索。因此這就引出了一系列對k-d樹算法的改進:BBF算法,和一系列M樹、VP樹、MVP樹等高維空間索引樹(下文2.6節kd樹近鄰搜索算法的改進:BBF算法,與2.7節球樹、M樹、VP樹、MVP樹)。
2.六、kd樹近鄰搜索算法的改進:BBF算法
我們順着上一節的思路,參考統計學習方法一書上的內容,再來總結下kd樹的最近鄰搜索算法:
輸入:以構造的kd樹,目標點x;
輸出:x 的最近鄰
算法步驟以下:
- 在kd樹種找出包含目標點x的葉結點:從根結點出發,遞歸地向下搜索kd樹。若目標點x當前維的座標小於切分點的座標,則移動到左子結點,不然移動到右子結點,直到子結點爲葉結點爲止。
- 以此葉結點爲「當前最近點」。
- 遞歸的向上回溯,在每一個結點進行如下操做:
(a)若是該結點保存的實例點比當前最近點距離目標點更近,則更新「當前最近點」,也就是說以該實例點爲「當前最近點」。
(b)當前最近點必定存在於該結點一個子結點對應的區域,檢查子結點的父結點的另外一子結點對應的區域是否有更近的點。具體作法是,檢查另外一子結點對應的區域是否以目標點位球心,以目標點與「當前最近點」間的距離爲半徑的圓或超球體相交:
若是相交,可能在另外一個子結點對應的區域內存在距目標點更近的點,移動到另外一個子結點,接着,繼續遞歸地進行最近鄰搜索;
若是不相交,向上回溯。
- 當回退到根結點時,搜索結束,最後的「當前最近點」即爲x 的最近鄰點。
若是實例點是隨機分佈的,那麼kd樹搜索的平均計算複雜度是O(NlogN),這裏的N是訓練實例樹。因此說,kd樹更適用於訓練實例數遠大於空間維數時的k近鄰搜索,當空間維數接近訓練實例數時,它的效率會迅速降低,一降降到「解放前」:線性掃描的速度。
也正由於上述k最近鄰搜索算法的第4個步驟中的所述:「回退到根結點時,搜索結束」,每一個最近鄰點的查詢比較完成過程最終都要回退到根結點而結束,而致使了許多沒必要要回溯訪問和比較到的結點,這些多餘的損耗在高維度數據查找的時候,搜索效率將變得至關之地下,那有什麼辦法能夠改進這個原始的kd樹最近鄰搜索算法呢?
從上述標準的kd樹查詢過程能夠看出其搜索過程當中的「回溯」是由「查詢路徑」決定的,並無考慮查詢路徑上一些數據點自己的一些性質。一個簡單的改進思路就是將「查詢路徑」上的結點進行排序,如按各自分割超平面(也稱bin)與查詢點的距離排序,也就是說,回溯檢查老是從優先級最高(Best Bin)的樹結點開始。
針對此BBF機制,讀者Feng&書童點評道:
- 在某一層,分割面是第ki維,分割值是kv,那麼 abs(q[ki]-kv) 就是沒有選擇的那個分支的優先級,也就是計算的是那一維上的距離;
- 同時,從優先隊列裏面取節點只在某次搜索到葉節點後才發生,計算過距離的節點不會出如今隊列的,好比1~10這10個節點,你第一次搜索到葉節點的路徑是1-5-7,那麼1,5,7是不會出如今優先隊列的。換句話說,優先隊列裏面存的都是查詢路徑上節點對應的相反子節點,好比:搜索左子樹,就把對應這一層的右節點存進隊列。
如此,就引出了本節要討論的kd樹最近鄰搜索算法的改進:BBF(Best-Bin-First)查詢算法,它是由發明sift算法的David Lowe在1997的一篇文章中針對高維數據提出的一種近似算法,此算法能確保優先檢索包含最近鄰點可能性較高的空間,此外,BBF機制還設置了一個運行超時限定。採用了BBF查詢機制後,kd樹即可以有效的擴展到高維數據集上。
僞代碼以下圖所示(圖取自圖像局部不變特性特徵與描述一書):
仍是以上面的查詢(2,4.5)爲例,搜索的算法流程爲:
- 將(7,2)壓人優先隊列中;
- 提取優先隊列中的(7,2),因爲(2,4.5)位於(7,2)分割超平面的左側,因此檢索其左子結點(5,4)。同時,根據BBF機制」搜索左/右子樹,就把對應這一層的兄弟結點即右/左結點存進隊列」,將其(5,4)對應的兄弟結點即右子結點(9,6)壓人優先隊列中,此時優先隊列爲{(9,6)},最佳點爲(7,2);而後一直檢索到葉子結點(4,7),此時優先隊列爲{(2,3),(9,6)},「最佳點」則爲(5,4);
- 提取優先級最高的結點(2,3),重複步驟2,直到優先隊列爲空。
如你在下圖所見到的那樣(話說,用鼠標在圖片上寫字着實很差寫):
2.七、球樹、M樹、VP樹、MVP樹
2.7.一、球樹
我們來針對上文內容總結回顧下,針對下面這樣一棵kd樹:
現要找它的最近鄰。
經過上文2.5節,總結來講,咱們已經知道:
一、爲了找到一個給定目標點的最近鄰,須要從樹的根結點開始向下沿樹找出目標點所在的區域,以下圖所示,給定目標點,用星號標示,咱們彷佛一眼看出,有一個點離目標點最近,由於它落在以目標點爲圓心以較小長度爲半徑的虛線圓內,但爲了肯定是否可能還村莊一個最近的近鄰,咱們會先檢查葉節點的同胞結點,然葉節點的同胞結點在圖中所示的陰影部分,虛線圓並不與之相交,因此肯定同胞葉結點不可能包含更近的近鄰。
二、因而咱們回溯到父節點,並檢查父節點的同胞結點,父節點的同胞結點覆蓋了圖中全部橫線X軸上的區域。由於虛線圓與右上方的矩形(KD樹把二維平面劃分紅一個一個矩形)相交...
如上,咱們看到,KD樹是可用於有效尋找最近鄰的一個樹結構,但這個樹結構其實並不完美,當處理不均勻分佈的數據集時便會呈現出一個基本衝突:既邀請樹有完美的平衡結構,又要求待查找的區域近似方形,但無論是近似方形,仍是矩形,甚至正方形,都不是最好的使用形狀,由於他們都有角。
什麼意思呢?就是說,在上圖中,若是黑色的實例點離目標點星點再遠一點,那麼勢必那個虛線圓會如紅線所示那樣擴大,以至與左上方矩形的右下角相交,既然相交了,那麼勢必又必須檢查這個左上方矩形,而實際上,最近的點離星點的距離很近,檢查左上方矩形區域已經是多餘。於此咱們看見,KD樹把二維平面劃分紅一個一個矩形,但矩形區域的角倒是個難以處理的問題。
解決的方案就是使用以下圖所示的球樹:
先從球中選擇一個離球的中心最遠的點,而後選擇第二個點離第一個點最遠,將球中全部的點分配到離這兩個聚類中心最近的一個上,而後計算每一個聚類的中心,以及聚類可以包含它全部數據點所需的最小半徑。這種方法的優勢是分裂一個包含n個殊絕點的球的成本只是隨n呈線性增長。
使用球樹找出給定目標點的最近鄰方法是,首先自上而下貫穿整棵樹找出包含目標點所在的葉子,並在這個球裏找出與目標點最靠近的點,這將肯定出目標點距離它的最近鄰點的一個上限值,而後跟KD樹查找同樣,檢查同胞結點,若是目標點到同胞結點中心的距離超過同胞結點的半徑與當前的上限值之和,那麼同胞結點裏不可能存在一個更近的點;不然的話,必須進一步檢查位於同胞結點如下的子樹。
以下圖,目標點仍是用一個星表示,黑色點是當前已知的的目標點的最近鄰,灰色球裏的全部內容將被排除,由於灰色球的中心點離的太遠,因此它不可能包含一個更近的點,像這樣,遞歸的向樹的根結點進行回溯處理,檢查全部可能包含一個更近於當前上限值的點的球。
球樹是自上而下的創建,和KD樹同樣,根本問題就是要找到一個好的方法將包含數據點集的球分裂成兩個,在實踐中,沒必要等到葉子結點只有兩個胡數據點時才中止,能夠採用和KD樹同樣的方法,一旦結點上的數據點打到預先設置的最小數量時,即可提早中止建樹過程。
也就是上面所述,先從球中選擇一個離球的中心最遠的點,而後選擇第二個點離第一個點最遠,將球中全部的點分配到離這兩個聚類中心最近的一個上,而後計算每一個聚類的中心,以及聚類可以包含它全部數據點所需的最小半徑。這種方法的優勢是分裂一個包含n個殊絕點的球的成本只是隨n呈線性增長(注:本小節內容主要來自參考條目19:數據挖掘實用機器學習技術,[新西蘭]Ian H.Witten 著,第4章4.7節)。
2.7.二、VP樹與MVP樹簡介
高維特徵向量的距離索引問題是基於內容的圖像檢索的一項關鍵技術,目前常常採用的解決辦法是首先對高維特徵空間作降維處理,而後採用包括四叉樹、kd樹、R樹族等在內的主流多維索引結構,這種方法的出發點是:目前的主流多維索引結構在處理維數較低的狀況時具備比較好的效率,但對於維數很高的狀況則顯得力不從心(即所謂的維數危機) 。
實驗結果代表當特徵空間的維數超過20 的時候,效率明顯下降,而可視化特徵每每採用高維向量描述,通常狀況下能夠達到10^2的量級,甚至更高。在表示圖像可視化特徵的高維向量中各維信息的重要程度是不一樣的,經過降維技術去除屬於次要信息的特徵向量以及相關性較強的特徵向量,從而下降特徵空間的維數,這種方法已經獲得了一些實際應用。
然而這種方法存在不足之處採用降維技術可能會致使有效信息的損失,尤爲不適合於處理特徵空間中的特徵向量相關性很小的狀況。另外主流的多維索引結構大都針對歐氏空間,設計須要利用到歐氏空間的幾何性質,而圖像的類似性計算極可能不限於基於歐氏距離。這種狀況下人們愈來愈關注基於距離的度量空間高維索引結構能夠直接應用於高維向量類似性查詢問題。
度量空間中對象之間的距離度量只能利用三角不等式性質,而不能利用其餘幾何性質。向量空間能夠看做由實數座標串組成的特殊度量空間,目前針對度量空間的高維索引問題提出的索引結構有不少種大體能夠做以下分類,以下圖所示:
其中,VP樹和MVP樹中特徵向量的舉例表示爲:
讀者點評:
- UESTC_HN_AY_GUOBO:如今主要是在kdtree的基礎上有了mtree或者mvptree,其實關鍵仍是pivot的選擇,以及度量空間中算法怎麼減小距離計算;
- mandycool:mvp-tree,是利用三角形不等式來縮小搜索區域的,不過mvp-tree的目標稍有不一樣,查詢的是到query點的距離小於某個值r的點;另外做者test的數據集只有20維,不知道上百維之後效果如何,而減小距離計算的一個思路是作embedding,經過不等式排除掉一部分點。
更多內容請參見論文1:DIST ANCE-BASED INDEXING FOR HIGH-DIMENSIONAL METRIC SP ACES,做者:Tolga Bozkaya & Meral Ozsoyoglu,及論文2:基於度量空間高維索引結構VP-tree及MVP-tree的圖像檢索,王志強,甘國輝,程起敏。
固然,若是你以爲上述論文還不夠知足你胃口的話,這裏有一大堆nearest neighbor algorithms相關的論文可供你看:http://scholar.google.com.hk/scholar?q=nearest+neighbor+algorithms&btnG=&hl=zh-CN&as_sdt=0&as_vis=1(其中,這篇能夠看下:Spill-Trees,An investigation of practical approximate nearest neighbor algorithms)。
第三部分、KD樹的應用:SIFT+KD_BBF搜索算法
3.一、SIFT特徵匹配算法
以前本blog內闡述過圖像特徵匹配SIFT算法,寫過五篇文章,這五篇文章分別爲:
不熟悉SIFT算法相關概念的能夠看上述幾篇文章,這裏再也不作贅述。與此同時,本文此部分也做爲十五個經典算法研究系列裏SIFT算法的九之四續。
OK,咱們知道,在sift算法中,給定兩幅圖片圖片,若要作特徵匹配,通常會先提取出圖片中的下列相關屬性做爲特徵點:
-
-
-
-
-
-
- struct feature
- {
- double x;
- double y;
- double a;
- double b;
- double c;
- double scl;
- double ori;
- int d;
- double descr[FEATURE_MAX_D];
- int type;
- int category;
- struct feature* fwd_match;
- struct feature* bck_match;
- struct feature* mdl_match;
- CvPoint2D64f img_pt;
- CvPoint2D64f mdl_pt;
- void* feature_data;
- char dense;
- };
然後在sift.h文件中定義兩個關鍵函數,這裏,咱們把它們稱之爲函數一,和函數二,以下所示:
函數一的聲明:
- extern int sift_features( IplImage* img, struct feature** feat );
函數一的實現:
- int sift_features( IplImage* img, struct feature** feat )
- {
- return _sift_features( img, feat, SIFT_INTVLS, SIFT_SIGMA, SIFT_CONTR_THR,
- SIFT_CURV_THR, SIFT_IMG_DBL, SIFT_DESCR_WIDTH,
- SIFT_DESCR_HIST_BINS );
- }
從上述函數一的實現中,咱們能夠看到,它內部實際上調用的是這個函數:_sift_features(..),也就是下面立刻要分析的函數二。
函數二的聲明:
- extern int _sift_features( IplImage* img, struct feature** feat, int intvls,
- double sigma, double contr_thr, int curv_thr,
- int img_dbl, int descr_width, int descr_hist_bins );
函數二的實現:
- int _sift_features( IplImage* img, struct feature** feat, int intvls,
- double sigma, double contr_thr, int curv_thr,
- int img_dbl, int descr_width, int descr_hist_bins )
- {
- IplImage* init_img;
- IplImage*** gauss_pyr, *** dog_pyr;
- CvMemStorage* storage;
- CvSeq* features;
- int octvs, i, n = 0,n0 = 0,n1 = 0,n2 = 0,n3 = 0,n4 = 0;
- int start;
-
-
- if( ! img )
- fatal_error( "NULL pointer error, %s, line %d", __FILE__, __LINE__ );
-
- if( ! feat )
- fatal_error( "NULL pointer error, %s, line %d", __FILE__, __LINE__ );
-
-
- start=GetTickCount();
- init_img = create_init_img( img, img_dbl, sigma );
- octvs = log( (float)(MIN( init_img->width, init_img->height )) ) / log((float)(2.0)) -5;
- gauss_pyr = build_gauss_pyr( init_img, octvs, intvls, sigma );
- dog_pyr = build_dog_pyr( gauss_pyr, octvs, intvls );
- fprintf( stderr, " creat the pyramid use %d\n",GetTickCount()-start);
-
- storage = cvCreateMemStorage( 0 );
- start=GetTickCount();
- features = scale_space_extrema( dog_pyr, octvs, intvls, contr_thr,
- curv_thr, storage );
- fprintf( stderr, " find the extrum points in DOG use %d\n",GetTickCount()-start);
-
- calc_feature_scales( features, sigma, intvls );
-
- if( img_dbl )
- adjust_for_img_dbl( features );
- start=GetTickCount();
- calc_feature_oris( features, gauss_pyr );
- fprintf( stderr, " get the main oritation use %d\n",GetTickCount()-start);
-
- start=GetTickCount();
- compute_descriptors( features, gauss_pyr, descr_width, descr_hist_bins );
- fprintf( stderr, " compute the descriptors use %d\n",GetTickCount()-start);
-
-
-
-
- n = features->total;
- *feat = (feature*)(calloc( n, sizeof(struct feature) ));
- *feat = (feature*)(cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ ));
-
- for( i = 0; i < n; i++ )
- {
- free( (*feat)[i].feature_data );
- (*feat)[i].feature_data = NULL;
- if((*feat)[i].dense == 4) ++n4;
- else if((*feat)[i].dense == 3) ++n3;
- else if((*feat)[i].dense == 2) ++n2;
- else if((*feat)[i].dense == 1) ++n1;
- else ++n0;
- }
-
-
-
- fprintf( stderr, "In the total feature points the extent4 points is %d\n",n4);
- fprintf( stderr, "In the total feature points the extent3 points is %d\n",n3);
- fprintf( stderr, "In the total feature points the extent2 points is %d\n",n2);
- fprintf( stderr, "In the total feature points the extent1 points is %d\n",n1);
- fprintf( stderr, "In the total feature points the extent0 points is %d\n",n0);
- cvReleaseMemStorage( &storage );
- cvReleaseImage( &init_img );
- release_pyr( &gauss_pyr, octvs, intvls + 3 );
- release_pyr( &dog_pyr, octvs, intvls + 2 );
-
- return n;
- }
說明:上面的函數二,包含了SIFT算法中幾乎全部函數,是SIFT算法的核心。本文不打算一一分析上面全部函數,只會抽取其中涉及到BBF查詢機制相關的函數。
3.二、K個最小近鄰的查找:大頂堆優先級隊列
上文中一直在講最近鄰問題,也就是說只找最近的那惟一一個鄰居,但若是現實中須要咱們找到k個最近的鄰居。該如何作呢?對的,以前blog內曾相近闡述過尋找最小的k個數的問題,顯然,尋找k個最近鄰與尋找最小的k個數的問題一模一樣。
回憶下尋找k個最小的數中關於構造大頂堆的解決方案:
「尋找最小的k個樹,更好的辦法是維護k個元素的最大堆,即用容量爲k的最大堆存儲最早遍歷到的k個數,並假設它們便是最小的k個數,建堆費時O(k)後,有k1<k2<...<kmax(kmax設爲大頂堆中最大元素)。繼續遍歷數列,每次遍歷一個元素x,與堆頂元素比較,x<kmax,更新堆(用時logk),不然不更新堆。這樣下來,總費時O(k+(n-k)*logk)=O(n*logk)。此方法得益於在堆中,查找等各項操做時間複雜度均爲logk。」
根據上述方法,我們來寫大頂堆最小優先級隊列相關代碼實現:
- void* minpq_extract_min( struct min_pq* min_pq )
- {
- void* data;
-
- if( min_pq->n < 1 )
- {
- fprintf( stderr, "Warning: PQ empty, %s line %d\n", __FILE__, __LINE__ );
- return NULL;
- }
- data = min_pq->pq_array[0].data;
- min_pq->n--;
- min_pq->pq_array[0] = min_pq->pq_array[min_pq->n];
- restore_minpq_order( min_pq->pq_array, 0, min_pq->n );
-
- return data;
- }
上述函數中,restore_minpq_order的實現以下:
- static void restore_minpq_order( struct pq_node* pq_array, int i, int n )
- {
- struct pq_node tmp;
- int l, r, min = i;
-
- l = left( i );
- r = right( i );
- if( l < n )
- if( pq_array[l].key < pq_array[i].key )
- min = l;
- if( r < n )
- if( pq_array[r].key < pq_array[min].key )
- min = r;
-
- if( min != i )
- {
- tmp = pq_array[min];
- pq_array[min] = pq_array[i];
- pq_array[i] = tmp;
- restore_minpq_order( pq_array, min, n );
- }
- }
3.三、KD樹近鄰搜索改進之BBF算法
原理在上文第二部分已經闡述明白,結合大頂堆找最近的K個近鄰思想,相關主體代碼以下:
-
- int kdtree_bbf_knn( struct kd_node* kd_root, struct feature* feat, int k,
- struct feature*** nbrs, int max_nn_chks )
- {
- struct kd_node* expl;
- struct min_pq* min_pq;
- struct feature* tree_feat, ** _nbrs;
- struct bbf_data* bbf_data;
- int i, t = 0, n = 0;
-
- if( ! nbrs || ! feat || ! kd_root )
- {
- fprintf( stderr, "Warning: NULL pointer error, %s, line %d\n",
- __FILE__, __LINE__ );
- return -1;
- }
-
- _nbrs = (feature**)(calloc( k, sizeof( struct feature* ) ));
- min_pq = minpq_init();
- minpq_insert( min_pq, kd_root, 0 );
-
-
- while( min_pq->n > 0 && t < max_nn_chks )
- {
-
-
-
- expl = (struct kd_node*)minpq_extract_min( min_pq );
- if( ! expl )
- {
- fprintf( stderr, "Warning: PQ unexpectedly empty, %s line %d\n",
- __FILE__, __LINE__ );
- goto fail;
- }
-
-
-
-
-
- expl = explore_to_leaf( expl, feat, min_pq );
- if( ! expl )
- {
- fprintf( stderr, "Warning: PQ unexpectedly empty, %s line %d\n",
- __FILE__, __LINE__ );
- goto fail;
- }
-
- for( i = 0; i < expl->n; i++ )
- {
-
- tree_feat = &expl->features[i];
- bbf_data = (struct bbf_data*)(malloc( sizeof( struct bbf_data ) ));
- if( ! bbf_data )
- {
- fprintf( stderr, "Warning: unable to allocate memory,"
- " %s line %d\n", __FILE__, __LINE__ );
- goto fail;
- }
- bbf_data->old_data = tree_feat->feature_data;
- bbf_data->d = descr_dist_sq(feat, tree_feat);
- tree_feat->feature_data = bbf_data;
-
-
- n += insert_into_nbr_array( tree_feat, _nbrs, n, k );
- }
- t++;
- }
-
- minpq_release( &min_pq );
- for( i = 0; i < n; i++ )
- {
- bbf_data = (struct bbf_data*)(_nbrs[i]->feature_data);
- _nbrs[i]->feature_data = bbf_data->old_data;
- free( bbf_data );
- }
- *nbrs = _nbrs;
- return n;
-
- fail:
- minpq_release( &min_pq );
- for( i = 0; i < n; i++ )
- {
- bbf_data = (struct bbf_data*)(_nbrs[i]->feature_data);
- _nbrs[i]->feature_data = bbf_data->old_data;
- free( bbf_data );
- }
- free( _nbrs );
- *nbrs = NULL;
- return -1;
- }
依據上述函數kdtree_bbf_knn從上而下看下來,注意幾點:
一、上述函數kdtree_bbf_knn中,explore_to_leaf的代碼以下:
- static struct kd_node* explore_to_leaf( struct kd_node* kd_node, struct feature* feat,
- struct min_pq* min_pq )
- {
- struct kd_node* unexpl, * expl = kd_node;
- double kv;
- int ki;
-
- while( expl && ! expl->leaf )
- {
- ki = expl->ki;
- kv = expl->kv;
-
- if( ki >= feat->d )
- {
- fprintf( stderr, "Warning: comparing imcompatible descriptors, %s" \
- " line %d\n", __FILE__, __LINE__ );
- return NULL;
- }
- if( feat->descr[ki] <= kv )
- {
- unexpl = expl->kd_right;
- expl = expl->kd_left;
- }
- else
- {
- unexpl = expl->kd_left;
- expl = expl->kd_right;
- }
-
-
- if( minpq_insert( min_pq, unexpl, ABS( kv - feat->descr[ki] ) ) )
- {
- fprintf( stderr, "Warning: unable to insert into PQ, %s, line %d\n",
- __FILE__, __LINE__ );
- return NULL;
- }
- }
- return expl;
- }
二、上述查找函數kdtree_bbf_knn中的參數k可調,表明的是要查找近鄰的個數,即number of neighbors to find,在sift特徵匹配中,k通常取2
- k = kdtree_bbf_knn( kd_root_0, feat, 2, &nbrs, KDTREE_BBF_MAX_NN_CHKS_0 );
- if( k == 2 )
三、上述函數kdtree_bbf_knn中「bbf_data->d = descr_dist_sq(feat, tree_feat); //計算兩個關鍵點描述器差平方和」,使用的計算方法是本文第一部分1.2節中所述的歐氏距離。
3.三、SIFT+BBF算法匹配效果
以前試了下sift + KD + BBF算法,用兩幅不一樣的圖片作了下匹配(固然,運行結果顯示是不匹配的),效果還不錯:http://weibo.com/1580904460/yDmzAEwcV#1348475194313。
「編譯的過程當中,直接用的VS2010 + opencv(並沒下gsl)。2012.09.24」。....
本文完整源碼有pudn賬號的朋友能夠前去這裏下載:http://www.pudn.com/downloads340/sourcecode/graph/texture_mapping/detail1486667.html(沒有pudn帳號的同窗請加羣:169056165,至羣共享下載,驗證信息:sift)。感謝諸位。
參考文獻及推薦閱讀
- 維基百科,http://en.wikipedia.org/wiki/K-nearest_neighbor_algorithm;
- 機器學習中的類似性度量,http://www.cnblogs.com/heaad/archive/2011/03/08/1977733.html;
- 傑卡德類似係數及距離,http://www.cnblogs.com/heaad/archive/2011/03/08/1977733.html;
- 統計學習方法,李航;
- 機率論與數理統計 第四版 盛驟等編,高教版;
- 《圖像局部不變特性特徵與描述》王永明 王貴錦 編著;
- 數據挖掘:實用機器學習技術,[新西蘭]Ian H.Witten 著,第4章4.7節;
- 模式分類,第4章 非參數技術,[美] IRichard O. Duda / Peter E. Hart / David G. Stork 著;
- http://underthehood.blog.51cto.com/2531780/687160;
- http://grunt1223.iteye.com/blog/921371;
- http://www.cnblogs.com/eyeszjwang/articles/2429382.html;
- http://blog.csdn.net/ijuliet/article/details/4471311;
- Rob Hess維護的sift庫,http://blogs.oregonstate.edu/hess/code/sift/;
- 酷殼,http://coolshell.cn/articles/8052.html;
- rubyist,http://segmentfault.com/q/1010000000094674;
- 皮爾遜相關係數維基百科頁面,http://t.cn/zjy6Gpg;
- 皮爾遜相關係數的一個應用:http://www.sobuhu.com/archives/567;
- http://blog.csdn.net/wsywl/article/details/5727327;
- 標準差,http://zh.wikipedia.org/wiki/%E6%A0%87%E5%87%86%E5%B7%AE;
- 協方差與相關性,http://t.cn/zjyXFRB;
- 電子科大kd樹電子課件:http://t.cn/zjbpXna;
- 編程藝術之尋找最小的k個數:http://blog.csdn.net/v_JULY_v/article/details/6403777;
- 機器學習那些事兒,http://vdisk.weibo.com/s/ix_9F;
- 大嘴巴漫談數據挖掘,http://vdisk.weibo.com/s/bUbzJ;
- http://www.codeproject.com/Articles/18113/KD-Tree-Searching-in-N-dimensions-Part-I;
- 一個庫:http://docs.pointclouds.org/trunk/group__kdtree.html;
- 3D上使用kd樹:http://pointclouds.org/;
- 編輯數學公式:http://webdemo.visionobjects.com/equation.html?locale=zh_CN;
- 基於R樹的最近鄰查找:http://blog.sina.com.cn/s/blog_72e1c7550101dsc3.html;
- 包含一個demo:http://www.leexiang.com/kd-tree;
- 機器學習相關降維算法,http://www.cnblogs.com/xbinworld/category/337861.html;
- Machine Learning相關topic,http://www.cnblogs.com/jerrylead/tag/Machine%20Learning/;
- 機器學習中的數學,http://www.cnblogs.com/LeftNotEasy/category/273623.html;
- 一堆概念性wikipedia頁面;
- 基於度量空間高維索引結構VP-tree及MVP-tree的圖像檢索,王志強,甘國輝,程起敏;
- Spill-Trees,An investigation of practical approximate nearest neighbor algorithms;
- DIST ANCE-BASED INDEXING FOR HIGH-DIMENSIONAL METRIC SP ACES,做者:Tolga Bozkaya & Meral Ozsoyoglu;
- 「Multidimensional Binary Search Trees Used for Associative Searching」,Jon Louis Bentley。
後記
從當天下午4點多一直寫,一直寫,直接寫到了隔日凌晨零點左右,或參考或直接間接引用了不少人的做品及一堆
wikipedia頁面(
固然,已經註明在上面參考文獻及推薦閱讀中前半部分條目,後半部分條目則爲行文以後看到的一些好的資料文章,能夠課外讀讀),但本文仍是有諸多地方有待修補完善,也歡迎廣大讀者不吝賜教 & 指正,感謝你們。完。July、二零一二年十一月二十一日零點五分。.
預告
本文以後,待寫的幾篇文章羅列以下:
- 數據挖掘中所需的機率論與數理統計知識,估計年末前完成(更新:已經完成,點擊左邊文字連接便可打開);
- 機器學習中相關的降維方法,如PCA/LDA等等,明年年初完成;
- 神經網絡入門學習導論,明年1月內完成;
- 程序員編程藝術第二十八章,時間待定;
- ...
歡迎期待。July、二零一二年十二月十七日。
(2):決策樹 Python代碼:http://blog.csdn.net/wyb_009/article/details/9191325
決策樹是個極其易懂的算法,建好模型後就是一連串嵌套的if..else...或嵌套的switch。
優勢:計算複雜度不高,輸出結果易於理解,對中間值的缺失不敏感,能夠處理不相關特徵數據;
缺點:可能會產生過分匹配的問題;
適用數據類型:數值型和標稱型。
決策樹的Python實現:
(一)先實現幾個工具函數:計算熵函數,劃分數據集工具函數,計算最大機率屬性;
(1)計算熵:熵表明集合的無序程度,集合越無序,熵越大;
- def entropy(dataset):
- from math import log
- log2 = lambda x:log(x)/log(2)
-
- results={}
- for row in dataset:
- r = row[len(row)-1]
- results[r] = results.get(r, 0) + 1
-
- ent = 0.0
- for r in results.keys():
- p = float(results[r]) / len(dataset)
- ent=ent-p*log2(p)
- return ent
-
(2)按屬性和值獲取數據集:
- def fetch_subdataset(dataset, k, v):
- return [d[:k]+d[k+1:] for d in dataset if d[k] == v]
這個函數只有短短一行,他的意義是:從dataset序列中取得第k列的值爲v的子集,並從得到的子集中去掉第k列。python的簡單優美顯現無遺。
(3)計算最大機率屬性。在構建決策樹時,在處理全部決策屬性後,還不能惟一區分數據時,咱們採用多數表決的方法來選擇最終分類:
- def get_max_feature(class_list):
- class_count = {}
- for cla in class_list:
- class_count[cla] = class_count.get(cla, 0) + 1
- sorted_class_count = sorted(class_count.items(), key=lambda d: d[1], reverse=True)
- return sorted_class_count[0][0]
(二)選取最優數據劃分方式函數:
選擇集合的最優劃分方式:以哪一列的值劃分集合,從而得到最大的信息增益呢?
- def choose_decision_feature(dataset):
- ent, feature = 100000000, -1
- for i in range(len(dataset[0]) - 1):
- feat_list = [e[i] for e in dataset]
- unq_feat_list = set(feat_list)
- ent_t = 0.0
- for f in unq_feat_list:
- sub_data = fetch_subdataset(dataset, i, f)
- ent_t += entropy(sub_data) * len(sub_data) / len(dataset)
-
- if ent_t < ent:
- ent, feature = ent_t, i
-
- return feature
(三)遞歸構建決策樹:
- def build_decision_tree(dataset, datalabel):
- cla = [c[-1] for c in dataset]
- if len(cla) == cla.count(cla[0]):
- return cla[0]
- if len(dataset[0]) == 1:
- return get_max_feature(dataset)
-
- feature = choose_decision_feature(dataset)
- feature_label = datalabel[feature]
- decision_tree = {feature_label:{}}
- del(datalabel[feature])
-
- feat_value = [d[feature] for d in dataset]
- unique_feat_value = set(feat_value)
- for value in unique_feat_value:
- sub_label = datalabel[:]
- decision_tree[feature_label][value] = build_decision_tree(\
- fetch_subdataset(dataset, feature, value), sub_label)
-
- return decision_tree
(四)使用決策樹
- def classify(decision_tree, feat_labels, testVec):
- label = decision_tree.keys()[0]
- next_dict = decision_tree[label]
- feat_index = feat_labels.index(label)
- for key in next_dict.keys():
- if testVec[feat_index] == key:
- if type(next_dict[key]).__name__ == 'dict':
- c_label = classify(next_dict[key], feat_labels, testVec)
- else:
- c_label = next_dict[key]
- return c_label
(五)決策樹持久化
(1)存儲
- def store_decision_tree(tree, filename):
- import pickle
- f = open(filename, 'w')
- pickle.dump(tree, f)
- f.close()
(2)讀取
- def load_decision_tree(filename):
- import pickle
- f = open(filename)
- return pickle.load(f)
(六)到了最後了,該回到主題了,給眼鏡男配眼鏡了。
下面的隱形眼鏡數據集來自UCI數據庫,它包含不少患者眼部情況的觀察條件以及醫生推薦的隱形眼鏡類型,隱形眼鏡類型包括硬材料、軟材料和不適合佩戴隱形眼鏡。
數據以下:
- young myope no reduced no lenses
- young myope no normal soft
- young myope yes reduced no lenses
- young myope yes normal hard
- young hyper no reduced no lenses
- young hyper no normal soft
- young hyper yes reduced no lenses
- young hyper yes normal hard
- pre myope no reduced no lenses
- pre myope no normal soft
- pre myope yes reduced no lenses
- pre myope yes normal hard
- pre hyper no reduced no lenses
- pre hyper no normal soft
- pre hyper yes reduced no lenses
- pre hyper yes normal no lenses
- presbyopic myope no reduced no lenses
- presbyopic myope no normal no lenses
- presbyopic myope yes reduced no lenses
- presbyopic myope yes normal hard
- presbyopic hyper no reduced no lenses
- presbyopic hyper no normal soft
- presbyopic hyper yes reduced no lenses
- presbyopic hyper yes normal no lenses
測試程序以下:
- def test():
- f = open('lenses.txt')
- lense_data = [inst.strip().split('\t') for inst in f.readlines()]
- lense_label = ['age', 'prescript', 'astigmatic', 'tearRate']
- lense_tree = build_decision_tree(lense_data, lense_label)
我這裏測試結果以下:
眼鏡男終於能夠買到合適的眼鏡啦。。。
全部代碼黏在下面:
- def entropy(dataset):
- from math import log
- log2 = lambda x:log(x)/log(2)
-
- results={}
- for row in dataset:
- r = row[len(row)-1]
- results[r] = results.get(r, 0) + 1
-
- ent = 0.0
- for r in results.keys():
- p = float(results[r]) / len(dataset)
- ent=ent-p*log2(p)
- return ent
-
- def fetch_subdataset(dataset, k, v):
- return [d[:k]+d[k+1:] for d in dataset if d[k] == v]
-
- def get_max_feature(class_list):
- class_count = {}
- for cla in class_list:
- class_count[cla] = class_count.get(cla, 0) + 1
- sorted_class_count = sorted(class_count.items(), key=lambda d: d[1], reverse=True)
- return sorted_class_count[0][0]
-
- def choose_decision_feature(dataset):
- ent, feature = 100000000, -1
- for i in range(len(dataset[0]) - 1):
- feat_list = [e[i] for e in dataset]
- unq_feat_list = set(feat_list)
- ent_t = 0.0
- for f in unq_feat_list:
- sub_data = fetch_subdataset(dataset, i, f)
- ent_t += entropy(sub_data) * len(sub_data) / len(dataset)
-
- if ent_t < ent:
- ent, feature = ent_t, i
-
- return feature
-
- def build_decision_tree(dataset, datalabel):
- cla = [c[-1] for c in dataset]
- if len(cla) == cla.count(cla[0]):
- return cla[0]
- if len(dataset[0]) == 1:
- return get_max_feature(dataset)
-
- feature = choose_decision_feature(dataset)
- feature_label = datalabel[feature]
- decision_tree = {feature_label:{}}
- del(datalabel[feature])
-
- feat_value = [d[feature] for d in dataset]
- unique_feat_value = set(feat_value)
- for value in unique_feat_value:
- sub_label = datalabel[:]
- decision_tree[feature_label][value] = build_decision_tree(\
- fetch_subdataset(dataset, feature, value), sub_label)
-
- return decision_tree
-
- def store_decision_tree(tree, filename):
- import pickle
- f = open(filename, 'w')
- pickle.dump(tree, f)
- f.close()
-
- def load_decision_tree(filename):
- import pickle
- f = open(filename)
- return pickle.load(f)
-
- def classify(decision_tree, feat_labels, testVec):
- label = decision_tree.keys()[0]
- next_dict = decision_tree[label]
- feat_index = feat_labels.index(label)
- for key in next_dict.keys():
- if testVec[feat_index] == key:
- if type(next_dict[key]).__name__ == 'dict':
- c_label = classify(next_dict[key], feat_labels, testVec)
- else:
- c_label = next_dict[key]
- return c_label
-
- def test():
- f = open('lenses.txt')
- lense_data = [inst.strip().split('\t') for inst in f.readlines()]
- lense_label = ['age', 'prescript', 'astigmatic', 'tearRate']
- lense_tree = build_decision_tree(lense_data, lense_label)
- return lense_tree
-
- if __name__ == "__main__":
- tree = test()
- print tree