浮点数误差分析及补偿


  在各种渲染引擎中使用浮点数几乎可以肯定一定会产生浮点数误差,而渲染引擎的大量计算量也不允许使用其他特别高精度的浮点数,因而需要一定程度的精度补偿。事实上,浮点数并不适合用于精确计算,而适合进行科学计算。

float:2^23 = 8388608,一共七位,这意味着最多能有7位有效数字,但绝对能保证的为6位,也即float的精度为6~7位有效数字;
double:2^52 = 4503599627370496,一共16位,同理,double的精度为15~16位。

相对误差与绝对误差

  度量浮点数误差有两种形式:绝对误差(absolute)和相对误差(relative)。如果我们使用了浮点数计算得到了结果 \tilde a ,我们可以使用实际值和计算差的绝对值来计算绝对误差:
\delta_a = |\tilde a − a|
  相对误差是指绝对误差和实际结果的比值,当 a\neq 0 时:
\delta_r=\lvert {\tilde a − a\over a } \rvert =\lvert {\delta_a\over a}\rvert
  定义了两种误差后,可以这样子表达:
a = a ± \delta_a = a(1 ± \delta_r)

前向误差分析

  下面尝试计算 a,b,c,d 四个数的和,如果我们这样计算:r = (((a + b) + c) + d),加上误差可以表示为:

\begin{aligned} &(((a\oplus b)\oplus c)\oplus d)\\ \subset &((((a+b)(1±\epsilon_m))+c)(1±\epsilon_m)+d)(1±\epsilon_m)\\ =&(a+b)(1±\epsilon_m)^3+c(1±\epsilon_m)^2+d(1±\epsilon_m) \end{aligned}
  因为\epsilon_m非常小,高阶的 \epsilon_m 可以表示为 \epsilon_m 的加法:
(1±\epsilon_m)^n\leq(1±(n+1)\epsilon_m)
  在实践中,(1\pm n\epsilon_m) 一般就可以表示误差的上界,因为高阶的 \epsilon_m 下降得非常快。
  所以,误差的方程可以改写成如下所示:
\begin{aligned} &(a+b)(1±4\epsilon_m)+c(1±3\epsilon_m)+d(1±2\epsilon_m)\\ =&a+b+c+d+[±4\epsilon_m(a+b)±3\epsilon_m c±2\epsilon_m d]\end{aligned}
  就可以得到绝对误差的上界:
4\epsilon_m\lvert a+b\rvert+3\epsilon_m\lvert c\rvert+2\epsilon_m\vert d\rvert
  因此,很容易得到 a,b,c,d 之和的误差上界。
  得到结果中,a+b 的误差贡献比 d 的误差贡献更大,这也说明了如果将大量浮点数加在一起,从小到大排通常比任意排序的误差更小。
  上述分析假定了编译器会根据表达式的顺序来计算和,因此必须要使编译器遵循运算顺序来最小化误差范围。
  如果改变一下相加顺序得到 float \, r = (a + b) + (c + d) 的误差,同样适用上述过程,可以得到误差上界为
3\epsilon_m\lvert a+b\rvert+3\epsilon_m\lvert c+d\rvert
  如果 |a + b| 相对比较大,此误差会比较小,但如果 |d| 相对较大,误差也可能更大。
  这种分析误差的方法称为前向误差分析(forward error analysis):给定输入,可以用相似的方法计算结果的误差上界。得到的误差上界一般会比实际误差更大,但是在实际运算中,误差项的符号往往是混合出现的,所以当相加时会产生误差抵消。
  所以需要一种新的误差分析方法——后向误差分析,将计算结果暂时视为精确,然后给出相同输入产生结果波动的上界。

后向误差分析

  (1±\epsilon_m)^n的保守误差上界(1±(n+1)\epsilon_m)在误差阶数较高时不太满足。Higham 得出了一种能够更加紧逼 (1±\epsilon_m)乘积误差项。
  如果有 (1±\epsilon_m)^n,那么可以值的上界为 1+ \theta_n,当 n\epsilon_m <1 时有 \lvert \theta_n\rvert \le{n\epsilon_m\over1-n\epsilon_m}
  在上式中分母比 1 小 n 个 \epsilon_m 误差,所以几乎没有增大 n\epsilon_m 就得到了一个保守上界。如果用 \gamma_n 来表示上界:
\gamma_n = {n\epsilon_m\over 1-n\epsilon_m}
  使用 γ 标记就可以得到四个值和的误差上界(这里转换更紧逼的 1\pm n\epsilon_m ),有:
\lvert a+b\rvert\gamma_3+\lvert c\rvert\gamma_2+\lvert d\rvert \gamma_1
  这种方法的一个优点是可以计算商的误差,如果有 (1\pm\epsilon_m)^m \over (1\pm \epsilon_m)^n,可以得到误差 1\pm\gamma_{m+n}
  如果计算两个误差的乘积,可以得到误差
a(1\pm\gamma_i)\otimes b(1\pm\gamma_j)\subset ab(1\pm \gamma_{i+j+1})
  故相对误差为
\lvert{ab\gamma_{i+j+1}\over ab}\rvert=\gamma_{i+j+1}
  这说明乘法的误差积累非常小。
  可是对于加减法来说,误差会被积累得很大。对于a(1\pm\gamma_i)\oplus b(1\pm\gamma_j),有误差为
\lvert a\rvert\gamma_{i+1}+\lvert b\rvert\gamma_{j+1}
  所以对于精度很敏感的部分,对加法一定要进行补偿。

PBRT 对浮点数精度的补偿

  在具体的实现中可以使用 constexpr 关键字来指定它为关键字常量:

  以下是 pbrt 中对光线变换坐标时做精度补偿的例子:

  使用了重载的小括号来计算变换后的原点 origin 的位置和浮点数绝对误差,计算误差的方式很简单。

  上述代码值得注意的是 gamma 函数的参数为 3,表达的意思也就是 3\epsilon_m


参考资料:Physically Based Rendering

About the Author

发表评论

电子邮件地址不会被公开。

Bitnami