利用K-Means LSH 加速求解格中的最短向量问题

金悦祺, 胡红钢

中国科学院 电磁空间信息重点实验室, 合肥230027

1 问题介绍

1.1 格密码

随着量子计算的深入研究, 越来越多的高性能量子算法被提出, 这使得某些传统的密码体系在未来的量子手段下不堪一击. 比如著名的Shor 算法就能够在多项式时间分解大整数, 从而攻破目前广泛应用的RSA 公钥体系. 基于格的密码体系具有对抗量子计算攻击的能力, 因此在近些年受到人们的广泛关注. 众所周知, 许多密码体系都是基于困难问题建立的. 在格中, 最短向量问题(shortest vector problem, SVP)就是一个很著名的NP-Hard 问题[1], 因此对它的研究至关重要.

1.2 寻找短向量

至今为止, 在寻找短向量方面已经有了很多种方法. 例如枚举[2,3],(启发式)筛法[4–6], 以及Micciancio 提出的构造Voronoi 网格法[7], 还有2014 年提出的基于高斯采样的方法[8]等等. 枚举法从字面意思理解就可以看出它具有很高的时间复杂度(2Θ(nlog(n))), 但是其空间复杂度很低(poly(n)). 筛法通常基于启发式假设, 具有相对较低的时间复杂度(2Θ(n)), 较高的空间复杂度(2Θ(n)), 以及较大的常数因子. 从理论上来讲, 在渐进意义下, 启发式筛法应该在时间上优于枚举法, 但是在目前的实践中, 相比于应用了高效剪枝的枚举, 筛法仍然要慢一些[9].

1.3 筛法

2008 年, Nguyen 与Vidick 提出了第一个启发式筛法(NV-Sieve)[5], 这是在SVP 问题中第一次引入启发式假设的解法. 虽然其效果不佳, 原始版本仅能解决不到50 维的SVP 问题, 但是它为解决SVP问题提供了一个新的思路. 2010 年, Micciancio 提出了一种新的筛法(GaussSieve)[6], 从实践角度来讲,它的表现非常出色, 以至于后续许多人的优化版本都是以该方法为基础的[10–14].

1.4 位置敏感哈希LSH

位置敏感哈希(locality-sensitive hashing,LSH)是在解决近邻搜索问题(nearest neighbor searching,NNS)[9]中提出的. 它的思路是保留敏感信息, 模糊其余信息, 使得具有相同敏感信息的向量获得同样的哈希值. 这里的敏感信息可以理解为距离、角度等等. 在搜索近邻的时候, 只需要搜索与目标向量具有相同哈希值的向量即可. 通过计算哈希值, 降低搜索成本, 这个思路在筛法中也是可以应用的.

1.5 筛法与LSH

2015 年, Laarhoven 提出了对于筛法的一种优化[15], 其中第一次引入了LSH. 在筛法中引入LSH,效果是显著的, 但是最初的Spherical LSH 似乎并不适用于GaussSieve, 只能用于NV-Sieve, 所以其实用性并不高. 后续他又使用Angular LSH[16]、Cross-Polytope LSH[17]来加速GaussSieve, 得到了更好的效果.

1.6 我们的结果

在本文中, 我们将K-Means LSH 引入GaussSieve, 对其进行了优化, 提高了其效率. 额外引入一个可调整的参数, 我们提高了其使用的灵活性.

2 预备知识

2.1 格

这里先介绍一下格的基础知识. 所谓格, 一句话描述, 就是Rn下离散的加法子群. 但是通常来说, 我们研究的是Zn下的加法子群, 也就是由整点构成的格. 我们定义L = L(B) 是由格基B ={b1,b2,··· ,bn} 生成的格, 即L(B) = {Bx|x ∈Zn}. 定义格的体积vol(L(B)) = det(B), 最短非零向量λ1(B)=λ1(L(B)).

2.2 近邻搜索问题

近邻搜索问题NNS 定义为: 给定向量集合L, 给定一个目标向量v /∈L, 希望找到向量v0∈L 满足d(v,v0) 最小. 如果在低维度的情况下, 向量集合L 的大小为N, 那么通过数据结构的优化与预处理, 可以在O(log(N)) 时间内进行查询. 但是放到高维空间, 随着维数n 增大, 似乎暴力搜索变成了复杂度最低的方法, 这个就是著名的“维数灾难(curse of dimensionality)” 现象[9].

