局部敏感哈希LSH，即matlab代码

本文主要是介绍局部敏感哈希LSH，即matlab代码，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

转自：http://blog.csdn.net/dudubird90/article/details/50907641

很早就想写一篇关于LSH的文章，后来发现前辈们已经写好了，容我这里再推荐一下该文。
Locality Sensitive Hashing(LSH)之随机投影法
http://www.strongczq.com/2012/04/locality-sensitive-hashinglsh%E4%B9%8B%E9%9A%8F%E6%9C%BA%E6%8A%95%E5%BD%B1%E6%B3%95.html

个人总结：这篇文章介绍了局部敏感哈希算法，局部敏感哈希是非监督的哈希算法。
算法的输入是实数域的特征向量，输出为一个binary vector。
利用哈希函数将数据点映射到不同的桶中是一种保形映射，使得数据点 i 和数据点 j 在原始空间的相似度 s 与映射后的在同一个桶的概率呈现正相关。之所以这么做，主要是避免exhausted search. 如果理想状态，每个桶中的元素数目大致相同，那么查询时的运算量将从原来的数据样本数目 m 个降低到 m/k 个，其中 k 为桶的数目。
由于输出是二值向量，设其长度为 L ，每个哈希值其实对应着一个桶，理想情况下每个桶中都有数据， k=2L 。
从原理上来说，代码实现是很简单的，matlab的版本的代码可见http://ttic.uchicago.edu/~gregory/download.html
这其实是一个比较完整的工具包

本文主要做关键部分的代码解析。

入口函数lsh

T1=lsh('lsh',20,24,size(patches,1),patches,'range',255);1
1

第一个参数是使用的算法的类型，包括两种类型，分别是lsh和e2lsh
生成一个range的参数，得到的[0 0 ,…0; 255 255 ,….,255]这样的形式

range = processRange(d,range);1
1

这个函数是用来产生lsh函数的。

Is = lshfunc(type,l,k,d,varargin{:});1
1

