direct variational solver sometimes yields different sign structure

[marscher]

With recent versions of NumPy/Scipy (1.15 / 1.0.0) one observes failing tests in variational/tests/solvers/tests/test_direct.py like this:

    def test_eig_corr(self):
        C0 = np.array([[1.0, 0.3, 0.2],
                       [0.3, 0.8, 0.5],
                       [0.2, 0.5, 0.9]])
        Ct_sym = np.array([[0.5, 0.1, 0.0],
                           [0.1, 0.3, 0.3],
                           [0.0, 0.3, 0.2]])
        Ct_nonsym = np.array([[0.5, 0.1, 0.3],
                              [0.1, 0.3, 0.3],
                              [0.0, 0.3, 0.2]])
        # reference solution
        import scipy
        for Ct in [Ct_sym, Ct_nonsym]:
            v0, R0 = scipy.linalg.eig(Ct, C0)
            v0, R0 = sort_by_norm_and_imag_sign(v0, R0)
            for method in ['QR', 'schur']:
                # Test correctness
                v, R = direct.eig_corr(C0, Ct, method=method)
                v, R = sort_by_norm_and_imag_sign(v, R)
>               np.testing.assert_allclose(v0, v, err_msg='eig vals not equal for method %s' % method)  # eigenvalues equal?
E               AssertionError: 
E               Not equal to tolerance rtol=1e-07, atol=0
E               eig vals not equal for method QR
E               (mismatch 66.66666666666666%)
E                x: array([ 0.447489-0.061667j,  0.447489+0.061667j, -0.094020+0.j      ])
E                y: array([ 0.447489+0.061667j,  0.447489-0.061667j, -0.094020+0.j      ])

Although in this unit test there is a utility function to sort the eigenvalues by sign and magnitude.

@fnueske Do you think this is an issue with reference solution obtained from SciPy? If so, do you think it is save to ignore the sign of the imaginary value, as it seems only these signs differ.

[fnueske]

Hey @marscher, are we already sorting the eigenvalues separately by real and imaginary parts in this test? It seems the results are fine, only the pair of conjugate eigenvalues is swapped.

[marscher]

On 09/17/2018 12:10 PM, Feliks Nüske wrote:

Hey @marscher https://github.com/marscher, are we already sorting the eigenvalues separately by real and imaginary parts in this test? Yes we are.

It seems the results are fine, only the pair of conjugate eigenvalues is swapped.

So is it sufficient to take the absolute value of the imaginary part into account?

[fnueske]

Yes.

[marscher]

https://github.com/marscher/PyEMMA/commit/239c336eaa9b3db660c7c1115b04f3dcb9139743 @fnueske I check against the complex conjugate in case the sign structure is different - could you please check that this approach is also correct?

[fnueske]

That should do it. Are the tests passing now?

[marscher]

yes they are.