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

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)

About this issue

Original URL
State: closed
Created 3 years ago
Comments: 16 (16 by maintainers)

Links to this issue

Mailman 3 [SciPy-Dev] ANN: SciPy 1.9.1 - SciPy-Dev - python.org

Commits related to this issue

BUG/TST: linalg: Check the results of 'schur' more carefully. The Schur decomposition is not unique, so testing the result of linalg.schur against hard-coded matrices is not reliable. The updated tes... — committed to WarrenWeckesser/scipy by WarrenWeckesser 2 years ago
BUG/TST: linalg: Check the results of 'schur' more carefully. The Schur decomposition is not unique, so testing the result of linalg.schur against hard-coded matrices is not reliable. The updated tes... — committed to tartopohm/scipy by WarrenWeckesser 2 years ago
BUG/TST: linalg: Check the results of 'schur' more carefully. The Schur decomposition is not unique, so testing the result of linalg.schur against hard-coded matrices is not reliable. The updated tes... — committed to tylerjereddy/scipy by WarrenWeckesser 2 years ago

Most upvoted comments

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.

WarrenWeckesser on Aug 10, 2022