求解两个向量之间的距离(或者说相似度)，广泛应用于数据挖掘，图像处理和深度学习中。例如自然语言处理中，根据文章中关键字的出现频率来组成一个特征向量，两个文章的特征向量的距离来表示文章主题的相似度。很多人脸识别的模型，在网络的最后一层都是计算出输入人脸的一个固定维数的特征向量，然后根据两个人脸对应的两个特征向量之间的距离来做人脸验证和识别。在真实的应用场景中，文章底库和人脸底库保存的特征向量数量很大(十万百万级)，计算特征向量之间的距离函数就会被调用很多次，所以优化这个搜索过程和最佳距离的计算函数就显得尤为重要。这里，我以优化在100万底库的场景下，找出底库中与给定512维向量之间余弦距离最近的那个样本的索引值为例，试验一下在通用PC平台上到底能优化多少(这里为了演示方便，实验中数据很多次取平均的过程没有展示，请注意)。
欢迎探讨，本文持续维护。
-
操作系统：Ubuntu 16.04 LTS
-
编译器：g++ (Ubuntu 5.4.0-6ubuntu1~16.04.10) 5.4.0 20160609
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.10) 5.4.0 20160609
-
硬件配置：
-
CPU Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz
-
Memory 8G DDR2
-
GPU N/A
-
最基础的实现方式也是最容易的方式，作为我们的起点，实验代码：
main函数和时间测量代码，main函数首先构造一些模拟数据，然后调用SearchBest做搜索，这里我们主要关心的就是SearchBest函数的耗时
SearchBest函数，这个函数调用Cosine_similarity计算余弦相似度，返回最接近的人脸的索引
Cosine_similarity函数，计算两个一维向量之间的余弦相似度
这里为了数据类型切换方便，我把SearchBest和Cosine_similarity函数写成了模板的形式。 测试结果：
耗用了3855587us
测试结果：
耗时793150us，仅仅加了一个-O3的编译选项，什么代码也没改，就加速了x4.86倍。
这里明显很多乘除开方之类的数学操作，如果再开激进一点的优化呢？
测试结果：
耗用了433704us，差不多又优化了加速了x1.82倍。不过，这里要注意一下，最佳人脸的index改变了，这是由于(1)我们在计算的时候都是用的浮点数数据类型，计算结果很精确，在开了-Ofast -ffast-math优化之后，牺牲了计算精度，来换取了计算速度。(2)我们的模拟数据也比较特殊，底库中有很多的特征向量，之间的余弦距离其实很近很密集，这也放大了精度损失带来的影响。
小结一下，由上面Step2和Step3的结果看来，开一下编译器的优化选项，就可以很快地收获到速度上的快速提升(差不多x10倍)，但是也要注意，对计算精度的影响会不会影响业务，有没有规避方式(例如对Feature做特别的预处理让他们变得稀疏)。 下面我们继续
C/C++里面浮点型有两种float和double，float占用字节小，计算速度快，我们把double改成float试试看，修改数据类型：
测试结果：
耗用了201012us，加速x2.16倍(可能是在我的电脑上，一个double是8字节，一个float是4字节，所以取数据时一个cache line填满刚好可以取2x的float数据)。注意，这里index有改变化，讨论和分析见Step2。
OpenMP是一种支持多平台共享存储器多处理器的C/C++编程规范和API，使用简单，只需要在需要做并行的地方加上编译制导语句，然后GCC的编译命令中加上-fopenmp选项即可。
现在的CPU都是多核心的了，我们试试：
测试结果：
耗时123842us，加速了1.62倍，还不错。
因为处理器有多个ALU，如果循环展开能充分利用这些ALU的话，速度应该还有提升。 但是，实验结果显示在循环展开后没有怎么变化，这里就不贴图了。原因的话我觉得可能是在OpenMP开了多线程之后，就已经用满了所有的计算单元吧。
在向量运算中，作用在向量元素上的操作都是一样的，只是数据不一样，这样就很适合用处理器的向量指令集来做加速了。这里我选的是AVX(CPU支持的指令集可以用lscpu或者cat /proc/cpuinfo命令来查看)，AVX指令每次最多可以操作256位的数据，也就是8个32位的浮点数，最理想的情况下，代码可以加速x8倍。代码如下：
辅助函数
计算相似度的函数
另外，因为SIMD指令需要内存对齐，所以在main函数里面申请内存的时候也要申请对齐了的内存。
测试结果：
搜索耗时122652us，和前面优化好像不是很明显。原因可能有(1)可能前面Step2,Step3,Step4已经做了向量化了;(2)AVX优势在提高CPU的计算资源的利用率上，但是也受限于内存带宽，如果内存带宽不够，在数据量大的情况下，很多CPU的时钟周期都浪费在等待数据从内存load进寄存器上浪费了。
一般来说，标准库提供的函数都是由高级的程序员打磨到很极致了的，但是，限于可移植性和跨平台要求，可能还有针对某些平台，某些应用优化的余地。试一下：
上图是QUAKE-III源代码里面计算平方根的倒数的快速计算方法，更进一步的讨论请看Chris Lomont论文
测试结果：
搜索耗时123637us，也不是很明显，一个原因是计算sqrt()这个函数不是热点，它的调用次数和for(int i = 0; i < step; i++){...}这个循环里面的_mm256_fmadd_ps函数调用次数来说不算很多，第二个原因可能数学库里面已经做了很好的优化了。
这是个取巧的办法，现在在特征数据库中的数据存的是100万裸的Feature向量，而在计算余弦距离的时候，每个特征向量都要计算它们的模(长度)，我们可以事前将这100万向量和它们的模一起保存(新入库的特征也一样，加上模一起入库保存)，增加了1/512的(如果加上数据库里面别的字段，这个数字还要稍微小一点)存储空间，但是这样理论上就可以减少1000,000x512次乘加和100,000次开方计算量了。很多查表法也是类似的套路，将一些复杂但是可预测的计算结果预计算，预存储，空间换时间。 实验如下：
测试结果：
速度121683us，变化不大，原因参见Step 7第二条。
余弦距离，只跟两个向量之间的夹角有关，和向量的模(长度)无关。如果我们在存特征的时候，就把特征向量做了模归一化，对新的特征向量也做同样的处理，那么在计算相似度的函数Cosine_similarity_avx里面就无需再计算向量的模了，这样在乘加数上算，理论上可以减少2/3的计算量，加速x3倍。
代码修改处如下，
首先模归一化为单位长度1：
然后把Cosine_similarity函数里面计算两向量模的代码注释掉：
测试一下：
耗时123436us，和前面的差别不大，没有达到预期的x3倍加速。感觉还是因为内存带宽的问题，从内存搬运数据到cache的速度，远小于CPU把cache上数据计算完的速度，所以CPU的计算能力已经冗余了。在我们的优化里面，只减少了CPU的计算量，但是内存访问没有减少（因为还有计算mult_add这一步，所以在内存访问vectorA[i]和vectorB[i]一点都没有少），所以，速度也没用什么提升。
特征向量里面有512维的特征，但是这512维特征的重要性不一定都是一样的，如果我们能有选择性地抽取里面信息量最大的特征，减少特征向量的维数，那么我们就减少了计算量，更重要的是，加大了内存访问数据量，增加了CPU的利用率。PCA就是一种在减少数据集的维度的同时，保持对方差贡献最大(信息量保持最多)的特征的特征选择压缩算法。 因为我们在这里不是探讨PCA算法的过程，我们就直接假设512维的特征向量压缩到了256维。 改动如下：
测试结果：
速度63201us，加速大约x2倍。乘加数减少2倍，内存效率提升2倍，加速2倍符合预期。
OpenBLAS是BLAS的跨平台开源项目，可以加速各种矩阵和线性代数的计算，广泛应用在机器学习应用和框架中，Caffe里面用的BLAS库就是OpenBLAS。
这里我们用OpenBLAS的矩阵与向量相乘的接口cblas_sgemv来做相似度计算。输入矩阵是底库中所有特征向量的集合，向量是需要搜索的那一个特征向量。
实验代码如下：
实验过程结果如下：
耗时63865us，改变不是很大。
在很多通用CPU上，特别是计算能力有限的嵌入式设备中，浮点数的运算要比定点数的计算慢。如果我们把float型的已经归一化到[0.0,1.0]了的浮点型特征值(因为余弦距离只与向量夹角有关，与模还有分量的缩放取值范围无关，所以很容易做归一化)，转换为[0,65536]的unsigned short定点类型，看看能不能有提升。
先做一点理论分析：
（1）在我的电脑上,sizeof(float)等于4，sizeof(short)等于2，如果用unsigned short代替float，内存读取效率可以提升2倍。
（2）unsigned short替换float误差分析
问题描述：两个N维浮点向量X,Y，已经归一化到模为1，即
，其中
，归一化后的余弦距离
假设转定点数时把[0.0, 1.0] 映射到 [0, M], 则分量最大损失为1/M，考虑到四舍五入则为0.5/M。
那么量化导入的误差为
，其中
，带入则有
，两遍加上取绝对值
.
N是特征向量的维数，N=256，M是unsigned short可以表示的最大整形值，M=65536，带进去可以得知相似度的量化误差在万分之三以内。
实际改动：
测试结果：
运行时间31298us，加速x2.02倍。主要原因在于乘加数没有增加，但是内存效率提升了2倍。如果换成4字节的unsigned int，实验证明速度没有什么提升。所以这里的定点化对速度提升的主要收益在于内存效率提升，而不在CPU的定点计算能力比浮点计算能力强。可能也是因为我是用的是PC电脑的通用CPU，浮点计算能力丝毫不逊色与定点计算，也许在嵌入式设备上会不一样。
注意，上面的分析基于特征值没有归一化到[0,1]区间，模没有归一化到1，而且量化误差对相似度影响不大的基础上的，另外，也需要注意溢出问题。
以下是几点心得：
- 不要抱怨你的硬件不好，考虑一下你自己有没有优化好你的程序是不是充分优化了，是不是充分利用了所有可用的计算资源(CPU,GPU,Memory)，是不是把石头里面最后一滴油压出来了;
- 优化无定法，一定要根据具体的问题来做来实验。要在优化前确定优化的目标和约束，优化是一个性能（程序运行时间），资源（CPU、GPU、内存占用，甚至手机电池电量），代码可维护性和跨平台与成本（开发测试时间、人力投入）多者的平衡，在开始动手前最好能对能牺牲什么，一定不能牺牲什么，优化的目标，有个基本的想法，不打无准备的仗;
- 算法优化要在代码优化前面，算法才是程序的核心灵魂，一辆汽车发动机不行，外观再流线型设计也跑不快，牛顿下降比梯度下降快，二分查找比线性搜索快，做代码优化前，先考虑一下算法是不是不能优化了？
- 优化要先找到热点，根据28法则，80%的时间都是被20%的代码消耗的，所以，我们应该放80%的注意力去重点优化那20%的代码。根据阿姆达尔定律，最终你程序的加速比，取决于你对关键热点的优化程度。特别是大的程序，找到程序中最耗时(或耗资源)的地方，下大力气挖潜力;放过非热点部分的优化;要不然就本末倒置，重拳击在棉花上，累死也无功;热点是变动的，一个热点优化了，可能原先排名第二的热点就成了现在排第一的了;
- 优化过程要先易后难，先吃肉后啃骨头，原因有三：(1)一般水平的程序员写的代码和算法或者学究训练出来的模型，都有很大的优化空间，往往开个编译器的-O3，改改数据类型，就能取得不错的效果。(2)从人的心理学上讲，先吃几块肉，看到点优化的效果，也能给自己一点啃骨头的信心和动力。(3)有的优化是重复的，烂程序各有各的烂法，但是优秀的优化，手段都是类似的，开-O3和手动写AVX代码，仔细调整流水线调整memory cache，效果可能是一样的，但是-O3成本多低。先开编译器优化选项，后做缓存优化和汇编，肯定没错的;
- 啃骨头的时候，要善于借助工具来指导优化，各个平台都有各自profiling的工具，也有各自的小实用程序(例如htop, gprof)或者自己开发的收集的小工具，找自己顺手的用，要依赖于自己的经验，分析，甚至直觉，但是要更相信工具和数字，数字不会撒谎;
- 在优化的过程中要试，要记下各个优化点和优化的时间大小，还要仔细检查优化没有影响输入和输出，分析输入输出不一样的原因，最好做个表格，万一修改没有取得好的效果，及时回滚;优化好了也要思考背后的道理，不能知其然不知其所以然，及时总结经验;
- 如果优化手段比较trick，做好记录和代码注释;
- 优化不必太早，在开发阶段，更重要的是代码的规范性，结果的正确性和开发进度;高神有言:
程序员浪费了太多的时间去思考和担忧程序中那些非关键部位的速度，而且考虑到调试和维护，这些为优化而进行的修改实际上是有很大负面影响的。我们应当忘记小的性能改善，97%的情况下，过早地优化都是万恶之源。 ——高德纳
- 优化也是有成本的，例如人力投入，资源占用等，而且优化往往会使代码难懂难维护，还可能引入新的bug，使代码产生对平台的依赖，不到万不得已，能不优化就不要优化
本文只是抛砖引玉，初步涉及了优化汪洋大海中的一点点。光优化的层次来分就有系统级的，应用级，代码级别的，优化还要区分是网络瓶颈，IO瓶颈，内存瓶颈还是CPU瓶颈，本文的这个小demo只涉及到代码级别的单个通用CPU和内存优化，还没有涉及到用上GPU(十分适合做大规模并行计算)，负责均衡。
| 优化步骤 | 时间(us) | 加速倍数(与上一步相比) | 总加速比(和基准线比) | 备注 |
|---|---|---|---|---|
| Step 1(基准线) | 3855587 | x1.0 | x1.0 | 基准线 |
| Step 2(-O3) | 793150 | x4.86 | x4.86 | |
| Step 3(-Ofast -ffast-math) | 433704 | x1.82 | x8.90 | |
| Step 4(double->float) | 201012 | x2.16 | x19.18 | |
| Step 5(OpenMP) | 123842 | x1.62 | x31.13 | |
| Step 6(循环展开) | N/A(与第五步相差不大) | ~x1.0 | ~x31.13 | |
| Step 7(SIMD) | 122652 | ~x1.0 | ~x31.13 | |
| Step 8(优化sqrt()) | 123637 | ~x1.0 | ~x31.13 | |
| Step 9(预计算) | 121683 | ~x1.0 | ~x31.13 | |
| Step 10(数据模归一化) | 123436 | ~x1.0 | ~x31.13 | |
| Step 11(特征选择) | 63201 | ~x1.96 | x61.01 | |
| Step 12(OpenBLAS) | 63865 | ~x1.0 | ~x61.01 | |
| Step 13(浮点转定点) | 31298 | ~x2.02 | x123.19 | 内存效率 |