为了解决这个问题, 我们希望将L 转化为哈希桶的结构, 使得对于特定的一个桶中的所有向量vi满足d(vi,v) <= r1, 而对于其余桶中的所有向量vi满足d(vi,v) >= (1+ϵ)r1, 这样我们就可以只遍历一个桶中的向量从而得到结果. 这个方法看上去很完美, 但是它只是理想情况下的方法. 实际上, 人们通常是使用LSH 得到近似结果.

2.3 LSH

一个LSH 是一个函数族, 类似于传统的hash 函数, 它也是将定义域中的向量映射为一个散列值. 不同于传统的hash 函数的地方在于, LSH 会把相似的向量映射为同一个散列值, 而把不相似的向量映射为不同的散列值. 不妨设函数D 是衡量向量相似程度的函数, 那么我们有下面的定义:

定义1[9]一个函数族H = h:Rn→U 定义为在相似函数D 下具有(r1,r2,p1,p2)-敏感的LSH,如果它满足对于任意的v,w ∈Rn, 下面两式成立:

(1) 如果D(v,w)≤r1, 那么Ph∈H[h(v)=h(w)]≥p1.

(2) 如果D(v,w)≥r2, 那么Ph∈H[h(v)=h(w)]≤p2.

根据定义我们也可以看出来, 为了明显地区别v 和w, 我们希望p1≫p2. 然而如果r1和r2比较接近, 那么我们很难满足上面的需求. 因此可以考虑使用多个函数进行比较, 这也是LSH 是一族函数而非一个函数的用意. 为了放大p1与p2的差距, 下面两个操作是至关重要的[9].

AND. k 个不同的哈希函数h1,h2,··· ,hk的AND 操作定义为一个函数h′:Rn→Uk, 其中h′(v)的第j 维坐标即为hj(v), 我们定义h′(v)=h′(u) 当且仅当对于任意i ∈{1,2,··· ,k} 有hi(v)=hi(u).这样任取k 个哈希函数就构成了一个新的LSH, 它是(d1,d2,pk1,pk2) 敏感的.

OR. t 个不同的哈希函数h1,h2,··· ,ht的OR 操作定义为一个函数h′: Rn→Ut, 其中h′(v) 的第j 维坐标即为hj(v), 我们定义h′(v) = h′(u) 当且仅当存在某个i ∈{1,2,··· ,t} 有hi(v) = hi(u).这样任取t 个哈希函数就构成了一个新的LSH, 它是(d1,d2,1 −(1 −p1)t,1 −(1 −p2)t) 敏感的.

在LSH 中任意选取kt 个不同的函数,通过AND 和OR 操作,就能得到一族(d1,d2,1−(1−pk1)t,1−(1 −pk2)t) 敏感的LSH, 注意到如果p1>p2, 那么就有1 −(1 −pk1)t>>1 −(1 −pk2)t, 这就得到了我们想要的LSH. 对于k 和t 的选取, 我们有下面的定理:

定理1[9]对于一个(r1,r2,p1,p2) 敏感的LSH 函数族H, 以及一个大小为N 的向量集合L, 令

对于给定的向量v, 下面两条结论很有可能至少有一条会满足:

(a) 能够找到向量w ∈L 满足D(v,w)≤r2;

(b) 不存在向量w ∈L 满足D(v,w)>r1.

其复杂度为:

(1) 对于列表L 的预处理时间: �O(kN(1+ρ));

(2) 预处理之后列表的空间复杂度: �O(N(1+ρ));

(3) 对于给定v 的一次查询: �O(Nρ).

2.4 K-Means LSH

本文中主要依赖的LSH 是K-Means LSH, 其中的K-Means 最早是机器学习中的最常用的聚类算法[18], 之后Loïc Paulevé 等人根据该想法设计了K-Means LSH[19]. 对于固定的一个K 值, 一个KMeans LSH 函数族H 中的每个函数h 可以描述为:

选取K 个不同的向量构成集合C ={ci|i=1,2,··· ,K}

