Skip to content

Commit a775479

Browse files
slivingstonslivingston
authored
Merge pull request #1234 from @kangwonlee
Add scipy fallbacks for generalized Lyapunov and discrete Sylvester equations
2 parents c7ea9c9 + 2dabe33 commit a775479

3 files changed

Lines changed: 360 additions & 28 deletions

File tree

‎benchmarks/scipy_fallback_bench.py‎

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
# scipy_fallback_bench.py - benchmarks for the SLICOT-free (scipy) fallbacks
2+
# KL, 1 Jul 2026
3+
#
4+
# This benchmark compares the pure scipy/numpy fallbacks against the SLICOT
5+
# (slycot) implementations for the matrix-equation routines that gained a
6+
# ``method`` argument:
7+
#
8+
# * generalized continuous Lyapunov lyap(A, Q, E=E)
9+
# * generalized discrete Lyapunov dlyap(A, Q, E=E)
10+
# * discrete Sylvester (Stein) dlyap(A, Q, C)
11+
#
12+
# The ``time_*`` methods time each (routine, size, method) combination. The
13+
# ``track_*`` method records the accuracy of the generalized-Lyapunov solution
14+
# as a function of cond(E). Every problem is constructed from a known solution
15+
# ``X`` so that both speed and accuracy are measured against ground truth; the
16+
# ``setup`` methods therefore build the matrices *outside* the timed region.
17+
#
18+
# When slycot is not installed the ``method='slycot'`` parameterizations are
19+
# skipped (asv treats NotImplementedError raised in setup() as "skip"), so the
20+
# suite runs with or without slycot.
21+
#
22+
# A single deterministic seed is used per problem, so runs are reproducible and
23+
# comparable across commits. (The tables discussed in PR #1234 were medians
24+
# over several seeds; the ratios here match, as asv's repeated sampling
25+
# averages the timing.)
26+
#
27+
# Run, e.g.:
28+
#
29+
# PYTHONPATH=`pwd` asv run --python=python --bench scipy_fallback
30+
#
31+
# or, since these are plain classes, call the methods directly to reproduce the
32+
# numbers without asv.
33+
34+
import numpy as np
35+
36+
import control as ct
37+
38+
# Fixed seed: deterministic, reproducible problems across runs and commits.
39+
SEED = 20260627
40+
41+
42+
def _slycot_available():
43+
try:
44+
return ct.slycot_check()
45+
except Exception:
46+
return False
47+
48+
49+
def _spd(rng, n):
50+
"""Return a symmetric positive-definite n-by-n matrix."""
51+
P = rng.standard_normal((n, n))
52+
return P @ P.T + n * np.eye(n)
53+
54+
55+
def _make_gen_cont_lyap(rng, n):
56+
# A X E' + E X A' + Q = 0, built from a known SPD solution X (A Hurwitz).
57+
E = np.eye(n) + 0.1 * rng.standard_normal((n, n))
58+
M = rng.standard_normal((n, n))
59+
S = M - (np.linalg.norm(M, 2) + 1.0) * np.eye(n)
60+
A = E @ S
61+
X = _spd(rng, n)
62+
Q = -(A @ X @ E.T + E @ X @ A.T)
63+
Q = 0.5 * (Q + Q.T)
64+
return ct.lyap, (A, Q), dict(E=E), X
65+
66+
67+
def _make_gen_disc_lyap(rng, n):
68+
# A X A' - E X E' + Q = 0, built from a known SPD solution X (A Schur).
69+
E = np.eye(n) + 0.1 * rng.standard_normal((n, n))
70+
M = rng.standard_normal((n, n))
71+
S = M / (np.linalg.norm(M, 2) + 1.0)
72+
A = E @ S
73+
X = _spd(rng, n)
74+
Q = -(A @ X @ A.T - E @ X @ E.T)
75+
Q = 0.5 * (Q + Q.T)
76+
return ct.dlyap, (A, Q), dict(E=E), X
77+
78+
79+
def _make_disc_sylvester(rng, n):
80+
# A X Q' - X + C = 0 (discrete Sylvester / Stein), from a known X.
81+
MA = rng.standard_normal((n, n))
82+
MQ = rng.standard_normal((n, n))
83+
A = MA / (np.linalg.norm(MA, 2) + 1.0)
84+
Q = MQ / (np.linalg.norm(MQ, 2) + 1.0)
85+
X = rng.standard_normal((n, n))
86+
C = X - A @ X @ Q.T
87+
return ct.dlyap, (A, Q, C), dict(), X
88+
89+
90+
_MAKERS = {
91+
'gen_cont_lyap': _make_gen_cont_lyap,
92+
'gen_disc_lyap': _make_gen_disc_lyap,
93+
'disc_sylvester': _make_disc_sylvester,
94+
}
95+
96+
97+
class MatrixEquationTiming:
98+
"""Time the scipy fallback against slycot for the ``method=`` routines."""
99+
100+
params = (
101+
['gen_cont_lyap', 'gen_disc_lyap', 'disc_sylvester'],
102+
[10, 50, 100, 200, 400],
103+
['scipy', 'slycot'],
104+
)
105+
param_names = ['routine', 'n', 'method']
106+
timeout = 120
107+
108+
def setup(self, routine, n, method):
109+
if method == 'slycot' and not _slycot_available():
110+
raise NotImplementedError("slycot not available")
111+
rng = np.random.default_rng(SEED)
112+
self.func, self.args, self.kwargs, X = _MAKERS[routine](rng, n)
113+
# Confirm the method actually solves the problem before timing it.
114+
Xhat = self.func(*self.args, method=method, **self.kwargs)
115+
relerr = np.linalg.norm(Xhat - X, 'fro') / np.linalg.norm(X, 'fro')
116+
assert relerr < 1e-6, f"{routine} {method} n={n}: relerr={relerr:.1e}"
117+
118+
def time_solve(self, routine, n, method):
119+
self.func(*self.args, method=method, **self.kwargs)
120+
121+
122+
class GenLyapAccuracy:
123+
"""Track generalized continuous Lyapunov accuracy versus cond(E).
124+
125+
Both the scipy and slycot paths require E nonsingular and degrade together
126+
as E becomes ill-conditioned (the problem is itself about cond(E)**2
127+
conditioned); this benchmark records that, rather than timing.
128+
"""
129+
130+
params = (
131+
[1e0, 1e2, 1e4, 1e6, 1e8, 1e10, 1e12],
132+
['scipy', 'slycot'],
133+
)
134+
param_names = ['cond_E', 'method']
135+
unit = "relative error"
136+
n = 100
137+
138+
def setup(self, cond_E, method):
139+
if method == 'slycot' and not _slycot_available():
140+
raise NotImplementedError("slycot not available")
141+
n = self.n
142+
rng = np.random.default_rng(SEED)
143+
U, _ = np.linalg.qr(rng.standard_normal((n, n)))
144+
V, _ = np.linalg.qr(rng.standard_normal((n, n)))
145+
M = rng.standard_normal((n, n))
146+
S = M - (np.linalg.norm(M, 2) + 1.0) * np.eye(n)
147+
X = _spd(rng, n)
148+
# E with prescribed condition number: singular values spanning cond_E.
149+
E = (U * np.logspace(0, -np.log10(cond_E), n)) @ V.T
150+
A = E @ S
151+
Q = -(A @ X @ E.T + E @ X @ A.T)
152+
self.A, self.Q, self.E, self.X = A, 0.5 * (Q + Q.T), E, X
153+
154+
def track_relerr(self, cond_E, method):
155+
import warnings
156+
with warnings.catch_warnings():
157+
# Ill-conditioned E deliberately triggers the accuracy warning.
158+
warnings.simplefilter("ignore")
159+
Xhat = ct.lyap(self.A, self.Q, E=self.E, method=method)
160+
return float(np.linalg.norm(Xhat - self.X, 'fro')
161+
/ np.linalg.norm(self.X, 'fro'))

