Repository: scipy/scipy
File: scipy/linalg/tests/test_solvers.py (Python)

Code before:
        assert_allclose(res, np.zeros_like(res), atol=atol)

    def test_infeasible(self):
        # An infeasible example taken from https://arxiv.org/abs/1505.04861v1
        A = np.triu(np.ones((3, 3)))
        A[0, 1] = -1
        B = np.array([[1, 1, 0], [0, 0, 1]]).T
        Q = np.full_like(A, -2) + np.diag([8, -1, -1.9])
        R = np.diag([-10, 0.1])
        assert_raises(LinAlgError, solve_continuous_are, A, B, Q, R)


def test_solve_generalized_continuous_are():
    cases = [
        # Two random examples differ by s term
        # in the absence of any literature for demanding examples.
        (np.array([[2.769230e-01, 8.234578e-01, 9.502220e-01],
                   [4.617139e-02, 6.948286e-01, 3.444608e-02],
                   [9.713178e-02, 3.170995e-01, 4.387444e-01]]),
         np.array([[3.815585e-01, 1.868726e-01],
                   [7.655168e-01, 4.897644e-01],
                   [7.951999e-01, 4.455862e-01]]),
         np.eye(3),
         np.eye(2),
         np.array([[6.463130e-01, 2.760251e-01, 1.626117e-01],
                   [7.093648e-01, 6.797027e-01, 1.189977e-01],
                   [7.546867e-01, 6.550980e-01, 4.983641e-01]]),
         np.zeros((3, 2)),
         None),
        (np.array([[2.769230e-01, 8.234578e-01, 9.502220e-01],
                   [4.617139e-02, 6.948286e-01, 3.444608e-02],
                   [9.713178e-02, 3.170995e-01, 4.387444e-01]]),
         np.array([[3.815585e-01, 1.868726e-01],
                   [7.655168e-01, 4.897644e-01],
                   [7.951999e-01, 4.455862e-01]]),
         np.eye(3),
         np.eye(2),
         np.array([[6.463130e-01, 2.760251e-01, 1.626117e-01],
                   [7.093648e-01, 6.797027e-01, 1.189977e-01],
                   [7.546867e-01, 6.550980e-01, 4.983641e-01]]),
         np.ones((3, 2)),
         None)
        ]

    min_decimal = (10, 10)

    def _test_factory(case, dec):
        """Checks if X = A'XA-(A'XB)(R+B'XB)^-1(B'XA)+Q) is true"""
        a, b, q, r, e, s, knownfailure = case
        if knownfailure:
            pytest.xfail(reason=knownfailure)

Review comment (bug):
This test does indeed fail when the source patch is reversed, and in a manner consistent with the original ticket, with `ValueError: Matrix e should be square.` So, that's a good sign. It is also good that no old test cases were modified, only new cases were added.

So, mechanically it looks pretty good to me. The CI failure is just https://github.com/scipy/scipy/issues/23689. That said, I've not been super active in the batching work apart from complaining about error handling things, so may make sense for someone else to give it a quick look.

Matt--how severe is this? Does it justify, on its own, another bug fix release for example? (assuming `1.17.0` is not coming out until early January 2026)

Code after:
        assert_allclose(res, np.zeros_like(res), atol=atol)

    def test_infeasible(self):
        # An infeasible example taken from https://arxiv.org/abs/1505.04861v1
        A = np.triu(np.ones((3, 3)))
        A[0, 1] = -1
        B = np.array([[1, 1, 0], [0, 0, 1]]).T
        Q = np.full_like(A, -2) + np.diag([8, -1, -1.9])
        R = np.diag([-10, 0.1])
        assert_raises(LinAlgError, solve_continuous_are, A, B, Q, R)


class TestSolveCommonAre:
    @pytest.mark.parametrize("solver", [solve_continuous_are, solve_discrete_are])
    def test_with_skipped_array_argument_gh23336(self, solver):
        # gh-23336 reported a failure when optional argument `e` was skipped
        A = np.array([[-0.9, 0.25], [0, -1.1]])
        B = np.array([[0.23], [0.45]])
        Q = np.eye(2)
        R = np.atleast_2d(0.45)
        E = np.eye(2)
        S = np.array([[0.1], [0.2]])

        res = solver(A, B, Q, R, s=S)
        ref = solver(A, B, Q, R, E, S)
        np.testing.assert_allclose(res, ref)


def test_solve_generalized_continuous_are():
    cases = [
        # Two random examples differ by s term
        # in the absence of any literature for demanding examples.
        (np.array([[2.769230e-01, 8.234578e-01, 9.502220e-01],
                   [4.617139e-02, 6.948286e-01, 3.444608e-02],
                   [9.713178e-02, 3.170995e-01, 4.387444e-01]]),
         np.array([[3.815585e-01, 1.868726e-01],
                   [7.655168e-01, 4.897644e-01],
                   [7.951999e-01, 4.455862e-01]]),
         np.eye(3),
         np.eye(2),
         np.array([[6.463130e-01, 2.760251e-01, 1.626117e-01],
                   [7.093648e-01, 6.797027e-01, 1.189977e-01],
                   [7.546867e-01, 6.550980e-01, 4.983641e-01]]),
         np.zeros((3, 2)),
         None),
        (np.array([[2.769230e-01, 8.234578e-01, 9.502220e-01],
                   [4.617139e-02, 6.948286e-01, 3.444608e-02],
                   [9.713178e-02, 3.170995e-01, 4.387444e-01]]),
         np.array([[3.815585e-01, 1.868726e-01],
                   [7.655168e-01, 4.897644e-01],
                   [7.951999e-01, 4.455862e-01]]),