python-control · kwlee2025cpp · Jun 13, 2026 · Jun 26, 2026 · Jun 29, 2026 · Jul 2, 2026
diff --git a/benchmarks/scipy_fallback_bench.py 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
@@ -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