scipy/linalg/tests/test_decomp.py::TestSchur::test_sort test failure

[mgorny]

scipy/linalg/tests/test_decomp.py::TestSchur::test_sort test fails in 1.7.1 (IIRC 1.7.0 too).

Error message:

_________________________________________________________ TestSchur.test_sort _________________________________________________________
scipy/linalg/tests/test_decomp.py:1918: in test_sort
    assert_array_almost_equal([[0.1134, 0.5436, 0.8316, 0.],
E   AssertionError: 
E   Arrays are not almost equal to 3 decimals
E   
E   Mismatched elements: 8 / 16 (50%)
E   Max absolute difference: 1.64897635
E   Max relative difference: 2.00019569
E    x: array([[ 0.113,  0.544,  0.832,  0.   ],
E          [-0.113, -0.825,  0.554,  0.   ],
E          [-0.821,  0.131,  0.026, -0.555],
E          [-0.547,  0.087,  0.018,  0.832]])
E    y: array([[-1.134e-01, -5.436e-01,  8.316e-01, -5.470e-15],
E          [ 1.134e-01,  8.245e-01,  5.544e-01, -4.318e-15],
E          [ 8.213e-01, -1.308e-01,  2.650e-02, -5.547e-01],
E          [ 5.475e-01, -8.718e-02,  1.767e-02,  8.321e-01]])
        a          = [[4.0, 3.0, 1.0, -1.0], [-4.5, -3.5, -1.0, 1.0], [9.0, 6.0, -4.0, 4.5], [6.0, 4.0, -3.0, 3.5]]
        s          = array([[-1.41421356,  0.14557215, 11.58158585,  7.71743518],
       [ 0.        , -0.5       , -9.44719662,  0.7184405...    [ 0.        ,  0.        ,  1.41421356, -0.14557215],
       [ 0.        ,  0.        ,  0.        ,  0.5       ]])
        sdim       = 2
        self       = <scipy.linalg.tests.test_decomp.TestSchur object at 0x7f88fcebe9d0>
        u          = array([[-1.13395338e-01, -5.43632173e-01,  8.31628257e-01,
        -5.47030881e-15],
       [ 1.13395338e-01,  8.24476...33e-02,
        -5.54700196e-01],
       [ 5.47521126e-01, -8.71829392e-02,  1.76652155e-02,
         8.32050294e-01]])

Scipy/Numpy/Python version information:

$  python3.8 -c 'import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)'
1.7.1 1.21.1 sys.version_info(major=3, minor=8, micro=11, releaselevel='final', serial=0)

[tupui]

Thanks for the report @mgorny. On which platform are you? I cannot reproduce this on macOS with master.

[mgorny]

Gentoo Linux, amd64.

I suppose this may be relevant:

lapack_opt_info:
system_info:
  NOT AVAILABLE

flame_info:
customize UnixCCompiler
  libraries flame not found in ['/usr/local/lib64', '/usr/local/lib', '/usr/lib64', '/usr/lib', '/usr/lib/']
  NOT AVAILABLE

lapack_info:
  FOUND:
    libraries = ['lapack']
    library_dirs = ['/usr/lib64']
    language = f77

blas_opt_info:
blas_info:
C compiler: x86_64-pc-linux-gnu-gcc-11.2.0 -march=znver2 --param l1-cache-size=32 --param l1-cache-line-size=64 -O2 -pipe -frecord-gcc-switches -fPIC

Fortran f77 compiler: x86_64-pc-linux-gnu-gfortran-11.2.0 -Wall -g -ffixed-form -fno-second-underscore -fPIC -march=znver2 --param l1-cache-size=32 --param l1-cache-line-size=64 -O2 -pipe -frecord-gcc-switches -fPIC -fallow-argument-mismatch
Fortran f90 compiler: x86_64-pc-linux-gnu-gfortran-11.2.0 -Wall -g -fno-second-underscore -fPIC -march=znver2 --param l1-cache-size=32 --param l1-cache-line-size=64 -O2 -pipe -frecord-gcc-switches -fPIC -fallow-argument-mismatch
Fortran fix compiler: x86_64-pc-linux-gnu-gfortran-11.2.0 -Wall -g -ffixed-form -fno-second-underscore -Wall -g -fno-second-underscore -fPIC -march=znver2 --param l1-cache-size=32 --param l1-cache-line-size=64 -O2 -pipe -frecord-gcc-switches -fPIC -fallow-argument-mismatch

[tupui]

How are you installing SciPy? Are you following this guide https://scipy.github.io/devdocs/dev/contributor/quickstart_ubuntu.html#quickstart-ubuntu?

[mgorny]

I'm trying to package the new version for Gentoo. It's based on the ebuild for 1.6.3:
https://gitweb.gentoo.org/repo/gentoo.git/tree/dev-python/scipy/scipy-1.6.3.ebuild

minus rm cythonize.dat || die line that I've removed.

[mgorny]

Ok, I think I've just accidentally discovered an important data point: the tests fail also on 1.6.3, and they did not fail before. Last time I've built scipy 1.6.3 was on 10th of July but i'm not sure if I've run tests then. I'm 90% sure this is due to some other package being upgraded…

[tupui]

FYI as for the infrastructure goes, we added support for Pythran and some wrapper around boost.

[tupui]

@mgorny do you still have an issue here?

[mgorny]

@mgorny do you still have an issue here?

I'm sorry, are you asking me to test git master, or retest the versions where I could reproduce it previously?

[tupui]

I am asking if you still have the problem. In the end did you manage to package SciPy? We did not change anything regarding this specific problem but you might have done something on your side and found a workaround.

[mgorny]

I've just tested 1.6.3 again and it still fails, so I presume 1.7.1 also still fails. To be honest, I've really forgotten about this problem (my plate's pretty full at the moment), and I'm afraid I still have no idea how to address it.

[ilayn]

For some reason the first two columns flipped signs. Annoying indeed but not a blocker at least. Might depend on a multitude of stuff.

[mgorny]

Ok, I've spent the whole day trying to figure this out and it seems to be caused by upgrade of lapack from 3.9.0 to 3.10.0.

[mgorny]

FTR, scipy/linalg/tests/test_solvers.py::test_solve_discrete_are failure is also caused by lapack-3.10.0.

[drew-parsons]

We're seeing this error (failure in TestSchur.test_sort) in the Debian build of scipy 1.7.1 for armhf,
https://ci.debian.net/data/autopkgtest/testing/armhf/s/scipy/16604596/log.gz

[rgommers]

This has started to fail in the macOS Meson CI job consistently now (see for example https://github.com/scipy/scipy/runs/7772683544?check_suite_focus=true), and I bet that it's due to OpenBLAS 0.3.21 becoming available on conda-forge (that CI job is the only one using conda-forge, and it does pull in 0.3.21).

[WarrenWeckesser]

It looks like the problem is actually the test. The Schur decomposition is not unique. As @ilayn pointed out in a previous comment, in the failing test, the signs of the first two columns of the unitary factor are flipped. One can see that some of the signs in the triangular factor are also flipped, and it turns out that these sign changes are consistent with the sign changes in the unitary factor.

Here's the test calculation when it passes:

In [25]: np.set_printoptions(suppress=True)
In [26]: from scipy.linalg import schur

In [27]: a = [[4., 3., 1., -1.],
    ...:      [-4.5, -3.5, -1., 1.],
    ...:      [9., 6., -4., 4.5],
    ...:      [6., 4., -3., 3.5]]

In [28]: s, u, sdim = schur(a, sort='lhp')

In [29]: s
Out[29]: 
array([[ -1.41421356,   0.14557215, -11.58158585,  -7.71743518],
       [  0.        ,  -0.5       ,   9.44719662,  -0.71844057],
       [  0.        ,   0.        ,   1.41421356,  -0.14557215],
       [  0.        ,   0.        ,   0.        ,   0.5       ]])

In [30]: u
Out[30]: 
array([[ 0.11339534,  0.54363217,  0.83162826, -0.        ],
       [-0.11339534, -0.82447635,  0.55441884, -0.        ],
       [-0.82128169,  0.13077441,  0.02649782, -0.5547002 ],
       [-0.54752113,  0.08718294,  0.01766522,  0.83205029]])

E is a reflection of the first two columns when it multiplies a matrix on the right:

In [31]: E = np.diag([-1, -1, 1, 1])

Compute another Schur decomposition from u and s; u2 and s2 look like the solutions computed in the failing test:

In [32]: u2 = u @ E

In [33]: s2 = E @ s @ E

Verify that this is a correct Schur decomposition:

In [34]: u2 @ s2 @ u2.T  # Should be `a`
Out[34]: 
array([[ 4. ,  3. ,  1. , -1. ],
       [-4.5, -3.5, -1. ,  1. ],
       [ 9. ,  6. , -4. ,  4.5],
       [ 6. ,  4. , -3. ,  3.5]])

In [35]: u2
Out[35]: 
array([[-0.11339534, -0.54363217,  0.83162826, -0.        ],
       [ 0.11339534,  0.82447635,  0.55441884, -0.        ],
       [ 0.82128169, -0.13077441,  0.02649782, -0.5547002 ],
       [ 0.54752113, -0.08718294,  0.01766522,  0.83205029]])

In [36]: s2
Out[36]: 
array([[-1.41421356,  0.14557215, 11.58158585,  7.71743518],
       [ 0.        , -0.5       , -9.44719662,  0.71844057],
       [ 0.        ,  0.        ,  1.41421356, -0.14557215],
       [ 0.        ,  0.        ,  0.        ,  0.5       ]])

Check that u2 is unitary:

In [38]: u2 @ u2.T
Out[38]: 
array([[ 1., -0.,  0., -0.],
       [-0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [-0.,  0.,  0.,  1.]])

So we have to fix the test.