Skip to content

Add scipy fallbacks for generalized Lyapunov and discrete Sylvester equations #1234

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
slivingston merged 4 commits into from
Aug 14, 2026
Merged

Add scipy fallbacks for generalized Lyapunov and discrete Sylvester equations #1234

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
4608123
Add scipy fallbacks for generalized Lyapunov and discrete Sylvester e…
kwlee2025cpp Jun 13, 2026
89b2f68
Correct SG03AD singular-E claim in lyap/dlyap docs
kwlee2025cpp Jun 26, 2026
3f6dbac
Correct ill-conditioned-E accuracy claim in lyap/dlyap
kwlee2025cpp Jun 29, 2026
2dabe33
Add benchmarks for scipy vs slycot matrix-equation fallbacks
kwlee2025cpp Jul 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
161 changes: 161 additions & 0 deletions benchmarks/scipy_fallback_bench.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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'))
151 changes: 143 additions & 8 deletions control/mateqn.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 "

@slivingston slivingston Aug 14, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A warning like this has too much scope. The warning can only appear when method='scipy', but it has a claim about method='slycot'. If the underlying SLICOT method is improved and this warning message is no longer true, will we remember to update this warning? Probably not because it is only on method='scipy' code paths. I recommend that either the claim is removed (i.e., the warning does not refer to slycot) or a qualifying word is used, like "method='slycot' may not be more accurate" (so, "may" instead of "is").

"measurably more accurate.", UserWarning, stacklevel=3)

#
# Lyapunov equation solvers lyap and dlyap
#
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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")

@slivingston slivingston Aug 14, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My comment about _warn_ill_conditioned_E() applies to this warning as well. Furthermore, "either method" requires that there are only 2 methods. In the future, there may be more than scipy and slycot methods, at which time this warning message text would need to change, whether or not a method='scipy' code path is affected.

_warn_ill_conditioned_E(E)
return sp.linalg.solve_continuous_lyapunov(At, -Qt)

# Make sure we have access to the write Slycot routine
try:
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

@slivingston slivingston Aug 14, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What does #running Ta @ Y indicate? If it is only a reminder that TaY = np.empty() is to allocate the space before filling it in the loop, then delete the comment because it does not help with reading the code. This purpose is obvious, and a similar call is made for Y but no such comment is there.

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)
Expand All @@ -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
Expand Down
Loading
Loading