点集C 通常称作中心点集. 那么通过构造不同的C 就可以得到哈希函数族H. 其具体的敏感度与点集C 的选取有关. 在机器学习中, 中心点集一般是通过随机选取数据集中的点, 然后进行迭代的方式得到的, 例如著名的SIFT 算法[20]. 在这里, 我们启发式地假设格点在空间中分布是均匀的, 因此我们只需要使C 中的点的分布尽可能均匀即可. 下面是启发式假设:

假设1 在格中均匀随机取两个向量v 和w, 其角度Θ(v,w) 与在单位球中均匀随机选取v 和w 得到的角度分布相同.

通俗地理解, 启发式假设意思是说Rn中的格点分布可以看做是均匀的. 然而, 这个均匀的意思与二三维空间并不一样. 事实上, 在高维空间中, 随机取两个点v 和w, 那么他们的角度Θ(v,w) 很有可能非常接近π/2, 而并不是通俗理解的[0,π] 中的均匀随机值. 随着维度的增高, 这个趋势会越发明显.

2.5 GaussSieve

2010 年, Micciancio 等人提出了一种新的筛法[6], 其思路类似于Gauss 在约简二维格基时采用的方法, 因此命名为GaussSieve. 该算法伪代码见算法1.

算法1 GaussSieve(B) [6]Input: 经过LLL 约简的格基B Output: 格中的短向量v 1 Function GaussSieve(B,µ):18 Function GaussReduce(v,L,S):2 Initialize L = {0},S = {},K = 0;3 while K ≤c do 19 while ∃vi ∈L :20 ∥vi∥≤∥v∥∧∥v −vi∥≤∥v∥do 4 if S ̸= ∅then 5 vnew = S.pop();6 else 7 vnew = SampleGaussian();21 v = v −vi;22 end 23 while ∃vi ∈L :24 ∥vi∥≥∥v∥∧∥vi −v∥≤∥vi∥do 8 end vnew = GaussReduce(vnew,L,S);10 if vnew = 0 then 9 11 K = K +1;12 else 25 L = L{vi};26 S.push(vi −v);27 end 28 Return v;29 End Function 13 L = L ∪{vnew};14 end 15 end 16 Return L 中最短向量;17 End Function

该算法的主要空间复杂度是栈S 和列表L 的大小, 主要空间复杂度是从L 中查找能约化向量的时间. 因此, 从减小时间复杂度的角度来考虑, 利用LSH 来加快查找速度是一个不错的想法.

对于两个向量u 和v, 如果∥u −v∥≤∥u∥或者∥u −v∥≤∥v∥, 那么显然有Θ(u,v)≤π/3, 反之则有Θ(u,v) > π/3. 看上去似乎不好区分两种情况, 但是事实上, 在高维空间中, 能够约化的向量对大概率满足Θ(u,v) = π/3, 不能约化的向量对大概率满足Θ(u,v) = π/2. 我们曾经做过一个实验, 对于50 维的空间, 随机选取向量对u 和v, 那么Θ(u,v) 有超过99% 的几率落在区间[85°, 95°] 中. 所以我们只需要选取形如(π/3,π/2,p1,p2) 的LSH 即可. 即我们考虑用

来代替

进行判定.

在上面的讨论中, 对于我们选取的参数为(π/3,π/2,p1,p2) 的LSH, 我们实际上只需要关心它的参数r1= π/3. 对于第二个参数r2= π/2, 实际上只是表明该LSH 能够区分夹角小于r1= π/3 与夹角接近r2= π/2 的向量对. 严格来说, 可以改写参数为r2= π/2 −δ 使我们的表述更加严谨, 其中δ << r2为一个常数. 对应的, p2也会增加一个可以忽略的常数, 这是因为高维空间中随机选取的向量对几乎都具有接近0 的内积, 即夹角会大于r2= π/2 −δ. 这些向量对是不能约化的, 因此利用具有r1= π/3 参数的LSH 能够帮助我们筛选出能够约化的向量对, 从而达到加速算法的目的.

3 应用LSH 的GaussSieve

3.1 算法描述

结合上面介绍的内容, 很自然的考虑用Θ(u,v) 来估算是否可以对u 和v 进行约简, 引入LSH 来估算Θ(u,v) 的大小. 于是考虑用一组哈希桶Ti代替GaussSieve 中的列表L, 用同一个哈希桶中的搜索取代原算法中对列表L 进行蛮力搜索的步骤, 得到算法2[16].

