从K近邻算法、距离度量谈到KD树、SIFT+BBF算法
前言
前两日,在微博上说:“到今天为止,我至少亏欠了3篇文章待写:1、KD树;2、神经网络;3、编程艺术第28章。你看到,blog内的文章与你于别处所见的任何都不同。于是,等啊等,等一台电脑,只好等待..”。得益于田,借了我一台电脑(借他电脑的时候,我连表示感谢,他说“能找到工作全靠你的博客,这点儿小忙还说,不地道”,有的时候,稍许感受到受人信任也是一种压力,愿我不辜负大家对我的信任),于是今天开始Top 10 Algorithms in Data Mining系列第三篇文章,即本文「从K近邻算法谈到KD树、SIFT+BBF算法」的创作。
一个人坚持自己的兴趣是比较难的,因为太多的人太容易为外界所动了,而尤其当你无法从中得到多少实际性的回报时,所幸,我能一直坚持下来。毕达哥拉斯学派有句名言:“万物皆数”,最近读完「微积分概念发展史」后也感受到了这一点。同时,从算法到数据挖掘、机器学习,再到数学,其中每一个领域任何一个细节都值得探索终生,或许,这就是“终生为学”的意思。
本文各部分内容分布如下:
- 第一部分讲K近邻算法,其中重点阐述了相关的距离度量表示法,
- 第二部分着重讲K近邻算法的实现--KD树,和KD树的插入,删除,最近邻查找等操作,及KD树的一系列相关改进(包括BBF,M树等);
- 第三部分讲KD树的应用:SIFT+kd_BBF搜索算法。
同时,你将看到,K近邻算法同本系列的前两篇文章所讲的决策树分类贝叶斯分类,及支持向量机SVM一样,也是用于解决分类问题的算法,
而本数据挖掘十大算法系列也会按照分类,聚类,关联分析,预测回归等问题依次展开阐述。
OK,行文仓促,本文若有任何漏洞,问题或者错误,欢迎朋友们随时不吝指正,各位的批评也是我继续写下去的动力之一。感谢。
第一部分、K近邻算法
1.1、什么是K近邻算法
何谓K近邻算法,即K-Nearest Neighbor algorithm,简称KNN算法,单从名字来猜想,可以简单粗暴的认为是:K个最近的邻居,当K=1时,算法便成了最近邻算法,即寻找最近的那个邻居。为何要找邻居?打个比方来说,假设你来到一个陌生的村庄,现在你要找到与你有着相似特征的人群融入他们,所谓入伙。
用官方的话来说,所谓K近邻算法,即是给定一个训练数据集,对新的输入实例,在训练数据集中找到与该实例最邻近的K个实例(也就是上面所说的K个邻居),这K个实例的多数属于某个类,就把该输入实例分类到这个类中。根据这个说法,咱们来看下引自维基百科上的一幅图:
如上图所示,有两类不同的样本数据,分别用蓝色的小正方形和红色的小三角形表示,而图正中间的那个绿色的圆所标示的数据则是待分类的数据。也就是说,现在,我们不知道中间那个绿色的数据是从属于哪一类(蓝色小正方形or红色小三角形),下面,我们就要解决这个问题:给这个绿色的圆分类。
我们常说,物以类聚,人以群分,判别一个人是一个什么样品质特征的人,常常可以从他/她身边的朋友入手,所谓观其友,而识其人。我们不是要判别上图中那个绿色的圆是属于哪一类数据么,好说,从它的邻居下手。但一次性看多少个邻居呢?从上图中,你还能看到:
- 如果K=3,绿色圆点的最近的3个邻居是2个红色小三角形和1个蓝色小正方形,少数从属于多数,基于统计的方法,判定绿色的这个待分类点属于红色的三角形一类。
- 如果K=5,绿色圆点的最近的5个邻居是2个红色三角形和3个蓝色的正方形,还是少数从属于多数,基于统计的方法,判定绿色的这个待分类点属于蓝色的正方形一类。
于此我们看到,当无法判定当前待分类点是从属于已知分类中的哪一类时,我们可以依据统计学的理论看它所处的位置特征,衡量它周围邻居的权重,而把它归为(或分配)到权重更大的那一类。这就是K近邻算法的核心思想。
1.2、近邻的距离度量表示法
上文第一节,我们看到,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)间的欧氏距离:
也可以用表示成向量运算的形式:
其上,二维平面上两点欧式距离,代码可以如下编写:
//unixfy:计算欧氏距离
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的曼哈顿距离为:,要注意的是,曼哈顿距离依赖座标系统的转度,而非系统在座标轴上的平移或映射。
(1)二维平面两点a(x1,y1)与b(x2,y2)间的曼哈顿距离(2)两个n维向量a(x11,x12,…,x1n)与 b(x21,x22,…,x2n)间的曼哈顿距离
- 3. 切比雪夫距离,若二个向量或二个点p 、and q,其座标分别为及,则两者之间的切比雪夫距离定义如下:,
(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)间的闵可夫斯基距离定义为:其中p是一个变参数。当p=1时,就是曼哈顿距离当p=2时,就是欧氏距离当p→∞时,就是切比雪夫距离根据变参数的不同,闵氏距离可以表示一类的距离。
-
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)间的标准化欧氏距离的公式:
-
6. 马氏距离(Mahalanobis Distance)
( 协方差矩阵中每个元素是各个矢量元素之间的协方差Cov(X,Y),Cov(X,Y) = E{ [X-E(X)] [Y-E(Y)]},其中E为数学期望)而其中向量Xi与Xj之间的马氏距离定义为:
(2)马氏距离的优缺点:量纲无关,排除变量之间的相关性的干扰。「 微博上的seafood高清版点评道:原来马氏距离是根据协方差矩阵演变,一直被老师误导了,怪不得看Killian在05年NIPS发表的LMNN论文时候老是看到协方差矩阵和半正定,原来是这回事」 - 7、巴氏距离(Bhattacharyya Distance),在统计中,Bhattacharyya距离测量两个离散或连续概率分布的相似性。它与衡量两个统计样品或种群之间的重叠量的Bhattacharyya系数密切相关。Bhattacharyya距离和Bhattacharyya系数以20世纪30年代曾在印度统计研究所工作的一个统计学家A. Bhattacharya命名。同时,Bhattacharyya系数可以被用来确定两个样本被认为相对接近的,它是用来测量中的类分类的可分离性。
(1)巴氏距离的定义对于离散概率分布 p和q在同一域 X,它被定义为:
其中:
是Bhattacharyya系数。对于连续概率分布,Bhattacharyya系数被定义为:
在 这两种情况下,巴氏距离 并没有服从三角不等式.(值得一提的是,Hellinger距离不服从三角不等式 )。对于多变量的高斯分布 ,
,
和是手段和协方差的分布 。需要注意的是,在这种情况下,第一项中的Bhattacharyya距离与马氏距离有关联。(2)Bhattacharyya系数Bhattacharyya系数是两个统计样本之间的重叠量的近似测量,可以被用于确定被考虑的两个样本的相对接近。计算Bhattacharyya系数涉及集成的基本形式的两个样本的重叠的时间间隔的值的两个样本被分裂成一个选定的分区数,并且在每个分区中的每个样品的成员的数量,在下面的公式中使用
考虑样品a 和 b ,n是的分区数,并且 , 被一个 和 b i的日分区中的样本数量的成员。更多介绍请参看: http://en.wikipedia.org/wiki/Bhattacharyya_coefficient。
- 8. 汉明距离(Hamming distance), 两个等长字符串s1与s2之间的汉明距离定义为将其中一个变为另外一个所需要作的最小替换次数。例如字符串“1111”与“1001”之间的汉明距离为2。应用:信息编码(为了增强容错性,应使得编码间的最小汉明距离尽可能大)。
或许,你还没明白我再说什么,不急,看下 上篇blog中第78题的第3小题整理的一道面试题目,便一目了然了。如下图所示:
//动态规划:
//f[i,j]表示s[0...i]与t[0...j]的最小编辑距离。
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) }
//分别表示:添加1个,删除1个,替换1个(相同就不用替换)。
与此同时, 面试官还可以继续问下去:那么,请问,如何设计一个比较两篇文章相似性的算法?( 这个问题的讨论可以看看这里: http://t.cn/zl82CAH,及这里关于simhash算法的介绍: http://www.cnblogs.com/linecong/archive/2010/08/28/simhash.html ),接下来,便引出了下文关于夹角余弦的讨论。(上篇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)
(1) 杰卡德相似系数两个集合A和B的交集元素在A,B的并集中所占的比例,称为两个集合的杰卡德相似系数,用符号J(A,B)表示。
杰卡德相似系数是衡量两个集合的相似度一种指标。(2) 杰卡德距离与杰卡德相似系数相反的概念是杰卡德距离(Jaccard distance)。杰卡德距离可用如下公式表示:
杰卡德距离用两个集合中不同元素占所有元素的比例来衡量两个集合的区分度。(3) 杰卡德相似系数与杰卡德距离的应用可将杰卡德相似系数用在衡量样本的相似度上。
举例:样本A与样本B是两个n维向量,而且所有维度的取值都是0或1,例如:A(0111)和B(1011)。我们将样本看成是一个集合,1表示集合包含该元素,0表示集合不包含该元素。M 11 :样本A与B都是1的维度的个数M0 1:样本A是0,样本B是1的维度的个数M10:样本A是1,样本B是0 的维度的个数M 00:样本A与B都是0的维度的个数依据上文给的杰卡德相似系数及杰卡德距离的相关定义,样本A与B的杰卡德相似系数J可以表示为:
这里M 11+ M 01+ M 10可理解为A与B的并集的元素个数,而M 11是A与B的交集的元素个数。而样本A与B的杰卡德距离表示为J':
- 11.皮尔逊系数(Pearson 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的相关系数,记为 )
- 当相关系数为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,那么两变量间的皮尔逊相关系数可通过以下公式计算:
- 公式一:
注:勿忘了上面说过,“皮尔逊相关系数定义为两个变量之间的协方差和标准差的商”,其中标准差的计算公式为:
(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.3、K值的选择
除了上述1.2节如何定义邻居的问题之外,还有一个选择多少个邻居,即K值定义为多大的问题。不要小看了这个K值选择问题,因为它对K近邻算法的结果会产生重大影响。如李航博士的一书「统计学习方法」上所说:
- 如果选择较小的K值,就相当于用较小的领域中的训练实例进行预测,“学习”近似误差会减小,只有与输入实例较近或相似的训练实例才会对预测结果起作用,与此同时带来的问题是“学习”的估计误差会增大,换句话说,K值的减小就意味着整体模型变得复杂,容易发生过拟合;
- 如果选择较大的K值,就相当于用较大领域中的训练实例进行预测,其优点是可以减少学习的估计误差,但缺点是学习的近似误差会增大。这时候,与输入实例较远(不相似的)训练实例也会对预测器作用,使预测发生错误,且K值的增大就意味着整体的模型变得简单。
- K=N,则完全不足取,因为此时无论输入实例是什么,都只是简单的预测它属于在训练实例中最多的累,模型过于简单,忽略了训练实例中大量有用信息。
第二部分、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.1、什么是KD树
Kd-树是K-dimension tree的缩写,是对数据点在k维空间(如二维(x,y),三维(x,y,z),k维(x1,y,z..))中划分的一种数据结构,主要应用于多维空间关键数据的搜索(如:范围搜索和最近邻搜索)。本质上说,Kd-树就是一种平衡二叉树。
首先必须搞清楚的是,k-d树是一种空间划分树,说白了,就是把整个空间划分为特定的几个部分,然后在特定空间的部分内进行相关搜索操作。想像一个三维(多维有点为难你的想象力了)空间,kd树按照一定的划分规则把这个三维空间划分了多个空间,如下图所示:
2.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)};
与此同时,经过对上面所示的空间划分之后,我们可以看出,点(7,2)可以为根结点,从根结点出发的两条红粗斜线指向的(5,4)和(9,6)则为根结点的左右子结点,而(2,3),(4,7)则为(5,4)的左右孩子(通过两条细红斜线相连),最后,(8,1)为(9,6)的左孩子(通过细红斜线相连)。如此,便形成了下面这样一棵k-d树:
k-d树的数据结构
针对上表给出的kd树的数据结构,转化成具体代码如下所示(注,本文以下代码分析基于Rob Hess维护的sift库):
/** a node in a k-d tree */
struct kd_node
{
int ki; /**< partition key index *///关键点直方图方差最大向量系列位置
double kv; /**< partition key value *///直方图方差最大向量系列中最中间模值
int leaf; /**< 1 if node is a leaf, 0 otherwise */
struct feature* features; /**< features at this node */
int n; /**< number of features */
struct kd_node* kd_left; /**< left child */
struct kd_node* kd_right; /**< right child */
};
也就是说,如之前所述,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 ); //n--number of features,initinalize root of tree.
expand_kd_node_subtree( kd_root ); //kd tree expand
return kd_root;
}
上面的涉及初始化操作的两个函数kd_node_init,及expand_kd_node_subtree代码分别如下所示:
static struct kd_node* kd_node_init( struct feature* features, int n )
{ //n--number of features
struct kd_node* kd_node;
kd_node = (struct kd_node*)(malloc( sizeof( struct kd_node ) ));
memset( kd_node, 0, sizeof( struct kd_node ) ); //0填充
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 )
{
/* base case: leaf node */
if( kd_node->n == 1 || kd_node->n == 0 )
{ //叶节点 //伪叶节点
kd_node->leaf = 1;
return;
}
assign_part_key( kd_node ); //get ki,kv
partition_features( kd_node ); //creat left and right children,特征点ki位置左树比右树模值小,kv作为分界模值
//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.3、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.4、KD树的删除
从K-D树中删除一个结点是代价很高的,很清楚删除子树的根受到子树中结点个数的限制。用TPL(T)表示树T总的路径长度。可看出树中子树大小的总和为TPL(T)+N。 以随机方式插入N个点形成树的TPL是O(N*log2N),这就意味着从一个随机形成的K-D树中删除一个随机选取的结点平均代价的上界是O(log2N) 。
2.5、KD树的最近邻搜索算法
现实生活中有许多问题需要在多维数据的快速分析和快速搜索,对于这个问题最常用的方法是所谓的kd树。在k-d树中进行数据的查找也是特征匹配的重要环节,其目的是检索在k-d树中与查询点距离最近的数据点。在一个N维的笛卡儿空间在两个点之间的距离是由下述公式确定:
2.5.1、k-d树查询算法的伪代码
k-d树查询算法的伪代码如下所示:
算法:k-d树最邻近查找
输入:Kd, //k-d tree类型
target //查询数据点
输出:nearest, //最邻近数据点
dist //最邻近数据点和查询点间的距离
1. If Kd为NULL,则设dist为infinite并返回
2. //进行二叉查找,生成搜索路径
Kd_point = &Kd; //Kd-point中保存k-d tree根节点地址
nearest = Kd_point -> Node-data; //初始化最近邻点
while(Kd_point)
push(Kd_point)到search_path中; //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取出一个节点指针; //从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; //如果target位于左子空间,就应进入右子空间
else
Kd_point = back_point -> left; //如果target位于右子空间,就应进入左子空间
将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、举例:查询点(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.3、举例:查询点(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.6、kd树近邻搜索算法的改进:BBF算法
咱们顺着上一节的思路,参考统计学习方法一书上的内容,再来总结下kd树的最近邻搜索算法:
输入:以构造的kd树,目标点x;输出:x 的最近邻
算法步骤如下:
- 在kd树种找出包含目标点x的叶结点:从根结点出发,递归地向下搜索kd树。若目标点x当前维的坐标小于切分点的坐标,则移动到左子结点,否则移动到右子结点,直到子结点为叶结点为止。
- 以此叶结点为“当前最近点”。
- 递归的向上回溯,在每个结点进行以下操作:
(a)如果该结点保存的实例点比当前最近点距离目标点更近,则更新“当前最近点”,也就是说以该实例点为“当前最近点”。
(b)当前最近点一定存在于该结点一个子结点对应的区域,检查子结点的父结点的另一子结点对应的区域是否有更近的点。具体做法是,检查另一子结点对应的区域是否以目标点位球心,以目标点与 http://blog.csdn.net/v_july_v/article/details/8203674/