在Rust中，什么可能导致UnitQuaternion连续乘法的数值错误？如何避免？

Question

假设定义以下函数

get_correction

，使得给定第一次旋转，它可以计算可以将 3d 向量旋转到“向上”方向的第二次连续旋转：

    pub fn get_correction(
        acc: Vector3<f32>,
        rotation: UnitQuaternion<f32>,
    ) -> Option<UnitQuaternion<f32>> {
        Self::get_correction_verified(&acc.to_owned(), &rotation.to_owned())
    }

    fn get_correction_raw(
        acc: &Vector3<f32>,
        rotation: &UnitQuaternion<f32>,
    ) -> Option<UnitQuaternion<f32>> {
        static UP_FRD: Vector3<f32> = Vector3::new(0.0, 0.0, -9.81);

        let uncorrected = rotation.transform_vector(&UP_FRD);
        let correction_opt = UnitQuaternion::rotation_between(&uncorrected, &acc);
        correction_opt
    }

    fn get_correction_verified(
        acc: &Vector3<f32>,
        rotation: &UnitQuaternion<f32>,
    ) -> Option<UnitQuaternion<f32>> {
        Self::get_correction_verified_internal(acc, rotation).0
    }

    fn get_correction_verified_internal(
        acc: &Vector3<f32>,
        rotation: &UnitQuaternion<f32>,
    ) -> (Option<UnitQuaternion<f32>>, f32) {
        let raw = Self::get_correction_raw(acc, rotation);
        match raw {
            Some(correction) => {
                //round trip verification
                let corrected = correction * rotation;

                let should_be_zero = Self::get_correction_raw(acc, &corrected);

                let zero = should_be_zero.unwrap().angle();
                if zero > 0.01 {
                    println!("residual={}", zero);
                    println!(
                        "compute: {}, {} => {}",
                        acc.transpose(),
                        rotation,
                        correction
                    );

                    // let again = Self::get_correction_verified_internal(acc, rotation);
                    // assert!((raw, zero) == again)
                }

                (raw, zero)
            }
            None => (raw, 0.0),
        }
    }

为了尽量减少数值误差，这个定义有一个详细的往返验证，以确保连续2次旋转后的散度接近0。

当使用随机向量和四元数生成器重复测试此函数时，我发现往返误差随着时间的推移逐渐恶化。它从 0 度开始，然后逐渐增加到几乎 180 度。奇怪的是，如果我用相同的参数重新启动测试程序，往返错误将再次变为0。以下是该函数的输入输出对的典型示例：

4分钟后，往返误差=1.6

      ┌                                     ┐
      │   3.0176868 -0.74084723     9.24847 │
      └                                     ┘
,
     UnitQuaternion angle: 2.262818 − axis: (-0.0015257122, -0.9227901, -0.38530007)

 => UnitQuaternion angle: 1.6769453 − axis: (-0.5872864, -0.79925746, 0.1276015)

使用相同的数据重新启动，往返错误=0

        compute:
      ┌                                     ┐
      │   3.0176868 -0.74084723     9.24847 │
      └                                     ┘
    , 
UnitQuaternion angle: 2.262818 − axis: (-0.0015257121, -0.92279005, -0.38530007)

 => UnitQuaternion angle: 0.87873745 − axis: (-0.67957145, -0.71493566, 0.16446783)

由于我故意创建了所有参数的深层副本，因此我认为这个问题不会是由两个参数的可变状态引起的。是什么导致了这种行为？

更新1。经过一些实验，我能够重现此错误：

#[test]
pub fn __accumulate_unit_quaternion() {
    use rand::Rng;

    let mut rng = rand::thread_rng();

    let mut qq = UnitQuaternion::from_euler_angles(0.0, 0.0, 0.0);

    let unit_vec = Vector3::new(1.0, 0.0, 0.0);

    loop {
        let u1 = rng.gen_range(-1.58f32..1.58);
        let u2 = rng.gen_range(-1.58f32..1.58);
        let u3 = rng.gen_range(-1.58f32..1.58);

        let random = UnitQuaternion::from_euler_angles(u1, u2, u3);

        qq = random * qq;

        let reconstruction_error = {
            let reconstructed = UnitQuaternion::from_axis_angle(&qq.axis().unwrap(), qq.angle());
            let error = (qq * unit_vec - reconstructed * unit_vec).norm();
            error
        };

        println!("error={}", reconstruction_error);

        assert!(reconstruction_error < 0.001);
    }
}

执行时，我们可以观察到UnitQuaternion的重构误差逐渐增大，直到断言失败。我的猜想是，线性代数库天真地认为 2 个 UnitQuaternion 可以直接相乘，得到一个 UnitQuaternion，无需归一化，这在理论上是正确的，但在数值计算中并不正确。

如何重写程序以最小化重构错误？如果上述情况属实，我们将需要对库进行短期旁路和长期修复。

Answer 1

根据库维护人员的说法，我需要定期调用：

qq.renormalize()

//OR

qq.renormalize_fast()

它们不会在每次乘法中自动调用，因为它们会产生很小的性能成本。