Incorrect selfadjoint_eigendecomposition for some matrices

[wisp3rwind]

Describe the bug I'm encountering incorrect results of selfadjoint_eigenvalues and selfadjoint_eigendecomposition. This is sensitive to the row/column ordering and the size of the matrix, as seen in the example below, where results change after re-ordering rows and columns.

To Reproduce Steps to reproduce the behavior:

Reproducer

```rust use faer::{ prelude::*, ComplexField, perm::PermRef, dyn_stack::{GlobalPodBuffer, PodStack}, reborrow::ReborrowMut, }; use faer_ext::IntoNdarray; use ndarray_linalg::{EigValsh, UPLO}; use itertools::iproduct; fn make_perm(n: usize, x: usize, y: usize, z: usize) -> faer::perm::Perm { let forward = iproduct!(0..n, 0..3) .map(|(i, a)| match a { 0 => 3 * i + x, 1 => 3 * i + y, 2 => 3 * i + z, _ => unreachable!(), }) .collect::>(); let inverse = forward.clone(); faer::perm::Perm::new_checked( forward, inverse, ) } fn permute_mat_in_place( mut m: MatMut<'_, E>, perm: PermRef<'_, I> ) { let req = faer::perm::permute_cols_in_place_req::(m.nrows(), m.ncols()) .or(faer::perm::permute_rows_in_place_req::(m.nrows(), m.ncols())); let mut buf = GlobalPodBuffer::new(req.unwrap()); let mut stack = PodStack::new(&mut buf); faer::perm::permute_cols_in_place(m.as_mut(), perm, stack.rb_mut()); faer::perm::permute_rows_in_place(m.as_mut(), perm, stack.rb_mut()); } pub fn main() { let n = 6; // Build example matrix, which has a similar 3x3 block structure as // encountered in the actual problem I'm solving let mut mx = Mat::zeros(3 * n, 3 * n); for j in 0..n { for b in 0..3 { for i in 0..n { let d = (i as f64 - j as f64).abs(); let f = 1.0 / (1.0 + (6.0 * (1.0 - d / 8.0)).exp()); mx[(3 * i + b, 3 * j + b)] = if b == 0 { -48.0 * f } else { 24.0 * f }; } mx[(3 * j + b, 3 * j + b)] = 137.0; } } // Build mz, my by permuting rows & columns of mx, which should preserve the spectrum let perm_xy = make_perm(n, 1, 0, 2); let perm_xz = make_perm(n, 2, 1, 0); let mut my = mx.to_owned(); permute_mat_in_place(my.as_mut(), perm_xy.as_ref()); let mut mz = mx.to_owned(); permute_mat_in_place(mz.as_mut(), perm_xz.as_ref()); // Compute eigenvalues using faer let mut ex = mx.selfadjoint_eigenvalues(faer::Side::Lower); ex.sort_by(|a, b| a.partial_cmp(b).unwrap()); let mut ey = my.selfadjoint_eigenvalues(faer::Side::Lower); ey.sort_by(|a, b| a.partial_cmp(b).unwrap()); let mut ez = mz.selfadjoint_eigenvalues(faer::Side::Lower); ez.sort_by(|a, b| a.partial_cmp(b).unwrap()); // Compute eigenvalues using ndarray_linalg let exn = mx.as_ref().into_ndarray().eigvalsh(UPLO::Lower).unwrap(); let eyn = my.as_ref().into_ndarray().eigvalsh(UPLO::Lower).unwrap(); let ezn = mz.as_ref().into_ndarray().eigvalsh(UPLO::Lower).unwrap(); println!("x"); dbg!(ex); dbg!(exn.as_slice().unwrap()); println!("y"); dbg!(ey); dbg!(eyn.as_slice().unwrap()); println!("z"); dbg!(ez); dbg!(ezn.as_slice().unwrap()); } ```

Output

``` x [src/bin/faer_debug.rs:81:5] ex = [ 134.35681418431935, 135.177783778977, 135.71494078360752, 136.83815974377552, 136.83815974534346, 136.90812033100278, 136.9081205343555, 136.91034135407693, 136.91034135407693, 136.91430263263163, 136.91916394454478, 137.1666010493402, 137.17949833298053, 137.18352329745693, 137.31880708244452, 138.70927332950305, 138.7947565653805, 139.25129195618356, ] [src/bin/faer_debug.rs:82:5] exn.as_slice().unwrap() = [ 130.56674808283776, 134.20758866645116, 134.20758866645124, 136.8381597394066, 136.83815973940665, 136.90812016196418, 136.90812016196426, 136.9103413540768, 136.91034135407685, 136.9191641195197, 136.91916411951976, 137.1616717609603, 137.1793172918463, 137.18375967607145, 137.32368052118662, 140.2166259585811, 140.21662595858112, 142.58482266709748, ] y [src/bin/faer_debug.rs:84:5] ey = [ 130.56674808283785, 134.2075886664512, 134.20758866645127, 136.8381597394067, 136.83815973940676, 136.90812016196435, 136.90812016196452, 136.91034135407682, 136.91034135407693, 136.91916411951973, 136.91916411951996, 137.16167176096053, 137.1793172918464, 137.18375967607173, 137.32368052118667, 140.21662595858103, 140.21662595858115, 142.5848226670974, ] [src/bin/faer_debug.rs:85:5] eyn.as_slice().unwrap() = [ 130.5667480828378, 134.20758866645127, 134.2075886664513, 136.8381597394067, 136.8381597394067, 136.90812016196423, 136.90812016196426, 136.9103413540767, 136.91034135407682, 136.91916411951982, 136.91916411951982, 137.16167176096044, 137.17931729184647, 137.18375967607162, 137.3236805211866, 140.21662595858106, 140.21662595858115, 142.58482266709746, ] z [src/bin/faer_debug.rs:87:5] ez = [ 130.56674808283793, 134.2075886664514, 134.20758866645144, 136.83815973940685, 136.83815973940696, 136.90812016196423, 136.90812016196435, 136.91034135407696, 136.91034135407702, 136.9191641195197, 136.91916411951988, 137.16167176096025, 137.17931729184622, 137.18375967607116, 137.32368052118647, 140.21662595858115, 140.21662595858123, 142.58482266709768, ] [src/bin/faer_debug.rs:88:5] ezn.as_slice().unwrap() = [ 130.56674808283768, 134.20758866645124, 134.20758866645133, 136.83815973940668, 136.83815973940673, 136.9081201619642, 136.90812016196432, 136.91034135407682, 136.91034135407685, 136.91916411951982, 136.91916411951982, 137.1616717609604, 137.1793172918463, 137.18375967607145, 137.32368052118653, 140.21662595858106, 140.21662595858112, 142.58482266709743, ] ```

Observe that the eigenvalues are not the same in all cases (for what it's worth, it appears that their sum (i.e. the trace of the matrix) is in fact identical in all cases).

Expected behavior Eigenvalues are identical in all cases.

Desktop (please complete the following information):

• OS: Fedora Linux
• Version: faer 0.19.0, Rust 1.80.0-nightly

[sarah-ek]

that looks like a bug, i'll look into it this weekend

[sarah-ek]

found the issue, seems to be due a comparison where im checking if values can be truncated to zero, and instead of comparing e < eps * sqrt(a * b), i was doing e < eps * sqrt(a + b)

[sarah-ek]

published a fix in 0.19.1 on crates.io. lemme know if this fixes the issue on your side

[wisp3rwind]

Yes, looks like this fixes both the minimized reproducer as well as the original code. Thanks for the quick fix!