贪心算法求快速平方根倒数算法中的“魔术数字”【含matlab源代码】
贪心算法求快速平方根倒数算法中的“魔术数字”【含matlab源代码】
一、快速平方根倒数算法简介及实现
1.1 算法简介
在计算平方根的倒数时，传统的计算方法是先计算a的平方根sqrt(a)，再计算它的倒数1/sqrt(a)。但在计算平方根时使用了牛顿迭代法，大量的浮点运算速度很慢。
而快速平方根倒数算法则将输入的32位浮点数a作为一个整数i，用“魔术数字”0x5f3759df减去i右移一位的值产生近似的估值y，再使用牛顿迭代法迭代一次，就得到了相当精度的计算结果。
此算法仅通过简单的位操作和整数减法就得到了一个近似的估计值，浮点运算只包含一次牛顿迭代，可以极大的提高运算速度。
其c/c++源代码如下所示：
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // what the fuck?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
return y;
}1.2 Matlab实现
此算法巧妙的利用了32位浮点数的编码格式。但通过指针将32位浮点数转化为32位整数的运算（以及其逆运算）很难在matlab中实现，但很容易通过c/c++实现。因此我们使用c++实现了float2int32和int32_2float这两个函数，它们将输入的浮点（整数）向量/矩阵中每一个元素转化为整数（浮点数）。编译为mex文件后它们可以被matlab直接调用。
程序列表
二、使用贪心算法计算“魔术数字”
2.1 计算“魔术数字”的最优化问题
理想的“魔术数字”（记为R）应当能够让计算结果的偏差最小化，也就是如下的非线性整数规划问题：
式中目标函数Cost(R)代表R作为魔术数字时算法的误差。为了较好的衡量计算误差，我们在正实数范围内选取1000个点，记为 as；通过传统的取平方根再求倒数方法计算出的精确值，记为rs。代码如下：
as=single(10.^linspace(-8,8,1000));
rs=1./sqrt(as);
再调用FastInvSqrtFloat函数计算出R作为魔术数字时的近似结果es，用es与rs的平均均方误差作为的值。
2.2 初始值的选取
由于自变量只有一个R，我们使用相对简单的贪心算法求解R的值。贪心算法的初始值最好位于全局最优解的邻域，否则往往会陷入局部最优解。快速平方根倒数速算法的思路表明通过R计算出的估值ei = R - bitshift(ai,-1)是比较好的近似结果，反而言之，用精确结果对应的整数ri代替ei应当可以计算出比较接近的R。因此我们分别将前面提到的as和ri它们转化为32位整数向量ai和ri，再计算出R的初始值。计算流程如下如下：
ai=float2int32(as);
ri=float2int32(rs);
R=int32(mean(ri+bitshift(ai,-1)));
得到R的初始值为1597311331，即0x5F350963，与参考代码中的0x5F3759DF（参考值）较为接近。记迭代初始值。在区间[-10^7+R0,10^7+R0]内均匀选择2001个R值计算对应损失函数值，得到R0邻域内R-Cost(R)的关系：
如上图所示，参考值和初始值非常接近最优解，可以通过简单的贪心算法迭代求得最优解R*。
3、贪心算法计算
我们假设迭代的步长为100，每次迭代时对当前的R值施加范围内的随机扰动：R(n+1)=R(n)+RAND, RAND∈[-100,100]。如果Cost[R(n+1)]＜Cost[R(n)]，则接受新的R作为当前值，否则放弃新的R，重新施加随机扰动。迭代直到连续10000个新的R值不被接受为止。
经过迭代，一共采纳过1428个新的值，误差由3.0944逐步下降至2.4536。最终得到最优解R* = 1597385922 = 0x5F362CC2。迭代结果如图所示（以更高的精度绘制了最优解邻域内误差函数的图像）。
三、检验
上面的太过抽象？下图为你展示了使用快速平方根倒数法计算4.0的-0.5次方的过程。
随机选取几个数，以上文所得的R* = 1597385922为“魔术数字”，用快速平方根倒数算法计算，结果如下：
可见使用该“魔术数字”的计算精度较高。
四、讨论与总结
4.1 计算出的最优解0x5F362CC2与参考源代码中的0x5F3759DF不同，因为两个数字的来源不同。前者源于贪心算法寻找得到的最优解，后者则来源未知，也许是理论计算所得。
4.2 最优的“魔术数字”的值与待计算的a值的分布有关。我认为对于特定用途（如光照渲染）的快速平方根倒数算法可以统计a值的概率分布（如需要正规化的向量二范数的分布），根据特定的a值分布来改进Cost函数，再通过最优化方法计算出特定用途下误差最小的“魔术数字”。
4.3 双精度浮点数同样可以采用该算法，只需将代码中的单精度浮点数换为双精度浮点数，32位整数换为64位整数即可。
本文涉及到的完整程序已上传至matlab编程爱好者Q群，如有需要的伙伴请在公众号中回复“QQ”加群领取，在群文件matlab爱好者公众号数据及程序文件夹下的《快速平方根倒数算法》。。
感谢TokiNobug给公众号投稿，欢迎更多爱好、喜欢matlab编程的朋友来稿，在公众号回复“投稿”了解投稿详情。
参考资料：
[1] https://www.zhihu.com/question/26287650
[2] 平方根倒数速算法_百度百科