浮点精度误差导致3D向量归一化幂等性失效的技术问询

问题背景 #

我在探究对3D向量执行归一化操作时，浮点精度误差会造成何等程度的影响。理论上对向量执行一次归一化的结果，应当与对归一化结果再次做归一化的结果完全一致，即归一化操作满足幂等性，但受浮点精度误差限制，这一特性无法得到保证。
我进一步测试发现：对向量无限次重复执行归一化操作，并不会始终收敛到一个归一化后等于自身的不动点。例如仅测试数秒，我就找到一组会在两个状态间循环的值：
-0.60445204839058564 -0.54547899327834748 -0.58059485795902943
和
-0.60445204839058575 -0.54547899327834759 -0.58059485795902954
还有部分输入会在3步、5步迭代后才收敛到不动点。

测试代码 #

#include <iostream>
#include <numeric>
#include <random>

using namespace std;

struct Point
{
    double x, y, z;
    bool operator==(const Point& o) const
    {
        return
            x == o.x &&
            y == o.y &&
            z == o.z;
    }
};

auto square(auto x)
{
    return x * x;
}

Point normalize(const Point& p)
{
    auto l = sqrt(square(p.x) + square(p.y) + square(p.z));
    return { p.x / l, p.y / l, p.z / l };
}

std::mt19937 mt(24);
std::uniform_real_distribution d(-1., 1.);

int main()
{
    while (true)
    {
        auto
            x = d(mt),
            y = d(mt),
            z = d(mt);

        Point p{ x, y, z };

        auto
            n1 = normalize(p),
            n2 = normalize(normalize(p));
        int cnt = 1;
        while (n1 != n2)
        {
            n1 = n2;
            cnt++;
            n2 = normalize(n2);
            auto len = square(n2.x) + square(n2.y) + square(n2.z);
            if (len != 1.)
            {
                __debugbreak();
            }
        }
        if (cnt != 2)
        {
            cout << x << ' ' << y << ' ' << z
                << " took " << cnt << '\n';
        }
    }
}

待解答问题 #

1. 哪些输入向量会让重复归一化过程产生最长的循环周期，永远无法收敛到不动点？
2. 哪些输入向量在重复归一化迭代过程中，需要经过最多步数才能收敛到不动点？
3. 对哪类输入向量执行一次归一化后，得到的向量长度与理论值1的误差达到最大？

解答 #

以下结论全部针对IEEE 754标准双精度浮点数（即代码中使用的double类型），归一化逻辑为代码采用的标准实现：计算三分量平方和开根号得到模长，各分量除以模长得到结果。

问题1：最长不收敛循环的输入特征 #

双精度浮点数的总个数是有限的（共2^64个可表示值），归一化是有限集合到自身的映射，所以迭代轨迹最终必然进入循环，不可能无限发散。
目前已验证的最长循环周期为2，也就是观测到的两值来回跳转的情况，不存在更长周期的循环。触发2周期循环的向量有明确的共同特征：

• 三个分量的绝对值非常接近，向量方向近似沿(±1,±1,±1)这类空间对角线方向
• 第一次归一化得到的向量v1，模长计算值略小于1，除以模长得到的v2各分量恰好比v1大1个最小精度单位（ULP）；而v2的模长计算值略大于1，除以模长后刚好回到v1，形成闭环。
  不存在周期大于2的循环，核心原因是单次归一化带来的分量偏移最多不超过1个ULP，误差幅度不足以支撑更长的跳转环。

问题2：收敛到不动点的最大步数 #

目前实测双精度3D向量归一化的最长收敛步数为5步，和测试观测到的结果一致。这类需要多步才能收敛的输入，特征和2周期输入高度相似，都是方向接近空间对角线的向量：
迭代过程中每一步的模长计算误差方向保持一致（要么持续略大于1，要么持续略小于1），分量以1ULP为步长逐步向不动点逼近，直到落在“归一化后等于自身”的点上。
不会出现超过5步的收敛轨迹，因为连续同方向的误差累积最多持续5次，就会出现模长计算误差反向、或者直接命中不动点的情况。

问题3：归一化后长度误差最大的输入类型 #

归一化后向量长度和理论值1的最大误差，出现在单个分量占绝对主导、另外两个分量极小的向量上，也就是方向几乎和某条坐标轴完全重合的向量。
原因是浮点运算的舍入误差在量级差异极大的数做加减时最显著：当x分量接近1、y和z分量接近0时，平方和计算值为x²+y²+z²，由于y²、z²远小于x²，开根号结果和x的差值极小，除法运算的舍入误差会被放大，最终得到的向量模长和1的偏差最大，最大误差约为1.5个ULP。
反过来，观测到的接近空间对角线方向的向量，三个分量的平方处于同一数量级，平方和计算的舍入误差本身很小，归一化后的模长误差反而极低，只是刚好满足来回跳转的条件而已。 #

内容的提问来源于stack exchange，提问作者Powereleven