dimforge / dimforge/nalgebra

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

Open
#1,109 4 comments 1 reaction 0 assignees View on GitHub
bug P-high
Dominant language
Rust
Stars
4.8k
Forks
565
PR merge metrics
No merged PRs in 30d

Description

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.

```rust
#[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
);
}
```

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.