局部敏感哈希matlab代码解读

很早就想写一篇关于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。

利用哈希函数将数据点映射到不同的桶中是一种保形映射，使得数据点 ii 和数据点 jj 在原始空间的相似度 ss 与映射后的在同一个桶的概率呈现正相关。之所以这么做，主要是避免exhausted search. 如果理想状态，每个桶中的元素数目大致相同，那么查询时的运算量将从原来的数据样本数目mm个降低到m/km/k个，其中kk为桶的数目。

由于输出是二值向量，设其长度为LL，每个哈希值其实对应着一个桶，理想情况下每个桶中都有数据，k=2Lk = 2^L。

从原理上来说，代码实现是很简单的，matlab的版本的代码可见http://ttic.uchicago.edu/~gregory/download.html

这其实是一个比较完整的工具包

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

入口函数lsh

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

第一个参数是使用的算法的类型，包括两种类型，分别是lsh和e2lsh

生成一个range的参数，得到的[0 0 ,…0; 255 255 ,….,255]这样的形式

range = processRange(d,range);

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

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

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

for j=1:l
% select random dimensions
I(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;
end

这里hash函数就是一个简单 阈值函数，将原始的400维的数据，随机选出k=24维，变为0到1，后文会有进一步说明。l为总共生成的哈希函数的数目，这里取值为20。

产生Is的变量的内容如下：



d是选择的维度下标，t是维度的阈值。

T = lshprep(type,Is,b);

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 table
T(j).Index = {};
T(j).verbose=1;

% set up secondary hash table for buckets
% max. index can be obtained by running lshhash on max. bucket
T(j).bhash = cell(lshhash(ones(1,k)*255),1); % lshhash是一个计算hash值的函数，将24维的二值向量映射为一个哈希值

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

T = lshins(T,x,ind);

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

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

它里面的主要采用了矩阵的比较，本质上就是用刚才生成的阈值函数做了一个二值化。

其中v是一个59500*24维的二值矩阵，每一行表示一个数据样本。

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

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

哈希向量

来简称这个二值向量。

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

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

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

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

对每一个单独的哈希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 buckets
if (isempty(ihash)) % nothing matches
isb = [];
else % may or may not match
isb = ihash(find(all(bsxfun(@eq,uniqBuck(b,:),T(j).buckets(ihash,:)),2)));
end

其中

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

是一种非常有效的写法，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 bucket
newBuckets=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);%该桶中存放的元素的下标
end

随后完成信息的更新

% 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 table
T(j).Index{isb}= newIndex;%重新为属于该桶的数据下标赋值
% update distinct element count
T(j).count = T(j).count + length(newIndex)-oldcount;
%新数目减去老数目为改变量，注意如果以前桶中有元素，是通过追加的方式添加上去的，在追加后再与T(j).B进行比较。作者这么做，就是为了保证桶中元素不会因为满了而倾向于保持老元素，新元素就加不进去了，所以先追加后然后再随机选择指定数目保留下来。当然这样做还是会造成桶中旧的元素更容易被扔掉这一情形。

运行分析

运行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

参数查看 lshstats

examine statistics of LSH data structure

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

例如；

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

输出为

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);

得到结果

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);

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});toc

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

时间已过 0.030697 秒。

下面来解析这个函数的实现

需要完成的任务是找到所有match这个query的tables。

步骤1 用哈希函数T(j)获取查询x0的映射的50维(维度为哈希函数中随机选定的位数的长度，即b)二值向量，由于加了128，所以范围是在[128,129]。

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

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

key = lshhash(buck);

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

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

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

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

步骤5

去除重复数据

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

步骤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);%返回小于指定距离的下标，基于iNN
iNN=iNN(sortind(ind));% 返回相似数据，这就完成了检索

else % random

rp=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;
end
end
end
iNN = choose;
end

end