‎control/mateqn.py‎

Lines changed: 143 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -47,15 +47,35 @@ def sb03md(n, C, A, U, dico, job='X', fact='N', trana='N', ldwork=None):
4747
try:
4848
from slycot import sb04qd
4949
except ImportError:
50-
sb0qmd = None
50+
sb04qd = None
5151

5252
try:
5353
from slycot import sg03ad
5454
except ImportError:
55-
sb04ad = None
55+
sg03ad = None
5656

5757
__all__ = ['lyap', 'dlyap', 'dare', 'care']
5858

59+
60+
def _warn_ill_conditioned_E(E):
61+
"""Warn that an ill-conditioned E costs accuracy.
62+
63+
The scipy generalized-Lyapunov fallback reduces the problem to a
64+
standard Lyapunov equation by inverting E, so a poorly conditioned E
65+
costs accuracy (continuous and discrete paths alike, regardless of
66+
whether the underlying scipy solve happens to warn). The generalized
67+
Lyapunov problem is itself ill-conditioned (about cond(E)**2) when E
68+
is, so method='slycot', though it does not form inv(E) explicitly, is
69+
not measurably more accurate in that regime.
70+
"""
71+
condE = np.linalg.cond(E)
72+
if condE > 1.0 / np.sqrt(finfo(float).eps):
73+
warnings.warn(
74+
f"E is ill-conditioned (cond(E) = {condE:.2g}); the generalized "
75+
"Lyapunov solution may have reduced accuracy. The problem itself "
76+
"is ill-conditioned for such E, so method='slycot' is not "
77+
"measurably more accurate.", UserWarning, stacklevel=3)
78+
5979
#
6080
# Lyapunov equation solvers lyap and dlyap
6181
#
@@ -103,6 +123,25 @@ def lyap(A, Q, C=None, E=None, method=None):
103123
X : 2D array
104124
Solution to the Lyapunov or Sylvester equation.
105125
126+
Notes
127+
-----
128+
For the generalized Lyapunov equation, method='slycot' uses the
129+
SLICOT routine SG03AD, based on the generalized Schur method of
130+
Penzl [1]_, which factors the matrix pencil without inverting E.
131+
With method='scipy', the equation is transformed to a standard
132+
Lyapunov equation by inverting E, which requires E to be nonsingular
133+
and loses accuracy when E is ill-conditioned (a UserWarning is then
134+
issued). The generalized Lyapunov problem is itself ill-conditioned
135+
(about cond(E)**2) when E is, so method='slycot', though it does not
136+
invert E, is not measurably more accurate in that case. Both methods
137+
require E nonsingular; a truly singular
138+
(descriptor) E is not currently handled by either.
139+
140+
References
141+
----------
142+
.. [1] Penzl, T., "Numerical solution of generalized Lyapunov
143+
equations", Advances in Computational Mathematics, 8:33-48, 1998.
144+
106145
"""
107146
# Decide what method to use
108147
method = _slycot_or_scipy(method)
@@ -162,8 +201,25 @@ def lyap(A, Q, C=None, E=None, method=None):
162201
_check_shape(E, n, n, square=True, name="E")
163202

164203
if method == 'scipy':
165-
raise ControlArgument(
166-
"method='scipy' not valid for generalized Lyapunov equation")
204+
# Transform to a standard Lyapunov equation by multiplying
205+
# from the left by inv(E) and from the right by inv(E).T:
206+
#
207+
# (E^-1 A) X + X (E^-1 A)^T + E^-1 Q E^-T = 0
208+
#
209+
# This requires E to be nonsingular. SG03AD (method='slycot',
210+
# Penzl's generalized Schur method) factors the pencil without
211+
# inverting E, but a truly singular E is not handled by either
212+
# method.
213+
try:
214+
At = solve(E, A)
215+
Qt = solve(E, solve(E, Q).T).T
216+
except np.linalg.LinAlgError:
217+
raise ControlArgument(
218+
"method='scipy' requires E to be nonsingular; "
219+
"a truly singular E (descriptor system) is not "
220+
"supported by either method")
221+
_warn_ill_conditioned_E(E)
222+
return sp.linalg.solve_continuous_lyapunov(At, -Qt)
167223

168224
# Make sure we have access to the write Slycot routine
169225
try:
@@ -229,6 +285,36 @@ def dlyap(A, Q, C=None, E=None, method=None):
229285
X : 2D array (or matrix)
230286
Solution to the Lyapunov or Sylvester equation.
231287
288+
Notes
289+
-----
290+
For the generalized Lyapunov equation, method='slycot' uses the
291+
SLICOT routine SG03AD, based on the generalized Schur method of
292+
Penzl [1]_, which factors the matrix pencil without inverting E.
293+
With method='scipy', the equation is transformed to a standard
294+
Lyapunov equation by inverting E, which requires E to be nonsingular
295+
and loses accuracy when E is ill-conditioned (a UserWarning is then
296+
issued). The generalized Lyapunov problem is itself ill-conditioned
297+
(about cond(E)**2) when E is, so method='slycot', though it does not
298+
invert E, is not measurably more accurate in that case. Both methods
299+
require E nonsingular; a truly singular
300+
(descriptor) E is not currently handled by either.
301+
302+
For the Sylvester equation, method='slycot' uses the
303+
Hessenberg-Schur method of the SLICOT routine SB04QD [2]_ and
304+
method='scipy' uses the Bartels-Stewart method [3]_; both reduce the
305+
coefficient matrices to (Hessenberg-)Schur form and solve the result
306+
by back-substitution, with O(n^3 + m^3) cost.
307+
308+
References
309+
----------
310+
.. [1] Penzl, T., "Numerical solution of generalized Lyapunov
311+
equations", Advances in Computational Mathematics, 8:33-48, 1998.
312+
.. [2] Golub, G.H., Nash, S., and Van Loan, C., "A Hessenberg-Schur
313+
method for the problem AX + XB = C", IEEE Trans. Automatic
314+
Control, AC-24, pp. 909-913, 1979.
315+
.. [3] Bartels, R.H. and Stewart, G.W., "Solution of the matrix
316+
equation AX + XB = C", Comm. ACM, 15(9), pp. 820-826, 1972.
317+
232318
"""
233319
# Decide what method to use
234320
method = _slycot_or_scipy(method)
@@ -279,8 +365,40 @@ def dlyap(A, Q, C=None, E=None, method=None):
279365
_check_shape(C, n, m, name="C")
280366

