Incorrect symmetric eigendecomposition for a 4x4 f64 matrix (swapped eigenvectors)
- 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.