l表示函数的个数，k表示一个函数中的位数，d表示数据的维度。

   for j=1:l% select random dimensionsI(j).d = include(unidrnd(length(include),1,k)); % 均匀分布的，随机选中k维% for each dimension select a threshold% hash key = [[ x(:,d)' >= t ]]t = unifrnd(0,1,1,k).*(range(2,I(j).d)-range(1,I(j).d)); %每一维都随机选中一个阈值位于0~255之间I(j).t = range(1,I(j).d)+t;I(j).k = k;end1
2
3
4
5
6
7
8
9
1
2
3
4
5
6
7
8
9

这里hash函数就是一个简单阈值函数，将原始的400维的数据，随机选出k=24维，变为0到1，后文会有进一步说明。l为总共生成的哈希函数的数目，这里取值为20。
产生Is的变量的内容如下：
这里写图片描述
d是选择的维度下标，t是维度的阈值。

T = lshprep(type,Is,b);1
1

T这个变量存储了哈希查找哈希值以及索引信息。

  T(j).type = type;T(j).Args = varargin;T(j).I = Is(j);T(j).B = B;T(j).count = 0;T(j).buckets = [];% prepare T's tableT(j).Index = {};T(j).verbose=1;% set up secondary hash table for buckets% max. index can be obtained by running lshhash on max. bucketT(j).bhash = cell(lshhash(ones(1,k)*255),1); % lshhash是一个计算hash值的函数，将24维的二值向量映射为一个哈希值1
2
3
4
5
6
7
8
9
10
11
12
13
1
2
3
4
5
6
7
8
9
10
11
12
13

随后的函数，将数据放入桶中，对T中变量进行赋值。

  T = lshins(T,x,ind);1
1

这个函数中有一些关键的处理，其中

  buck = findbucket(T(j).type,x,T(j).I);%这是一个将数据转化为二值向量的函数1
1

它里面的主要采用了矩阵的比较，本质上就是用刚才生成的阈值函数做了一个二值化。
其中v是一个59500*24维的二值矩阵，每一行表示一个数据样本。

 v = x(I.d,:)' <= repmat(I.t,size(x,2),1);v = uint8(v+128);1
2
1
2

但注意，输出的d维二值向量每一维并不是[0， 1]，而在区间[128 129]，这可能是要用于后文二次哈希的计算方便。为了后文方便说明，我们用哈希向量来简称这个二值向量。

这里一个桶buck对应着一个哈希向量，但是桶的数目非常多，直接来进行比较是很费时间的。

  [uniqBuck,ib,bID] = unique(buck,'rows');keys = lshhash(uniqBuck);%返回每个桶的哈希key值1
2
1
2

例如，对j=1这个哈希函数而言，总共有14615个不同的桶（新分配空间为14615*24），如果要查找一个桶就需要14615次比较非常费时。作者的优化方案是进行二次哈希，让多个哈希向量映射为一个整型的hash-key值，用lshhash函数完成此功能。

  % allocate space for new buckets -- possibly excessiveT(j).buckets=[T(j).buckets; zeros(length(ib),T(j).I.k,'uint8')];1
2
1
2

对每一个单独的哈希key值ib(b)

    % find which data go to bucket uniqBuck(b)thisBucket = find(bID==bID(ib(b)));% find out if this bucket already has anything% first, which bucket is it? 该hash函数T(j)下的，对应于哈希key值keys(b)的桶是否已经存在ihash = T(j).bhash{keys(b)}; % possible matching bucketsif (isempty(ihash)) % nothing matchesisb = [];else % may or may not matchisb = ihash(find(all(bsxfun(@eq,uniqBuck(b,:),T(j).buckets(ihash,:)),2)));end1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11

其中

      isb = ihash(find(all(bsxfun(@eq,uniqBuck(b,:),T(j).buckets(ihash,:)),2)));1
1

是一种非常有效的写法，bsxfun(@eq ,a,b)这种形式会得到两个向量之间的逐位比较，它matlab内部的实现是通过循环来实现的。通过all在水平方向上进行判别，就相当于比较两个向量是否相等。这一步是比较在T(j).bhash中存放的哈希向量中是否已经存在当前的获得的哈希向量，即是否已经记录了当前的桶，这样我们就可以分情况讨论是往这个桶里添加新的数据，还是要先创建一个桶再添加新的数据。

  if (~isempty(isb)) % 如果isb不为空，那么即该bucket已经存在% adding to an existing bucket.oldcount=length(T(j).Index{isb}); % # elements in the bucket prior% to addition 添加前桶中元素的数目，主要是方便统计newIndex = [T(j).Index{isb}  ind(thisBucket)];else% creating new bucketnewBuckets=newBuckets+1;oldcount=0;isb = oldBuckets+newBuckets;T(j).buckets(isb,:)=uniqBuck(b,:);%为什么用128 129表示T(j).bhash{keys(b)} = [T(j).bhash{keys(b)}; isb];%根据hash-key值来映射桶序号newIndex = ind(thisBucket);%该桶中存放的元素的下标end1
2
3
4
5
6
7
8
9
10
11
12
13
14
1
2
3
4
5
6
7
8
9
10
11
12
13
14

随后完成信息的更新

    % if there is a bound on bucket capacity, and the bucket is full,% keep a random subset of B elements (note: we do this rather than% simply skip the new elements since that could introduce bias% towards older elements.)% There is still a bias since older elements have more chances to get% thrown out.if (length(newIndex) > T(j).B)rp=randperm(length(newIndex));newIndex = newIndex(rp(1:T(j).B));% 如果超过的了桶的容量限制，那么随机选定T(j).B个数据end% ready to put this into the tableT(j).Index{isb}= newIndex;%重新为属于该桶的数据下标赋值% update distinct element countT(j).count = T(j).count + length(newIndex)-oldcount;%新数目减去老数目为改变量，注意如果以前桶中有元素，是通过追加的方式添加上去的，在追加后再与T(j).B进行比较。作者这么做，就是为了保证桶中元素不会因为满了而倾向于保持老元素，新元素就加不进去了，所以先追加后然后再随机选择指定数目保留下来。当然这样做还是会造成桶中旧的元素更容易被扔掉这一情形。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

运行分析

运行lsh函数会得到：

Table 5 adding 13852 buckets (now 13852)
Table 5: 59500 elements
12619 distinct buckets
Table 6 adding 12619 buckets (now 12619)
Table 6: 59500 elements
11936 distinct buckets
Table 7 adding 11936 buckets (now 11936)
Table 7: 59500 elements
15997 distinct buckets
1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9
10

参数查看 lshstats

examine statistics of LSH data structure

[mi,ma,me]=lshstats(T,B,xref,xtst,minNN)1
1

例如；

lshstats(T1(1:5),'test',patches,patches(:,1:1000),2);1
1

输出为
Table 1: 59500 in 13404 bkts, med 1, max 4288, avg 813.19
Table 2: 59500 in 12661 bkts, med 1, max 2646, avg 544.55
Table 3: 59500 in 16147 bkts, med 1, max 4057, avg 751.01
Table 4: 59500 in 11627 bkts, med 1, max 4989, avg 864.60
Table 5: 59500 in 13630 bkts, med 1, max 3528, avg 601.55

这表示table1有13404 个桶，平均容量是每个桶1个数据，最大容量为4288，期望容量为813.19

Running test…10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
# of comparisons: mean 980.14, max 8122, failures: 54

这里使用了5个哈希函数，它的含义是对前1000个样本进行查找，平均每次查找需要比较980个样本，但是同时失败次数为54次

如果增加哈希函数的数目，会得到不同的结果，根据参考文献中的分析，如果增加哈希函数的数目，那么会需要更长的查找时间，但是同时recall将会增加，例如这里我们用全部的20个哈希函数来做实验。

 lshstats(T1,'test',patches,patches(:,1:1000),2);1
1

得到结果
Running test…10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
# of comparisons: mean 2957.24, max 13120, failures: 2
可以发现平均查找所需的时间变长了，但是recall相应的变高的(几乎没有错误)。

lshlookup

下面是查找第50个样本，在这之前，首先增加二值向量的长度，即引用文献中的b的长度，这会减少平均每个桶中的元素数目

lshstats(T2(1:10),'test',patches,patches(:,1:1000),2);1
1

Table 1: 59500 in 33066 bkts, med 1, max 1829, avg 146.51
Table 2: 59500 in 34018 bkts, med 1, max 1638, avg 160.95
Table 3: 59500 in 34077 bkts, med 1, max 1386, avg 156.09
Table 4: 59500 in 35716 bkts, med 1, max 2813, avg 210.50
Table 5: 59500 in 34492 bkts, med 1, max 1470, avg 194.75
Table 6: 59500 in 34659 bkts, med 1, max 1543, avg 156.86
Table 7: 59500 in 33033 bkts, med 1, max 1232, avg 146.30
Table 8: 59500 in 33923 bkts, med 1, max 1955, avg 152.32
Table 9: 59500 in 34032 bkts, med 1, max 1718, avg 176.25
Table 10: 59500 in 32402 bkts, med 1, max 2862, avg 226.41

注意avg变小了

tic; [nnlsh,numcand]=lshlookup(patches(:,50),patches,T2,'k',11,'distfun','lpnorm','distargs',{1});toc1
1

算法运行结果结果实现检索一个数据所需的时间：

时间已过 0.030697 秒。

下面来解析这个函数的实现
需要完成的任务是找到所有match这个query的tables。
步骤1 用哈希函数T(j)获取查询x0的映射的50维(维度为哈希函数中随机选定的位数的长度，即b)二值向量，由于加了128，所以范围是在[128,129]。

  buck = findbucket(T(j).type,x0,T(j).I); 1
1

步骤2 将该向量转化成哈希key，这一步不是一一映射，而是多对一的映射，主要目的是为了提升向量的检索速度。

 key = lshhash(buck);1
1

步骤3 根据哈希key值获取所有的哈希向量，一个哈希key值对应着多个bucket

 ihash = T(j).bhash{key}; % possible matching buckets1
1

步骤4 进一步查找到该哈希向量，即找到对应的桶

 if (~isempty(ihash)) % nothing matchesb = ihash(find(all(bsxfun(@eq,buck,T(j).buckets(ihash,:)),2)));if (~isempty(b))iNN = [iNN T(j).Index{b}]; %把该桶中的数据union起来，因为不同的哈希函数会有不同的结果endend1
2
3
4
5
6
1
2
3
4
5
6

步骤5
去除重复数据

[iNN,iu]=unique(iNN);
cand = length(iNN);1
2
1
2

步骤6
这一步主要是将相似列表中的数据做个排序返回。用于CBIR检索很合适。

if (~isempty(iNN))if (strcmp(sel,'best'))D=feval(distfun,x0,Xsel(x,iNN),distargs{:});% 即比较这些桶中的最近邻数据和query的距离[dist,sortind]=sort(D);ind = find(dist(1:min(k,length(dist)))<=r);%返回小于指定距离的下标，基于iNNiNN=iNN(sortind(ind));% 返回相似数据，这就完成了检索else % randomrp=randperm(cand);choose=[];for i=1:length(rp)d = feval(distfun,x0,Xsel(x,iNN(rp(i))),distargs{:});if (d <= r) choose = [choose iNN(rp(i))];if (length(choose) == k)break;endendendiNN = choose;endend