算法2 GaussSieve-with-LSH(B) [16]Input: 经过LLL 约简的格基B Output: 格中的短向量v 1 Function GaussSieve-with-LSH(B,µ):21 Function GaussReduce-with-LSH(v,L,S):2 Initialize L = {0},S = {},K = 0;3 Initialize k ∗t random hash functions hi,j ∈H;4 Initialize t empty hash tables Ti;5 while K ≤c do 22 C = ∪t i=1 Ti[hi(v)];23 while ∃vi ∈C :24 ∥vi∥≤∥v∥∧∥v −vi∥≤∥v∥do 6 if S ̸= ∅then 7 vnew = S.pop();8 else 25 v = v −vi;26 end 27 while ∃vi ∈C :28 ∥vi∥≥∥v∥∧∥vi −v∥≤∥vi∥do vnew = SampleGaussian();10 end 11 vnew=GaussReduce−with−LSH(vnew,L,S);12 if vnew = 0 then 9 13 K = K +1;14 else 29 L = L{vi};30 Delete v from all hash tables;31 S.push(vi −v);32 end 33 Return v;34 End Function 15 L = L ∪{vnew};16 Insert v to all hash tables Ti;17 end 18 end 19 Return L 中最短向量;20 End Function

其实相比于以往的此类算法, 他们的框架(算法2) 都是相同的, 区别只在于LSH 的具体形式.

3.2 K-Means LSH

根据2.4 节中介绍的内容明显能看出来, 我们希望中心点集的分布尽可能均匀, 例如K =2n 的时候,取中心点集为C = {±ei} 或其一个旋转. 但是事实上, 对于高维空间, 这是一个名叫Tammes 的困难问题[21], 考虑到我们求解的格至少有50 维, 在这样的维度下求解Tammes 问题并不可行. 那么我们退而求其次, 我们随机选取K 个点组成中心点集C, 对于某个固定的角度α, 下式成立

上式通俗的理解就是希望C 中的点尽量相互远离. 为此, 我们设计了一个很简单的算法可以实现这个目标, 即算法3.

算法3 很简单, 而且作为筛法的预处理步骤, 我们并不需要过多的关注它的复杂度, 这一切似乎都很完美. 然而事实上, 如果假设α = π/2, 明显最优解是|C| = 2n, 但是如果使用算法3, 甚至连n 个点都很难取出来. 万幸的是, 虽然对于任意维度, 在α = π/2 下的最优解都是|C| = 2n, 但是随着α 的减小, 例如α = 80◦,70◦, 那么通过算法3 得到的|C| 与维度n 至少是指数级的关系, 这个在文献[12] 中有详细的描述. 特别的, 我们如果取α=π/3, 在n=10 的情况下, 利用算法3 得到的|C|≥104. 对于α=π/3 的情况, 文献[5] 的结论告诉我们至多存在20.2075n个两两夹角大于α 的向量, 而在实际使用中, 我们只需要取K = kn 个中心点集. 因此, 在本文的主要算法中, 使用这个简单的取样算法是合理的. 得到中心点集后, 就可以根据(1)式计算出hash 值以完成算法2.

算法3 Central point set sampling Input: 最小角度α Output: 中心点集C 1 Const 最大碰撞次数m;2 Initialize C = {}, 当前碰撞次数k = 0;3 while k ≤m do k = k+1;7 else 4 p = Simple_on_Unitsphere();5 if ∃p0 ∈C : Θ(p,p0) ≤α then 6 8 Insert p into C;k = 0;10 end 11 end 12 Return C;9

3.3 与CP-Sieve 对比

先回忆一下CP-Sieve 中使用的Cross-Polytope LSH[22], 其哈希函数为

其正负号等于绝对值最大坐标的符号. 然后将这个哈希函数拓展为一族哈希函数

其中A 是满足每一个元素ai,j∼N[0,1] 的矩阵, N[0,1] 表示[0,1] 上的正态分布. 现在让我们从几何的角度理解一下这一族函数: 显然h 的取值跟v 的模长无关, 所以不妨将其看成是单位球面上的哈希函数,A 可以看做是旋转矩阵. 因此其几何含义便是找到C = {±ei} 中距离目标v 最近的点, 其函数族即为对C 进行了旋转A−1.

