possible bug, hanging on getting eigenvalues?

[hwchen]

Hi,

I’m trying to find the roots of a polynomial, and nalgebra has overall worked great.

Recently, I’ve found that my calculations hang on certain inputs, for which numpy immediately returns a value.

You can see my code for finding roots at https://github.com/hwchen/roots.

A pathological input is in this test case, which you can also see in the repo's unit tests:

   #[test]
    fn hanging() {
        let input = &[
            -770.9700059294701,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            434.3959909081459,
            0.0,
        ];

        let input: Vec<f64> = input.into_iter().rev().cloned().collect();
        dbg!(&input);
        let rs = roots(input.as_slice());
    }

Thanks for any help.

[sebcrozet]

Hi! This is a bug. It seems like the following matrix makes the Schur decomposition loop indefinitely:

  ┌                                                                                                                                                                                               ┐
  │  0 0 0 0 0 0 0 0 0 0 1.7748092111017977 │
  │  1 0 0 0 0 0 0 0 0 0 │
  │  0 1 0 0 0 0 0 0 0 0 │
  │  0 0 1 0 0 0 0 0 0 0 │
  │  0 0 0 1 0 0 0 0 0 0 │
  │  0 0 0 0 1 0 0 0 0 0 │
  │  0 0 0 0 0 1 0 0 0 0 │
  │  0 0 0 0 0 0 1 0 0 0 │
  │  0 0 0 0 0 0 0 1 0 0 │
  │  0 0 0 0 0 0 0 0 1 0 │

For what it's worth, note that you can use Schur::try_new to specify a maximum number of iterations. This does not fix your problem, but at least this gives you a way to avoid the infinite hanging (though the result will be None in those cases).

[hwchen]

Thanks! try_new will get me through for now, looking forward to the bugfix. Let me know if there's anything I can do to help.

[TimsHisProjects]

I can reproduce this issue with multiple matrices. I also tested your solution @hwchen which works great on most occasions. However, there are instances in which the function complex_eigenvalues() returns NaN while eigenvalues() returns correct eigenvalues. One such case is with:

fn test_real_roots_1() {
    dbg!(roots(&vec![-225., 0., 224., 0., 0.]));
}

This yields:

roots(p) = [
    Complex {
        re: 0.0,
        im: NaN
    },
    Complex {
        re: 0.0,
        im: NaN
    },
    Complex {
        re: 0.0,
        im: 0.0
    },
    Complex {
        re: 0.0,
        im: 0.0
    }
]

While just calling eigenvalues() yields:

Matrix {
    data: VecStorage {
        data: [
            0.9977753031397177,
            -0.9977753031397177
        ],
        nrows: Dynamic {
            value: 2
        },
        ncols: U1
    }
}

I might be wrong here but complex_eigenvalues should return all real eigenvalues in complex form right? Thanks for any help.