sarah-ek / faer-rs

Linear algebra foundation for the Rust programming language
https://faer-rs.github.io
MIT License
1.75k stars 56 forks source link

`selfadjoint_eigendecomposition` sometimes produces incorrect results #124

Closed jbncode closed 3 months ago

jbncode commented 3 months ago

With certain matrices, Mat::selfadjoint_eigendecomposition appears to produce incorrect results (faer 0.18.2). In the cases I've looked at it seems like the order of the eigenvalues is not consistent with the order of the eigenvectors.

Below is code that reproduces the bug with the simplest matrix I could find that demonstrates the effect. Changing the value of n makes the issue appear and disappear (e.g. n = 32 and n = 34 seem fine, but not n = 33).

use faer::{Mat, Side, solvers::SolverCore};

fn main() {
    let n = 33;

    let mut a = Mat::zeros(n, n);
    for i in 0 .. n {
        a[(i, i)] = (n - i) as f64;
    }

    let eig = a.selfadjoint_eigendecomposition(Side::Upper);
    let a_reconstructed = eig.reconstruct();

    for i in 0 .. n {
        println!("{:20} {:20} {:20}",
            a[(i, i)],
            a_reconstructed[(i, i)],
            a_reconstructed[(i, i)] - a[(i, i)]);
    }
}
sarah-ek commented 3 months ago

thanks a lot for finding this! i pushed a commit to main that should fix the issue

jbncode commented 3 months ago

Thank you!