Incorrect symmetric eigendecomposition for a 4x4 f64 matrix (swapped eigenvectors)

[acshi]

For a certain 4x4 symmetric matrix, nalgebra does not compute a correct eigendecomposition. The maximum relative error in the recomposed matrix is 3.4 (so 340% off). As far as I can tell, this is because two of the eigenvectors (or eigenvalues, equivalently) are in the wrong order. If we swap them back into the right order, we get a maximum relative error of 4.1e-11. If we instead use peroxide to compute the eigendecomposition, we get a relative error of 1.9e-15.

I've tried looking at the code for the eigendecomposition, but I'm really not familiar enough with the specific algorithms (some kind of implicit QR iteration with various optimizations?) to know how to debug this.

#[test]
#[rustfmt::skip]
fn symmetric_eigen_issue_4x4() {
    let m = Matrix4::new(
        -19884.07f64, -10.07188, 11.277279, -188560.63,
        -10.07188, 12.518197, 1.3770627, -102.97504,
        11.277279, 1.3770627, 14.587362, 113.26099,
        -188560.63, -102.97504, 113.26099, -1788112.3,
    );

    let swap_wrong_columns = false;
    let use_peroxide = false;

    let (evals, evecs) = if !use_peroxide {
        let mut eig = m.clone().symmetric_eigen();
        if swap_wrong_columns {
            let (eval1, eval2) = (eig.eigenvalues[1], eig.eigenvalues[2]);
            eig.eigenvalues[1] = eval2;
            eig.eigenvalues[2] = eval1;
        }

        (eig.eigenvalues, eig.eigenvectors)
    } else {
        let eig = eigen(&peroxide::structure::matrix::Matrix {
            data: m.iter().copied().collect(),
            row: 4,
            col: 4,
            shape: peroxide::structure::matrix::Shape::Col,
        });
        let evals = nalgebra::Vector4::from_iterator(eig.eigenvalue.iter().copied());
        let evecs =
            nalgebra::Matrix4::from_iterator(eig.eigenvector.data.iter().copied());
        (evals, evecs)
    };

    eprintln!("evecs: {evecs:.4e}");
    eprintln!("evals: {:.4e}", nalgebra::Matrix::from_diagonal(&evals));

    let recomp = evecs * nalgebra::Matrix::from_diagonal(&evals) * evecs.transpose();
    let diff = (m.lower_triangle() - recomp.lower_triangle()).abs().component_div(&m.lower_triangle().abs().add_scalar(1e-12));
    eprintln!("relative error: {diff}");
    eprintln!("maximum relative error: {}", diff.max());

    assert_relative_eq!(
        m.lower_triangle(),
        recomp.lower_triangle(),
        epsilon = 1e-5,
        max_relative = 1e-6
    );
}

[Andlon]

Thanks for reporting this @acshi. Can you confirm your version of nalgebra, please?

I also rely on eigen decompositions, so I also have a personal interest in making sure they are reliable. However, I unfortunately don't have the bandwidth right now to look at this, but hope to be able to look into it in the future. Perhaps @sebcrozet already has some thoughts on what could be going wrong.

In the meantime, any help from the community for looking into this would be really great!

[acshi]

I found this bug while using 0.30.1 and then confirmed this behavior on the latest dev branch of the repository.

[acshi]

A perhaps more convenient work-around (than converting between types to use peroxide) is to use SymmetricEigen from nalgebra-lapack.

[JulianKnodt]

Also just ran into this (nalgebra = "0.32.1") when adjusting eigenvalues for many 4x4 symmetric matrices, does the lapack version work for all matrices?