2.2 通用软件优化方法
在这一节我们介绍下面几类通用优化方法，分别是：
·循环结构优化——针对具有多层嵌套的循环执行结构的优化。
·时间空间平衡——增加内存空间换取运算时间降低或者反过来的优化方法。
·运算精度和性能平衡——通过降低运算精度来提升运算速度并降低存储空间的方法。
·底层运算的快速实现算法——对于底层的乘除运算用等效的位操作提速。
·内存使用优化——通过调整运算次序或者复用内存，降低程序对嵌入式系统的最大内存需求。
这几类优化方法除了应用于嵌入式机器学习领域外，也能用于其他嵌入式应用，属于相对通用的优化方法。
2.2.1 循环结构优化
循环结构优化可以通过对比代码清单2-8中给出的三段功能相同的代码来说明，其中第一段和第二段代码对应的循环结构最内层的运算次数相同，都是109次，但对于第一种写法calc1，最内层的循环变量k需要重复加载108次，而第二种写法calc2的循环变量k只需要重复加载4次，降低了循环变量k上的重复操作次数，提高了运算效率。一般而言，如果嵌套循环的循环变量次序可以更换的话，我们会倾向于将循环次数少的循环变量放置到外层。
在calc3中，最内层的循环k被展开，程序行数增加了，但完全消除了额外的循环变量k的运算，因此运行效率会有所提升。
代码清单2-8 不同循环结构的代码示例
void calc1()
{
for (unsigned i=0; i<10000; i++)
for (unsigned j=0; j<10000; j++)
for (unsigned k=0; k<4; k++)
data_proc(i,j,k);
}
void calc2()
{
for (unsigned k=0; k<4; k++)
for (unsigned i=0; i<10000; i++)
for (unsigned j=0; j<10000; j++)
data_proc(i,j,k);
}
void calc3()
{
for (unsigned i=0; i<10000; i++)
for (unsigned j=0; j<10000; j++)
{
data_proc(i,j,0);
data_proc(i,j,1);
data_proc(i,j,2);
data_proc(i,j,3);
}
}
下面再看一下代码清单2-9所列出的矩阵乘法代码中循环结构的优化，考虑不同循环结构的两种矩阵乘法实现代码。
代码清单2-9 不同循环结构的矩阵乘法实现代码
void mat_mul1(int **a, int **b, int **c)
{
for (int i=0; i<100; i++)
for (int j=0; j<100; j++)
for (int k=0; k<100; k++)
c[i][j] = c[i][j] + a[i][k] * b[k][j];
}
void mat_mul2(int **a, int **b, int **c)
{
for(int i=0; i<100; i++)
for(int k=0; k<100; k++)
for(int j=0; j<100; j++)
c[i][j] = c[i][j] + a[i][k] * b[k][j];
}
如果观察最内层循环对应的数组访问，可以看到函数mat_mul1有一个连续内存访问，而mat_mul2为两个。函数mat_mul1的最内层循环变量k使得对数据a[i][k]的访问是连续的，但对数据b[k][j]的访问是跳跃的，如图2-17所示。
图2-17 函数mat_mul1的最内层循环对数据b[k][j]的访问模式
函数mat_mul2改变了循环结构，使得最内层循环变量j涉及的两个数据访问都是连续的，如图2-18所示。
图2-18 函数mat_mul2的最内层循环对数据b[k][j]和c[i][j]的访问模式
即循环变量j使得对b[k][j]和c[i][j]的访问都是连续的。在高性能嵌入式系统中，通常对于跳跃的内存的访问效率比连续内存的低，这是因为CPU内核一般通过Cache访问数据。而Cache通过内存控制器间歇地从外部存储器读取地址连续的一批数据放入Cache。对于连续地址访问，CPU通过一次Cache填充就能够得到所需要访问的多个连续地址数据，相比之下，如果跳跃地访问，CPU需要频繁地启动数据加载硬件填充Cache，运行效率低。
2.2.2 时间空间平衡
在嵌入式系统中，由于成本和功耗约束，对软件的运算速度和内存也有约束，实际编程过程中可以通过时间换取空间或者空间换取时间的方式进行优化。考虑下面的代码，它计算输入数据二进制形式下的1的个数，可以使用代码清单2-10列出的几种方案实现。
代码清单2-10 统计二进制数据x中位1的个数的几种代码实现
int count1(unsigned long x)
{
int c=0;
for (unsigned i=0; i<32; i++)
if (x&(1UL<<i))
c++;
return c;
}
int count2(unsigned long x)
{
int c=0;
int lut[]={0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};
for (unsigned i=0; i<8; i++)
{
c+=lut[x&0xf];
x>>=4;
}
return c;
}
int count3(unsigned long x)
{
x = (x&0x55555555) + ((x>>1)&0x55555555);
x = (x&0x33333333) + ((x>>2)&0x33333333);
x = (x&0x0f0f0f0f) + ((x>>4)&0x0f0f0f0f);
x = (x&0x00ff00ff) + ((x>>8)&0x00ff00ff);
x = (x&0x0000ffff) + ((x>>16)&0x0000ffff);
return (int)x;
}
其中算法count2通过32次循环，对1的个数累加得到答案，而count2通过8次循环得到结果，其中每一次分析x中的4个位，而这4个位中1的数量通过查表lut直接获得。函数count2的代码尺寸和存储空间增加了，但换取了运算速度的提升。
代码片段count3还给出了基于“并行位相加”的方案，该方案利用CPU内部的32位无符号整数加法器的特性进行计算，其中第一行
x = (x&0x55555555) + ((x>>1)&0x55555555);
实现x中序号为奇数的位和序号为偶数的位相加，结果依次存在连续2位空间，如图2-19所示。
图2-19 每两个相邻位相加，运算结果存放于相邻的两个位中
接下来的代码
x = (x&0x33333333) + ((x>>2)&0x33333333);
实现对相邻的两组2位数据相加，结果依次存于连续4位空间，如图2-20所示。
图2-20 每相邻两组2位数据相加，运算结果依次存放于连续4位空间中
后续的代码
x = (x&0x0f0f0f0f) + ((x>>4)&0x0f0f0f0f); x = (x&0x00ff00ff) + ((x>>8)&0x00ff00ff); x = (x&0x0000ffff) + ((x>>16)&0x0000ffff);
分别实现：相邻的两组4位数据相加，结果依次存于连续8位的空间；相邻的两组8位数据相加，结果依次存于连续16位的空间；相邻的两组16位数据相加，结果存于连续32位的空间，如图2-21所示。
图2-21 通过邻两组数据相加统计数据中非零位数目的示意图
上述运算结构巧妙地应用了CPU内部高数据位宽的加法器的特性实现并行加法，快速得到位1的计数结果。
2.2.3 运算精度和性能平衡
不少对速度要求很高的应用，允许牺牲运算精度来换取速度的提升。比如代码清单2-11中通过查表计算三角函数cos的近似结果。
代码清单2-11 基于查表得到余弦函数的实现代码
#define _USE_MATH_DEFINES
#include <math.h>
const float LUT[]=
{
1.0000f, 0.9848f, 0.9396f, 0.8660f,
0.7660f, 0.6427f, 0.5000f, 0.3420f,
0.1736f, 0.0000f, -0.1736f,-0.3420f,
-0.5000f,-0.6427f, -0.7660f,-0.8660f,
-0.9396f,-0.9848f, -1.0000f,-0.9848f,
-0.9396f,-0.8660f, -0.7660f,-0.6427f,
-0.5000f,-0.3420f, -0.1736f,-0.0000f,
0.1736f, 0.3420f, 0.5000f, 0.6427f,
0.7660f, 0.8660f, 0.9396f, 0.9848f, 1.0f
};
float lut_cos_deg(unsigned x)
{
x%=360;
unsigned u=x/10,v=x%10;
float a=LUT[u], b=LUT[u+1];
return (float)(a+(b-a)*((float)v)/10.0f);
}
函数lut_cos_deg计算输入数据x的余弦值，其中x是单位为度（°）的非负整数。将通过查表和线性插值结合得到的结果与进行余弦精确计算的结果的差别小于0.004。图2-22给出了精确计算余弦值和上述近似运算的结果的差别，两种运算方式对应的曲线几乎重合。
下面再给出双曲正切（tanh）运算的近似算法。tanh运算在神经网络的激活函数中会被用到，它的定义为
它和另一种常用的激活函数sigmoid关系紧密，sigmoid函数定义为
图2-22 精确计算的余弦值和使用查表近似算法计算的余弦值之间的差别比较。图a为两个运算结果的比较；图b为计算结果的差的曲线
tanh和sigmoid函数的关系为
传统的tanh运算需要进行至少一次指数运算，运算效率低，但我们可以用代码清单2-12给出的近似算法实现。
代码清单2-12 基于近似公式计算tanh的C程序代码
#include <math.h>
float tanh(float x)
{
if (x>(float) 3.4f) return 1;
if (x<(float)-3.4f) return -1;
x*= 1.0f/3.4f;
x*= fabsf(x)-2.0f;
return x*(fabsf(x)-2.0f);
}
这一运算需要进行3次乘法，运算量低于精确计算tanh，但付出的代价是运算结果有最大0.018的误差。图2-23是精确计算的tanh和上面的近似算法得到的tanh的曲线对比，可见两者相差无几。
最后我们介绍以2为底的浮点数对数运算近似算法[1]。由于浮点数本身是以指数加上尾数的形式存储的，对它取以2为底的对数可以直接用到浮点数表示形式的特点。下面是IEEE754格式的单精度浮点数的格式说明。
图2-23 精确计算的tanh和使用近似算法计算的tanh之间的差别比较。图a为两个运算结果的比较；图b为计算结果的差的曲线
图2-24中的指数位E可以看作无符号整数，上述格式的单精度浮点数代表的数值为
图2-24 IEEE754单精度浮点数格式示意图
对于正浮点数，s=0，代表的数值简化为
对该正浮点数求对数，并应用近似得到
根据上述近似可以给出求单精度正浮点数的近似算法，如代码清单2-13所示。
代码清单2-13 单精度浮点数的近似对数计算
float log2_approx(float x)
{
unsigned long v = *(unsigned long*)&x;
return (float)((v >> 23) & 0xFF)-127.0 + \
(float)(v & 0x7FFFFF) / (float)0x800000;
}
我们在应用中除了计算2为底的对数外，也会计算其他数为底的对数，比如计算
，这可以通过
实现，即使用给出的近似计算算法得到log2x后再乘以常数
即可。
2.2.4 底层运算的快速实现算法
对于很多整数运算，可以通过其二进制的位结构的特点来优化计算。当嵌入式处理器没有高效乘法和除法硬件时，我们可以用加减法以及移位实现很多运算。比如考虑计算下面的乘法运算：
其中x和y都是整数（C语言的long类型），这一运算可以通过代码清单2-14给出的程序片段实现。
代码清单2-14 整数乘法运算的快速实现程序
long b=(x<<3)-x; // b=x*(0111)2 y=(b<<4)+b;
如上所示代码中，乘数(01110111)2被拆成相同的两部分(0111)2，因此计算出b=x×(0111)2后，通过让程序中的b左移4位再和它自己相加，就得到结果了。而b=x×(0111)2可以根据(0111)2=(1000)2-1简化，即
上述过程将乘法分解为移位相加运算。对于特定常数的运算，通过观察可以得到最优的分解方案，也可以利用软件搜索最优分解方案。这一部分内容会在第4章详细介绍。
对于除法运算，我们同样可以使用上述优化方案，这里的关键是“化除为乘”，即把计算x/y的运算改成
。比如计算x/11，我们选取N=8，可以得到
当x为整数类型时，可以用
计算，进一步用移位取代乘除法得到下面的C代码片段：
((x<<4)+(x<<3)-x)>>8
上述运算的精度和选取的位宽N有关（比如在式（2-9）中N=8），N的取值越大，精度越高，但相应的运算越复杂（将乘法转成移位加减后，运算量可能增加）。
下面考虑浮点数和整数的乘法。由于数据移位运算不能直接应用在浮点数上，因此在浮点数和整数的乘法中，需要将移位操作用数据相加实现，比如计算y=x×(01110111)2，但这次x是浮点数，运算过程需要修改成代码清单2-15所示的形式。
代码清单2-15 浮点数和整数相乘的快速实现程序
float b,c; b=x+x; b+=b; b+=b; // b=x*(1000)2 b-=x; // b=x*(0111)2 c=b+b; c+=c; c+=c; c+=c; // c=b*(10000)2 y=c+b;
上述移位操作转成数据相加操作，对于左移n位需要数据反复和自己相加n次，当n的取值较大时，运算效率低。另一个可选的方案是根据浮点数格式的特点，直接将移位值作用到浮点数的指数域上。对于IEEE754标准定义的单精度浮点数（即C语言的float类型）由32位数据表示，其中第23～30位存放浮点数的指数部分，对这个区域加或者减n等效为乘以或者除以2n。根据这一特性可以得到代码清单2-16所示的计算浮点数乘以整数(01110111)2=119的代码。
代码清单2-16 基于IEEE754数据格式实现浮点数和整数快速相乘的代码
// 计算x*119，其中119的二进制形式为01110111
float mul_119(float x)
{
float y, b = x;
*(unsigned long*)&b += (3 << 23); // b=(x<<3)
b -= x; // b=(x<<3)-x
y = b;
*(unsigned long*)&y += (4 << 23); // y=(b<<4)
return y + b; // 返回(b<<4)+b;
}
在没有浮点乘法硬件的CPU平台上，上述代码能够快速实现浮点数乘法，但使用这一代码时需要注意数据溢出问题，实际应用时需要补充异常检测和处理过程。读者可以结合第4章内容理解上述代码的含义。
我们最后给出用加法实现浮点数乘法的近似算法。这一算法用了上一节介绍的浮点数近似对数计算的思想。首先，两个单精度浮点数相乘可以写成下面的形式：
x和y的浮点数表示格式中，指数和尾数域的内容分别用
和
表示，根据上一节近似对数计算方法，有
于是
令
，并且z也用单精度浮点数表示，即
，于是有
对比式（2-13）和式（2-14）得到
注意
上面的公式隐含假设
，如果该条件不满足，则会产生额外的误差。
基于上面给出的z的指数和尾数，可以根据单精度浮点数的二进制格式重新拼接成为单精度浮点数z，整个运算过程通过代码清单2-17给出的程序实现。
代码清单2-17 单精度浮点数的近似乘法运算
float mul_approx(float x, float y)
{
long z=*(long*)&x + *(long*)&y -(127 << 23);
return *(float*)& z;
}
上述代码中隐含了多种运算技巧：首先，式（2-15）中浮点数x和y的指数和尾数相加可以用两个浮点数的二进制形式直接当成整数相加实现；其次，上述运算考虑了乘法结果的符号位可以直接通过两个被乘数的符号位直接相加并丢弃进位得到。代码中-(127<<23)对应式（2-15）中减去127的运算，考虑单精度浮点数的指数位在第23～30位，因此将其左移23位。需要注意的是，在实际应用上述代码时，还应补充额外代码来检测数据溢出问题并进行处理。
2.2.5 内存使用优化
1.缓冲区分配
嵌入式系统只有有限的内存空间，需要通过合理地复用内存来降低存储压力。比如考虑之前的计算图，我们之前计算了它在计算过程中的最大数据暂存量，在实现过程中，需要通过存储器复用策略达到最低内存占用量。
回顾图2-12，根据每个运算环节的输入和输出数据（张量）尺寸，可以看到总共有下面几种尺寸的张量：100×100×3，50×50×8，20×20×8，20×20×64，128，10×10。但在实际运行时，我们可以只分配3块数据缓存区：
·1号数据缓冲区——存储尺寸是25 600（20×20×64）。
·2号数据缓冲区——存储尺寸是30 000（100×100×3）。
·3号数据缓冲区——存储尺寸是3200（20×20×8）。
图2-25中给出了特定运算时间分配方案下的缓冲区的分配方案。
图2-25 计算图计算期间的数据缓冲区分配示意图
图2-25中灰色底色的矩形块代表数据缓冲区，里面数字分别对应1、2、3号缓冲区。虽然各个运算需要的输出缓冲区有大有小，但只要分配时确保输出数据小于缓冲区尺寸即可保证算法正常运行。上述分配方案需要对应的数据缓冲区的总存储量是58 800个数据。代码清单2-18给出了上述执行流程的伪代码。
代码清单2-18 计算图执行过程中缓冲区复用的伪代码
float buffer_1[25600];
float buffer_2[30000];
float buffer_3[3200];
while (1)
{
得到图像数据(buffer_2); // 输入数据存储于buffer_1
执行卷积A(buffer_2, buffer_1); // 输入buffer_2，输出buffer_1
执行卷积B(buffer_2, buffer_3); // 输入buffer_2，输出buffer_3
执行卷积C(buffer_1, buffer_2); // 输入buffer_1，输出buffer_2
执行数据拼接和卷积D (buffer_3, buffer_1); // 输入buffer_3，
// 输出buffer_1
执行全连接层E(buffer_2, buffer_3); // 输入buffer_2，输出buffer_3
执行卷积层F(buffer_1, buffer_2); // 输入buffer_1，输出buffer_2
输出结果();
}
2.原址运算
前面讨论的缓冲区内存分配是假设了每个运算执行期间，输入和输出数据缓存区不能相同，因此每个运算占用的空间是输入和输出数据缓存区的总和，对于内存受限的嵌入式系统，可以通过将部分运算改成“原址计算”的模式进一步降低内存占用量，比如比较代码清单2-19所列出的两个运算（神经网络的ReLU运算）的实现方案的内存使用情况。
代码清单2-19 不使用原址计算和使用原址计算的ReLU激活函数代码
void relu1(float *in, float *out, int size)
{ for (int n=0; n<size; n++) out[n]=max(in[n],0); }
void relu2(float *inout, int size)
{ for (int n=0; n<size; n++) inout[n]=max(inout[n],0); }
代码中，函数relu1的输入和输出是分离的两个缓存区in和out，而函数relu2的输入和输出使用相同的缓冲区inout。可见relu2在运算期间的内存占用总量降低到原先的50%。
3.内存“折叠”复用
原址运算并不能适用于所有的运算环节，但有一大类运算可以部分实现原址计算。下面以卷积神经网络中常用的深度可分离卷积（depthwise separable convolution）为例进行说明，该卷积运算作用于C个输入的二维特征图（矩阵），对输入特征图逐一进行卷积得到对应的C个输出特征图。比如考虑在输入3个特征图的情况下，图2-26a是输入和输出数据存放的内存分布，图2-26b是原始的内存访问模式，图2-26c是改进的内存访问模式。
图2-26 进行深度可分离卷积计算时内存访问模式示意图。图a为原始内存访问模式下占用的总内存；图b为运算的三个阶段实际占用的内存；图c为通过内存复用，将写数据“折叠”到低地址，降低总的内存占用
在图2-26b的内存访问模式中，灰色对应需要保留的内存空间，虚线对应暂时不需要的内存空间。可见在实际运算期间，低地址内存可以依次释放，而高地址内存被依次占用，每一时刻需要占用的内存量小于输入和输出占用的内存之和。根据这一特性，我们可以按图2-26c所示模式复用内存。图2-26c和图2-26b相比，可见由于我们在计算的第2阶段和第3阶段复用了之前释放的低地址内存，因此需要的总内存量小于图2-26a所示的内存访问方案。
上面给出的内存优化方式可以看作“地址折叠”访问方式的一个例子，即输出数据存储时，若写地址超出分配的内存空间上限，则折返到低位地址。代码清单2-20为地址折叠缓存区复用的代码片段，显示了基于地址折叠方式计算并存放输出数据的具体编程模式。
代码清单2-20 地址折叠缓存区复用的代码示例
// buffer是缓冲区，存放输入和输出数据， buffer[0]是第1个输入数据
// offset是运算结果的存放位置，即buffer[offset]对应第1个输出数据
// buffer_size是缓冲区buffer的总尺寸
void calc(float* buffer, int buffer_size)
{
float result;
// 循环处理每个输入数据
for (int n=0; n<N; n++)
{
// 运行数据处理
result=data_processing(buffer+n);
// 存放处理结果，注意数组下标用%buffer_size运算实现存储地址"折叠"
buffer[(offset+n)%buffer_size]=result;
}
}
上述程序运行前后的内存结构如图2-27所示。
图2-27 内存“折叠”访问模式下程序运行前后的缓冲区占用情况示意图。图a为程序运行前输入数据连续存放于缓存区“底部”；图b为程序运行后，输出数据存放位置
由图2-27可以看到缓冲区内部数据的存放位置。在程序运行期间，输出数据由offset依次向上存放，到达缓存区边界后“折返”到缓存器底部。这一访问模式的缺点是输出数据有可能被“切分”成两段，进行后续运算时要根据这一结构访问数据。
图2-28给出能够使用“内存折叠”的代码的内存访问模式的特点。如图2-28a所示，阴影部分代表程序运行期间使用的内存地址范围。程序执行期间，内存占用区间不断改变，呈现图中所示的“斜条带”状，对于这一类代码，内存按时间分配使用方式可以改成图2-28b所示模式，运算执行到后面阶段时内存访问被“折叠”到之前空闲下来的空间。
图2-28 内存“折叠”访问模式示意图，阴影部分代表程序运行期间使用的内存地址范围。图a为优化前的内存占用随时间变化的规律；图b为通过地址“折叠”降低内存占用