diff --git a/benchmarks/scipy_fallback_bench.py b/benchmarks/scipy_fallback_bench.py new file mode 100644 index 000000000..7934a3b5d --- /dev/null +++ b/benchmarks/scipy_fallback_bench.py @@ -0,0 +1,161 @@ +# scipy_fallback_bench.py - benchmarks for the SLICOT-free (scipy) fallbacks +# KL, 1 Jul 2026 +# +# This benchmark compares the pure scipy/numpy fallbacks against the SLICOT +# (slycot) implementations for the matrix-equation routines that gained a +# ``method`` argument: +# +# * generalized continuous Lyapunov lyap(A, Q, E=E) +# * generalized discrete Lyapunov dlyap(A, Q, E=E) +# * discrete Sylvester (Stein) dlyap(A, Q, C) +# +# The ``time_*`` methods time each (routine, size, method) combination. The +# ``track_*`` method records the accuracy of the generalized-Lyapunov solution +# as a function of cond(E). Every problem is constructed from a known solution +# ``X`` so that both speed and accuracy are measured against ground truth; the +# ``setup`` methods therefore build the matrices *outside* the timed region. +# +# When slycot is not installed the ``method='slycot'`` parameterizations are +# skipped (asv treats NotImplementedError raised in setup() as "skip"), so the +# suite runs with or without slycot. +# +# A single deterministic seed is used per problem, so runs are reproducible and +# comparable across commits. (The tables discussed in PR #1234 were medians +# over several seeds; the ratios here match, as asv's repeated sampling +# averages the timing.) +# +# Run, e.g.: +# +# PYTHONPATH=`pwd` asv run --python=python --bench scipy_fallback +# +# or, since these are plain classes, call the methods directly to reproduce the +# numbers without asv. + +import numpy as np + +import control as ct + +# Fixed seed: deterministic, reproducible problems across runs and commits. +SEED = 20260627 + + +def _slycot_available(): + try: + return ct.slycot_check() + except Exception: + return False + + +def _spd(rng, n): + """Return a symmetric positive-definite n-by-n matrix.""" + P = rng.standard_normal((n, n)) + return P @ P.T + n * np.eye(n) + + +def _make_gen_cont_lyap(rng, n): + # A X E' + E X A' + Q = 0, built from a known SPD solution X (A Hurwitz). + E = np.eye(n) + 0.1 * rng.standard_normal((n, n)) + M = rng.standard_normal((n, n)) + S = M - (np.linalg.norm(M, 2) + 1.0) * np.eye(n) + A = E @ S + X = _spd(rng, n) + Q = -(A @ X @ E.T + E @ X @ A.T) + Q = 0.5 * (Q + Q.T) + return ct.lyap, (A, Q), dict(E=E), X + + +def _make_gen_disc_lyap(rng, n): + # A X A' - E X E' + Q = 0, built from a known SPD solution X (A Schur). + E = np.eye(n) + 0.1 * rng.standard_normal((n, n)) + M = rng.standard_normal((n, n)) + S = M / (np.linalg.norm(M, 2) + 1.0) + A = E @ S + X = _spd(rng, n) + Q = -(A @ X @ A.T - E @ X @ E.T) + Q = 0.5 * (Q + Q.T) + return ct.dlyap, (A, Q), dict(E=E), X + + +def _make_disc_sylvester(rng, n): + # A X Q' - X + C = 0 (discrete Sylvester / Stein), from a known X. + MA = rng.standard_normal((n, n)) + MQ = rng.standard_normal((n, n)) + A = MA / (np.linalg.norm(MA, 2) + 1.0) + Q = MQ / (np.linalg.norm(MQ, 2) + 1.0) + X = rng.standard_normal((n, n)) + C = X - A @ X @ Q.T + return ct.dlyap, (A, Q, C), dict(), X + + +_MAKERS = { + 'gen_cont_lyap': _make_gen_cont_lyap, + 'gen_disc_lyap': _make_gen_disc_lyap, + 'disc_sylvester': _make_disc_sylvester, +} + + +class MatrixEquationTiming: + """Time the scipy fallback against slycot for the ``method=`` routines.""" + + params = ( + ['gen_cont_lyap', 'gen_disc_lyap', 'disc_sylvester'], + [10, 50, 100, 200, 400], + ['scipy', 'slycot'], + ) + param_names = ['routine', 'n', 'method'] + timeout = 120 + + def setup(self, routine, n, method): + if method == 'slycot' and not _slycot_available(): + raise NotImplementedError("slycot not available") + rng = np.random.default_rng(SEED) + self.func, self.args, self.kwargs, X = _MAKERS[routine](rng, n) + # Confirm the method actually solves the problem before timing it. + Xhat = self.func(*self.args, method=method, **self.kwargs) + relerr = np.linalg.norm(Xhat - X, 'fro') / np.linalg.norm(X, 'fro') + assert relerr < 1e-6, f"{routine} {method} n={n}: relerr={relerr:.1e}" + + def time_solve(self, routine, n, method): + self.func(*self.args, method=method, **self.kwargs) + + +class GenLyapAccuracy: + """Track generalized continuous Lyapunov accuracy versus cond(E). + + Both the scipy and slycot paths require E nonsingular and degrade together + as E becomes ill-conditioned (the problem is itself about cond(E)**2 + conditioned); this benchmark records that, rather than timing. + """ + + params = ( + [1e0, 1e2, 1e4, 1e6, 1e8, 1e10, 1e12], + ['scipy', 'slycot'], + ) + param_names = ['cond_E', 'method'] + unit = "relative error" + n = 100 + + def setup(self, cond_E, method): + if method == 'slycot' and not _slycot_available(): + raise NotImplementedError("slycot not available") + n = self.n + rng = np.random.default_rng(SEED) + U, _ = np.linalg.qr(rng.standard_normal((n, n))) + V, _ = np.linalg.qr(rng.standard_normal((n, n))) + M = rng.standard_normal((n, n)) + S = M - (np.linalg.norm(M, 2) + 1.0) * np.eye(n) + X = _spd(rng, n) + # E with prescribed condition number: singular values spanning cond_E. + E = (U * np.logspace(0, -np.log10(cond_E), n)) @ V.T + A = E @ S + Q = -(A @ X @ E.T + E @ X @ A.T) + self.A, self.Q, self.E, self.X = A, 0.5 * (Q + Q.T), E, X + + def track_relerr(self, cond_E, method): + import warnings + with warnings.catch_warnings(): + # Ill-conditioned E deliberately triggers the accuracy warning. + warnings.simplefilter("ignore") + Xhat = ct.lyap(self.A, self.Q, E=self.E, method=method) + return float(np.linalg.norm(Xhat - self.X, 'fro') + / np.linalg.norm(self.X, 'fro')) diff --git a/control/mateqn.py b/control/mateqn.py index 9d1349b0c..b23e42def 100644 --- a/control/mateqn.py +++ b/control/mateqn.py @@ -47,15 +47,35 @@ def sb03md(n, C, A, U, dico, job='X', fact='N', trana='N', ldwork=None): try: from slycot import sb04qd except ImportError: - sb0qmd = None + sb04qd = None try: from slycot import sg03ad except ImportError: - sb04ad = None + sg03ad = None __all__ = ['lyap', 'dlyap', 'dare', 'care'] + +def _warn_ill_conditioned_E(E): + """Warn that an ill-conditioned E costs accuracy. + + The scipy generalized-Lyapunov fallback reduces the problem to a + standard Lyapunov equation by inverting E, so a poorly conditioned E + costs accuracy (continuous and discrete paths alike, regardless of + whether the underlying scipy solve happens to warn). The generalized + Lyapunov problem is itself ill-conditioned (about cond(E)**2) when E + is, so method='slycot', though it does not form inv(E) explicitly, is + not measurably more accurate in that regime. + """ + condE = np.linalg.cond(E) + if condE > 1.0 / np.sqrt(finfo(float).eps): + warnings.warn( + f"E is ill-conditioned (cond(E) = {condE:.2g}); the generalized " + "Lyapunov solution may have reduced accuracy. The problem itself " + "is ill-conditioned for such E, so method='slycot' is not " + "measurably more accurate.", UserWarning, stacklevel=3) + # # Lyapunov equation solvers lyap and dlyap # @@ -103,6 +123,25 @@ def lyap(A, Q, C=None, E=None, method=None): X : 2D array Solution to the Lyapunov or Sylvester equation. + Notes + ----- + For the generalized Lyapunov equation, method='slycot' uses the + SLICOT routine SG03AD, based on the generalized Schur method of + Penzl [1]_, which factors the matrix pencil without inverting E. + With method='scipy', the equation is transformed to a standard + Lyapunov equation by inverting E, which requires E to be nonsingular + and loses accuracy when E is ill-conditioned (a UserWarning is then + issued). The generalized Lyapunov problem is itself ill-conditioned + (about cond(E)**2) when E is, so method='slycot', though it does not + invert E, is not measurably more accurate in that case. Both methods + require E nonsingular; a truly singular + (descriptor) E is not currently handled by either. + + References + ---------- + .. [1] Penzl, T., "Numerical solution of generalized Lyapunov + equations", Advances in Computational Mathematics, 8:33-48, 1998. + """ # Decide what method to use method = _slycot_or_scipy(method) @@ -162,8 +201,25 @@ def lyap(A, Q, C=None, E=None, method=None): _check_shape(E, n, n, square=True, name="E") if method == 'scipy': - raise ControlArgument( - "method='scipy' not valid for generalized Lyapunov equation") + # Transform to a standard Lyapunov equation by multiplying + # from the left by inv(E) and from the right by inv(E).T: + # + # (E^-1 A) X + X (E^-1 A)^T + E^-1 Q E^-T = 0 + # + # This requires E to be nonsingular. SG03AD (method='slycot', + # Penzl's generalized Schur method) factors the pencil without + # inverting E, but a truly singular E is not handled by either + # method. + try: + At = solve(E, A) + Qt = solve(E, solve(E, Q).T).T + except np.linalg.LinAlgError: + raise ControlArgument( + "method='scipy' requires E to be nonsingular; " + "a truly singular E (descriptor system) is not " + "supported by either method") + _warn_ill_conditioned_E(E) + return sp.linalg.solve_continuous_lyapunov(At, -Qt) # Make sure we have access to the write Slycot routine try: @@ -229,6 +285,36 @@ def dlyap(A, Q, C=None, E=None, method=None): X : 2D array (or matrix) Solution to the Lyapunov or Sylvester equation. + Notes + ----- + For the generalized Lyapunov equation, method='slycot' uses the + SLICOT routine SG03AD, based on the generalized Schur method of + Penzl [1]_, which factors the matrix pencil without inverting E. + With method='scipy', the equation is transformed to a standard + Lyapunov equation by inverting E, which requires E to be nonsingular + and loses accuracy when E is ill-conditioned (a UserWarning is then + issued). The generalized Lyapunov problem is itself ill-conditioned + (about cond(E)**2) when E is, so method='slycot', though it does not + invert E, is not measurably more accurate in that case. Both methods + require E nonsingular; a truly singular + (descriptor) E is not currently handled by either. + + For the Sylvester equation, method='slycot' uses the + Hessenberg-Schur method of the SLICOT routine SB04QD [2]_ and + method='scipy' uses the Bartels-Stewart method [3]_; both reduce the + coefficient matrices to (Hessenberg-)Schur form and solve the result + by back-substitution, with O(n^3 + m^3) cost. + + References + ---------- + .. [1] Penzl, T., "Numerical solution of generalized Lyapunov + equations", Advances in Computational Mathematics, 8:33-48, 1998. + .. [2] Golub, G.H., Nash, S., and Van Loan, C., "A Hessenberg-Schur + method for the problem AX + XB = C", IEEE Trans. Automatic + Control, AC-24, pp. 909-913, 1979. + .. [3] Bartels, R.H. and Stewart, G.W., "Solution of the matrix + equation AX + XB = C", Comm. ACM, 15(9), pp. 820-826, 1972. + """ # Decide what method to use method = _slycot_or_scipy(method) @@ -279,8 +365,40 @@ def dlyap(A, Q, C=None, E=None, method=None): _check_shape(C, n, m, name="C") if method == 'scipy': - raise ControlArgument( - "method='scipy' not valid for Sylvester equation") + # Solve the discrete-time Sylvester equation + # + # A X Q^T - X + C = 0 + # + # by the Bartels-Stewart method, matching the complexity of + # the Hessenberg-Schur algorithm of the SLICOT routine + # SB04QD used by method='slycot' (Golub, Nash, and Van + # Loan, 1979): with complex Schur forms A = U Ta U^H and + # Q^T = V Tq V^H and Y = U^H X V, the transformed equation + # Ta Y Tq - Y + U^H C V = 0 is solved column by column, + # each column requiring one triangular solve. O(n^3 + m^3) + # flops overall. + Ta, U = sp.linalg.schur(A, output='complex') + Tq, V = sp.linalg.schur(Q.T, output='complex') + Ct = U.conj().T @ C @ V + # Solvability requires lam_A * lam_Q != 1 for all pairs of + # eigenvalues (the diagonals of the triangular factors) + if np.min(np.abs(np.outer(np.diag(Tq), np.diag(Ta)) - 1.)) \ + < finfo(float).eps * max( + 1., np.abs(np.diag(Ta)).max() + * np.abs(np.diag(Tq)).max()): + raise ControlArgument( + "A and Q have a pair of eigenvalues whose product " + "is (almost) equal to 1; the discrete-time " + "Sylvester equation is singular") + Y = np.empty((n, m), dtype=complex) + TaY = np.empty((n, m), dtype=complex) # running Ta @ Y + In = np.eye(n) + for k in range(m): + rhs = -Ct[:, k] - TaY[:, :k] @ Tq[:k, k] + Y[:, k] = sp.linalg.solve_triangular( + Tq[k, k] * Ta - In, rhs) + TaY[:, k] = Ta @ Y[:, k] + return np.real(U @ Y @ V.conj().T) # Solve the Sylvester equation by calling Slycot function sb04qd X = sb04qd(n, m, -A, Q.T, C) @@ -292,8 +410,25 @@ def dlyap(A, Q, C=None, E=None, method=None): _check_shape(E, n, n, square=True, name="E") if method == 'scipy': - raise ControlArgument( - "method='scipy' not valid for generalized Lyapunov equation") + # Transform to a standard Lyapunov equation by multiplying + # from the left by inv(E) and from the right by inv(E).T: + # + # (E^-1 A) X (E^-1 A)^T - X + E^-1 Q E^-T = 0 + # + # This requires E to be nonsingular. SG03AD (method='slycot', + # Penzl's generalized Schur method) factors the pencil without + # inverting E, but a truly singular E is not handled by either + # method. + try: + At = solve(E, A) + Qt = solve(E, solve(E, Q).T).T + except np.linalg.LinAlgError: + raise ControlArgument( + "method='scipy' requires E to be nonsingular; " + "a truly singular E (descriptor system) is not " + "supported by either method") + _warn_ill_conditioned_E(E) + return sp.linalg.solve_discrete_lyapunov(At, Qt) # Solve the generalized Lyapunov equation by calling Slycot # function sg03ad diff --git a/control/tests/mateqn_test.py b/control/tests/mateqn_test.py index 77bf553bf..5cc34c2ab 100644 --- a/control/tests/mateqn_test.py +++ b/control/tests/mateqn_test.py @@ -90,19 +90,22 @@ def test_lyap_sylvester(self, method): X_scipy = lyap(A, B, C, method='scipy') assert_array_almost_equal(X_scipy, X) - @pytest.mark.slycot - def test_lyap_g(self): + @pytest.mark.parametrize('method', + ['scipy', + pytest.param('slycot', marks=pytest.mark.slycot)]) + def test_lyap_g(self, method): A = array([[-1, 2], [-3, -4]]) Q = array([[3, 1], [1, 1]]) E = array([[1, 2], [2, 1]]) - X = lyap(A, Q, None, E) + X = lyap(A, Q, None, E, method=method) # print("The solution obtained is ", X) assert_array_almost_equal(A @ X @ E.T + E @ X @ A.T + Q, zeros((2,2))) - # Make sure that trying to solve with SciPy generates an error - with pytest.raises(ControlArgument, match="'scipy' not valid"): - X = lyap(A, Q, None, E, method='scipy') + # Compare methods + if method == 'slycot': + X_scipy = lyap(A, Q, None, E, method='scipy') + assert_array_almost_equal(X_scipy, X) @pytest.mark.parametrize('method', ['scipy', @@ -125,39 +128,73 @@ def test_dlyap(self, method): X_scipy = dlyap(A,Q, method='scipy') assert_array_almost_equal(X_scipy, X) - @pytest.mark.slycot - def test_dlyap_g(self): + @pytest.mark.parametrize('method', + ['scipy', + pytest.param('slycot', marks=pytest.mark.slycot)]) + def test_dlyap_g(self, method): A = array([[-0.6, 0],[-0.1, -0.4]]) Q = array([[3, 1],[1, 1]]) E = array([[1, 1],[2, 1]]) - X = dlyap(A, Q, None, E) + X = dlyap(A, Q, None, E, method=method) # print("The solution obtained is ", X) assert_array_almost_equal(A @ X @ A.T - E @ X @ E.T + Q, zeros((2,2))) - # Make sure that trying to solve with SciPy generates an error - with pytest.raises(ControlArgument, match="'scipy' not valid"): - X = dlyap(A, Q, None, E, method='scipy') + # Compare methods + if method == 'slycot': + X_scipy = dlyap(A, Q, None, E, method='scipy') + assert_array_almost_equal(X_scipy, X) - @pytest.mark.slycot - def test_dlyap_sylvester(self): + @pytest.mark.parametrize('method', + ['scipy', + pytest.param('slycot', marks=pytest.mark.slycot)]) + def test_dlyap_sylvester(self, method): A = 5 B = array([[4, 3], [4, 3]]) C = array([2, 1]) - X = dlyap(A,B,C) + X = dlyap(A, B, C, method=method) # print("The solution obtained is ", X) assert_array_almost_equal(A * X @ B.T - X + C, zeros((1,2))) A = array([[2, 1], [1, 2]]) B = array([[1, 2], [0.5, 0.1]]) C = array([[1, 0], [0, 1]]) - X = dlyap(A, B, C) + X = dlyap(A, B, C, method=method) # print("The solution obtained is ", X) assert_array_almost_equal(A @ X @ B.T - X + C, zeros((2,2))) - # Make sure that trying to solve with SciPy generates an error - with pytest.raises(ControlArgument, match="'scipy' not valid"): - X = dlyap(A, B, C, method='scipy') + # Compare methods + if method == 'slycot': + X_scipy = dlyap(A, B, C, method='scipy') + assert_array_almost_equal(X_scipy, X) + + @pytest.mark.parametrize("cdlyap", [lyap, dlyap]) + def test_lyap_g_singular_E_scipy(self, cdlyap): + """Generalized Lyapunov with singular E raises on scipy path""" + A = array([[-1, 2], [-3, -4]]) + Q = array([[3, 1], [1, 1]]) + E = array([[1, 2], [2, 4]]) # singular + with pytest.raises(ControlArgument, match="E to be nonsingular"): + cdlyap(A, Q, None, E, method='scipy') + + @pytest.mark.parametrize("cdlyap", [lyap, dlyap]) + def test_lyap_g_ill_conditioned_E_scipy(self, cdlyap): + """Generalized Lyapunov with ill-conditioned E warns on scipy path""" + A = array([[-1, 2], [-3, -4]]) + Q = array([[3, 1], [1, 1]]) + E = array([[1, 1], [1, 1 + 1e-12]]) # nonsingular, cond ~ 4e12 + with pytest.warns(UserWarning, match="ill-conditioned"): + cdlyap(A, Q, None, E, method='scipy') + + def test_dlyap_sylvester_singular_scipy(self): + """Discrete Sylvester with an eigenvalue product of 1 is singular""" + # eig(A) = {2, 3}, eig(Q) = {0.5, 0.7}; the product 2 * 0.5 = 1 + # makes the discrete-time Sylvester operator singular + A = array([[2., 0.], [0., 3.]]) + Q = array([[0.5, 0.], [0., 0.7]]) + C = array([[1., 0.], [0., 1.]]) + with pytest.raises(ControlArgument, match="singular"): + dlyap(A, Q, C, method='scipy') @pytest.mark.parametrize('method', ['scipy', @@ -274,7 +311,6 @@ def test_dare(self, method): lam = eigvals(A - B @ G) assert_array_less(abs(lam), 1.0) - @pytest.mark.slycot def test_dare_compare(self): A = np.array([[-0.6, 0], [-0.1, -0.4]]) Q = np.array([[2, 1], [1, 0]])