浮点数误差分析及补偿

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

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 关键字来指定它为关键字常量：

inline constexpr Float gamma(int n) {
     return (n * MachineEpsilon) / (1 - n * MachineEpsilon);
}

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

inline Ray Transform::operator()(const Ray &r) const {
    Vector3f oError;
    Point3f o = (*this)(r.o, &oError);
    Vector3f d = (*this)(r.d);
    // Offset ray origin to edge of error bounds and compute _tMax_
    Float lengthSquared = d.LengthSquared();
    Float tMax = r.tMax;
    if (lengthSquared > 0) {
        Float dt = Dot(Abs(d), oError) / lengthSquared;
        o += d * dt;
        tMax -= dt;
    }
    return Ray(o, d, tMax, r.time, r.medium);
}

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

template <typename T>
inline Point3<T> Transform::operator()(const Point3<T> &p,
                                       Vector3<T> *pError) const {
    T x = p.x, y = p.y, z = p.z;
    // Compute transformed coordinates from point _pt_
    T xp = m.m[0][0] * x + m.m[0][1] * y + m.m[0][2] * z + m.m[0][3];
    T yp = m.m[1][0] * x + m.m[1][1] * y + m.m[1][2] * z + m.m[1][3];
    T zp = m.m[2][0] * x + m.m[2][1] * y + m.m[2][2] * z + m.m[2][3];
    T wp = m.m[3][0] * x + m.m[3][1] * y + m.m[3][2] * z + m.m[3][3];

    // Compute absolute error for transformed point
    T xAbsSum = (std::abs(m.m[0][0] * x) + std::abs(m.m[0][1] * y) +
                 std::abs(m.m[0][2] * z) + std::abs(m.m[0][3]));
    T yAbsSum = (std::abs(m.m[1][0] * x) + std::abs(m.m[1][1] * y) +
                 std::abs(m.m[1][2] * z) + std::abs(m.m[1][3]));
    T zAbsSum = (std::abs(m.m[2][0] * x) + std::abs(m.m[2][1] * y) +
                 std::abs(m.m[2][2] * z) + std::abs(m.m[2][3]));
    *pError = gamma(3) * Vector3<T>(xAbsSum, yAbsSum, zAbsSum);
    CHECK_NE(wp, 0);
    if (wp == 1)
        return Point3<T>(xp, yp, zp);
    else
        return Point3<T>(xp, yp, zp) / wp;
}

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

参考资料：Physically Based Rendering

About the Author

Oscar