局部敏感哈希算法

这篇文章介绍了局部敏感哈希算法，局部敏感哈希是非监督的哈希算法。 
算法的输入是实数域的特征向量，输出为一个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);

第一个参数是使用的算法的类型，包括两种类型，分别是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