Issue with Schur decomposition?

[lesurp]

Hi,

Code

    let p = nalgebra::Matrix2::new(2.5000002, -1.4999999, -1.4999999, 2.4999998);
    let schur = p.schur();
    println!("{}", schur.unpack().0);
    println!("{}", schur.complex_eigenvalues());
    match schur.eigenvalues() {
        Some(d) => println!("{}", d),
        None => println!("None"),
    }

Result

  ┌                                         ┐
  │  0.7071068283270011  0.7071067340460907 │
  │ -0.7071067340460907  0.7071068283270011 │
  └                                         ┘

  ┌                         ┐
  │ 2.4999999999999996+NaNi │
  │ 2.4999999999999996+NaNi │
  └                         ┘

None

Expected

An equivalent cpp code:

#include <Eigen/Eigen>
#include <Eigen/src/Eigenvalues/RealSchur.h>
#include <iostream>
#include <iomanip>

int main(int argc, char *argv[])
{
    Eigen::Matrix2d p;
    p << 2.5000002, -1.4999999, -1.4999999, 2.4999998;

    auto schur = Eigen::RealSchur<Eigen::Matrix2d>(p);
    std::cout << std::setprecision(17);
    std::cout << schur.matrixU() << std::endl;
    std::cout << schur.matrixT() << std::endl;
}

produces:

 0.70710682832700111  0.70710673404609059
-0.70710673404609059  0.70710682832700111

3.9999999000000122                  0
                 0 1.0000000999999865

If I replace the 2.* with an exact 2.5 the schur decomposition is correct. Note that the computed eigenvectors are correct - only the computed eigenvalues are wrong.

Am I misusing this decomposition? Of course for a 22 I don't need* to use schur so it's not exactly a problem but still...

[daingun]

Hi,

I'm having a similar problem calculating the eigenvalues of a 2x2 matrix. I think that both issues are the same.

let a = DMatrix::from_row_slice(2, 2, &[-2., 0., 3., -7.]);
let schur = Schur::new(a);
dbg!(&schur);
let poles = schur.complex_eigenvalues();
dbg!(poles);

Here the eigenvalues must be equal to -2. and -7., however the result is Complex {re: -4.499999999999999, im: NaN} for both eigenvalues.

The issue is that Schur decomposition for 2x2 matrix is performed calculating the eigenvalues inside the function compute_2x2_eigvals, whose result is correct, than the initial matrix is rotated to get the t matrix. Here the precision is lost and the following matrix is returned

t: Matrix {
    data: VecStorage {
        data: [
            -1.9999999999999996,
            -0.0000000000000004440892098500626,
            -2.999999999999999,
            -6.999999999999999,
        ],
    ...
}

When complex_eigenvalues() is called the computation is delegated to do_complex_eigenvalues which checks for exact zero equality the lower left element of the t matrix, since it fails, it assumes that the eigenvalues are complex. The calculation of the eigenvalues is differs from the one in compute_2x2_eigvals, the discriminant is positive so the sqrt returns NaN.

I think that matrix rotations to get t should preserve the zero in the lower left corner of the 2x2 matrix and the precision of eigenvalues computed in the first place.

[sebcrozet]

Closing this as it should be fixed by #668. Please reopen if this problem persists.