281367
if method == 'scipy':
282-
raise ControlArgument(
283-
"method='scipy' not valid for Sylvester equation")
368+
# Solve the discrete-time Sylvester equation
369+
#
370+
# A X Q^T - X + C = 0
371+
#
372+
# by the Bartels-Stewart method, matching the complexity of
373+
# the Hessenberg-Schur algorithm of the SLICOT routine
374+
# SB04QD used by method='slycot' (Golub, Nash, and Van
375+
# Loan, 1979): with complex Schur forms A = U Ta U^H and
376+
# Q^T = V Tq V^H and Y = U^H X V, the transformed equation
377+
# Ta Y Tq - Y + U^H C V = 0 is solved column by column,
378+
# each column requiring one triangular solve. O(n^3 + m^3)
379+
# flops overall.
380+
Ta, U = sp.linalg.schur(A, output='complex')
381+
Tq, V = sp.linalg.schur(Q.T, output='complex')
382+
Ct = U.conj().T @ C @ V
383+
# Solvability requires lam_A * lam_Q != 1 for all pairs of
384+
# eigenvalues (the diagonals of the triangular factors)
385+
if np.min(np.abs(np.outer(np.diag(Tq), np.diag(Ta)) - 1.)) \
386+
< finfo(float).eps * max(
387+
1., np.abs(np.diag(Ta)).max()
388+
* np.abs(np.diag(Tq)).max()):
389+
raise ControlArgument(
390+
"A and Q have a pair of eigenvalues whose product "
391+
"is (almost) equal to 1; the discrete-time "
392+
"Sylvester equation is singular")
393+
Y = np.empty((n, m), dtype=complex)
394+
TaY = np.empty((n, m), dtype=complex) # running Ta @ Y
395+
In = np.eye(n)
396+
for k in range(m):
397+
rhs = -Ct[:, k] - TaY[:, :k] @ Tq[:k, k]
398+
Y[:, k] = sp.linalg.solve_triangular(
399+
Tq[k, k] * Ta - In, rhs)
400+
TaY[:, k] = Ta @ Y[:, k]
401+
return np.real(U @ Y @ V.conj().T)
284402