因此我们可以看出, Cross-Polytope LSH 是K-Means LSH 在α = π/2 时的最优解, 可以看做是K-Means LSH 的一种特殊情况. 相比于Cross-Polytope LSH, K-Means LSH 可以通过更大的C 以及更小的α 来提高其精确度.

对于Cross-Polytope LSH, 文献[23] 中给出了其单次的精确程度, 如下引理

引理1 记H 为Cross-Polytope LSH 函数族, 记θ =Θ(u,v) 为两个向量u 和v 的夹角, 那么对于较大的n, 有下式成立

虽然这个引理看上去结果很漂亮, 但是事实上, 在我们的应用范围内, 即维数n 相对比较小的时候,O(log log n) 对于性能的影响是相当大的. 例如在n=50 的时候, CP-LSH 的实验结果显示其p1≈3.5%,p2≈1%, 这与引理1 去掉O(log log n) 的结果相去甚远. 在实践中, K-Means LSH 的表现相对优秀一些,针对这方面的对比, 我们会在后面进行展示.

从计算hash 值的复杂性上对比, CP-LSH 的计算过程是旋转O(n2) 与寻找最大坐标O(n), 通过寻找特殊形式的旋转矩阵A, 可以将旋转的复杂度降为O(n log n). K-Means LSH 的计算过程是计算K 个内积O(Kn) 与寻找最大值O(K), 考虑到K 的大小一般正比于n, 因此其复杂度也是O(n2), 这与未经优化的CP-LSH 复杂度相同. K-Means LSH 在复杂度上还有优化空间, 这是一个开放性问题.

另外, Thijs Laarhoven 提出了一种步进式的优化[14], 能够提高框架(算法2) 的性能, 也可以对本算法进一步优化.

4 实验结果

本文中提到的K-Means-LSH 与机器学习中聚类分析使用的同名LSH 并不完全相同, 这是因为聚类分析中是根据已有的固定的样本点通过迭代的方式生成中心点集, 而本文针对的是球面上均匀分布的点确定中心点集, 相比前者, 该方法更像对球面的均分, 这与文献中提到的Angular-LSH 和CP-LSH 的目的是一样的. 为了观察K-Means-LSH 能否适应这个目的, 本节中前两小节先对该LSH 进行一些实验测试.

4.1 不同K 值对应LSH 的性能对比

本节中, 我们观察不同K 值对于LSH 性能的影响, 从而在后面的分析中选择更合适的LSH. 在这里, 我们只考虑单个LSH 的结果. 对于单个LSH 而不是LSH 函数族, CP-LSH 的复杂度是O(n), 因为不需要进行旋转, 从均匀取样的角度来看, 只需要寻找最大的坐标分量即可, 然而K-Means-LSH 的复杂度是O(Kn). 受限于计算能力, 我们只在中维度(n ∈[30,70]) 进行观察. 为了结果的准确性, 我们取了r2=85°. 结果如图1.

由于该LSH 的具体性能受中心点集位置的影响, 而均匀布点又是困难问题, 所以这里只能通过多次试验给出数值模拟, 具体的理论值有待研究.

图1 记录了该实验的结果. 看起来似乎违背了我们的直觉, 随着K 的增大, 参数p1有明显的的下降趋势. 事实上, 回到定义, 实验中的p1表示的是下式的值:

那么我们考虑哈希值不同且夹角小于π/3 的向量对, 这样的向量对会出现在两个中心点“控制范围” 的边界附近一定范围内(中心点p 控制范围指的是集合{u : Θ(u,p) ≤minp′̸=pΘ(u,p′)}, 其边界附近的一定区域可以用集合∆p= {u : minp′̸=pΘ(u,p′)−δp≤Θ(u,p) ≤minp′̸=pΘ(u,p′)} 来表示, 其中δp是常数). 随着K 的增大, 这些区域的并集∪p∆p会具有更大的面积, 从而导致p1的下降.

4.2 K-Means-LSH 与CP-LSH 的性能对比

本节中, 我们来对比两种不同的LSH 在不同维度下的p1与p2值, 以此反映两种LSH 的性能. 考虑到实际应用, 我们取算法3 中的α=π/3, 结果如图2(a).

从实践结果上来看, K-Means LSH 的效果要明显优于CP-LSH. 当K = N 时, 值p1与p1/p2都是K-Means-LSH 更优. 当K = 2N 时, 两种LSH 具有相似的p2值, 这很好理解, 因为CP-LSH 就是均匀情况下的K-Means-LSH. 然而, K-Means-LSH 的p1值远大于CP-LSH. 这种情况随着维数的增加表现的更明显, 因为在高维, 随机取两个向量, 它们之间的夹角有很大概率接近π/2, 在这种情况下, 由于中心点集的不均匀会造成某些局部中心点比较密集, 这个子中心点集能够过滤接近但是小于π/2 的点对, 从而提高p1值. 随着K 的增大, K-Means-LSH 的p1会降低, p2也会略微降低, 但是p1/p2会升高(如图2(b)).

4.3 K-Means Sieve 性能测试

本节中, 我们将对使用了K-Means-LSH 的GaussSieve 进行测试. 本算法涉及三个参数: AND 操作次数k、OR 操作次数t、中心点集大小K. 根据定理1 与文献[17], 考虑到两种LSH 的参数p1与p1/p2对比, 我们取t = 30, 对不同的K, 在维度[35,70] 之间测试其渐进复杂度. 在这里, 为了排除计算机性能对于算法的干扰, 我们先同文献[17] 一样, 考虑Operations 参数(这里Operations = Comparisons +Hashes, 第一项表示比较两个向量能否约减的次数, 第二项表示计算hash 值的次数). 测试结果如图3(图中的小数表示直线的斜率):

从图3 可以看出, 我们的算法在具体复杂度的数值上略低于CP-Sieve, 这是因为CP-LSH 在计算量上相当于K = N 的K-Means LSH, 但是由于中心点集分布不均匀, 在K = N 的时候效率比CPSieve 差. 但是从渐进复杂度的角度来看, 我们的算法是有优势的, 随着K 的增加, 渐进复杂度是优于CP-Sieve 并且还有下降空间的. 更重要的是, 上面采用Operations 作为比较参数其实并不是很合适, 因为计算一次Comparison 需要进行O(n) 次运算(加法与乘法), 而计算一次hash 需要进行O(n2) 次运算(CP-LSH 需要计算一次旋转, K-Means LSH 需要计算K = kn 次内积), 因此这里我们考虑使用参数Complexity=Comparisons+Hashes2来衡量其复杂度, 测试结果如图4(图中的小数表示直线的斜率):

从图4 中可以看出, 在Complexity 参数下, 无论是在渐进复杂度意义下还是具体数值下, 我们的算法都是有优势的.

综合来看, 我们的算法相比于现有的CP-Sieve 是有优势的, 不仅仅是因为复杂度上的优势, 更因为我们的算法多一个参数K, 使得在实际应用中有更大的选择空间, 对于较低维度的空间, 可以不考虑渐进复杂度以获得更优的时间性能, 而对于高维空间, 可以通过增大K 以获得更好的Operation 意义下的渐进复杂度(这点从图中能看出, 对于固定的k, 增大参数K 的值能明显优化其渐进复杂度).

由于Tammes 问题的复杂性, 对于高维空间, 我们甚至无法对均匀布点的个数进行估计, 从算法3 测试的结果来看, 它至少是指数级别的. 从实验结果上来看, 对于略高维度的空间, 增大K 的效果是明显的,这就说明中心点集与维度n 之间的关系在最优情况下很可能不是线性的. 因此, 中心点集的大小或许可以进一步优化. 除此之外, 它还有两个很大的优化空间: 一是能否通过选取特定形式而又不失均匀性的中心点集, 以达到降低计算hash 值复杂度的目的; 二是优化中心点集的均匀性, 以提高其性能.

总的来说, CP-Sieve 可以看做是本方法的特殊情况. 受限于Tammes 问题, K-Means Sieve 的潜力并没有被完全地挖掘出来. 考虑到K-Means LSH 性能优秀, 我们完全有理由相信, 当中心点集足够均匀的情况下, K-Means Sieve 会有更好的结果.