285403
# Solve the Sylvester equation by calling Slycot function sb04qd
286404
X = sb04qd(n, m, -A, Q.T, C)
@@ -292,8 +410,25 @@ def dlyap(A, Q, C=None, E=None, method=None):
292410
_check_shape(E, n, n, square=True, name="E")
293411

294412
if method == 'scipy':
295-
raise ControlArgument(
296-
"method='scipy' not valid for generalized Lyapunov equation")
413+
# Transform to a standard Lyapunov equation by multiplying
414+
# from the left by inv(E) and from the right by inv(E).T:
415+
#
416+
# (E^-1 A) X (E^-1 A)^T - X + E^-1 Q E^-T = 0
417+
#
418+
# This requires E to be nonsingular. SG03AD (method='slycot',
419+
# Penzl's generalized Schur method) factors the pencil without
420+
# inverting E, but a truly singular E is not handled by either
421+
# method.
422+
try:
423+
At = solve(E, A)
424+
Qt = solve(E, solve(E, Q).T).T
425+
except np.linalg.LinAlgError:
426+
raise ControlArgument(
427+
"method='scipy' requires E to be nonsingular; "
428+
"a truly singular E (descriptor system) is not "
429+
"supported by either method")
430+
_warn_ill_conditioned_E(E)
431+
return sp.linalg.solve_discrete_lyapunov(At, Qt)
297432

298433
# Solve the generalized Lyapunov equation by calling Slycot
299434
# function sg03ad

0 commit comments

Comments
 (0)