Skip to content

Commit c1ab92c

Browse files
kwlee2025cppclaude
kwlee2025cpp
and
claude
committed
Add scipy fallbacks for generalized Lyapunov and discrete Sylvester equations
lyap and dlyap previously raised ControlArgument for method='scipy' (explicit, or auto-selected when slycot is absent) on three cases that SLICOT handles. Add pure scipy/numpy fallbacks for all three: - Generalized Lyapunov, continuous and discrete (E != I): congruence transform by inv(E) to a standard Lyapunov equation, then scipy's solve_continuous/discrete_lyapunov. Requires E nonsingular; SLICOT sg03ad (Penzl's generalized Schur method) also handles singular E and remains the method='slycot' path. A nonsingular-E failure raises a clear ControlArgument. Because the transform inverts E, an ill-conditioned E yields reduced accuracy; this now emits a UserWarning recommending method='slycot' (the continuous path was previously silent, while the discrete path warned only incidentally). - Discrete Sylvester (A X Q^T - X + C = 0): Bartels-Stewart method via complex Schur factors of A and Q^T with column-by-column triangular solves, O(n^3 + m^3), matching the Hessenberg-Schur cost of SLICOT sb04qd. Includes the solvability check that no eigenvalue pair of A and Q^T has product (almost) equal to 1. Also fix two incorrect slycot import-fallback aliases (sb0qmd -> sb04qd, sb04ad -> sg03ad) so the names resolve to None, not NameError, when slycot is absent; and correct the dlyap Notes, which described the scipy Sylvester path as a Kronecker O((nm)^3) solve (it is now Bartels-Stewart, O(n^3 + m^3)). Tests parametrize method=[None, 'scipy', 'slycot'] and cross-check that the scipy and slycot solutions agree, and cover the nonsingular-E requirement, the ill-conditioned-E warning, and the singular discrete-Sylvester case. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
1 parent 7d5019e commit c1ab92c

2 files changed

Lines changed: 185 additions & 28 deletions

File tree

‎control/mateqn.py‎

Lines changed: 129 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -47,15 +47,33 @@ 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 the inv(E) congruence transform will lose 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). SLICOT sg03ad
67+
(method='slycot') avoids inverting E and is preferable in that case.
68+
"""
69+
condE = np.linalg.cond(E)
70+
if condE > 1.0 / np.sqrt(finfo(float).eps):
71+
warnings.warn(
72+
f"E is ill-conditioned (cond(E) = {condE:.2g}); the "
73+
"method='scipy' generalized Lyapunov solution may have reduced "
74+
"accuracy. Use method='slycot' (SLICOT sg03ad) for a more "
75+
"robust solution.", UserWarning, stacklevel=3)
76+
5977
#
6078
# Lyapunov equation solvers lyap and dlyap
6179
#
@@ -103,6 +121,21 @@ def lyap(A, Q, C=None, E=None, method=None):
103121
X : 2D array
104122
Solution to the Lyapunov or Sylvester equation.
105123
124+
Notes
125+
-----
126+
For the generalized Lyapunov equation, method='slycot' uses the
127+
SLICOT routine SG03AD, based on the generalized Schur method of
128+
Penzl [1]_, which also handles singular E. With method='scipy', the
129+
equation is transformed to a standard Lyapunov equation by inverting
130+
E, which requires E to be nonsingular and loses accuracy when E is
131+
ill-conditioned (a UserWarning is then issued); method='slycot' does
132+
not invert E and is preferable in that case.
133+
134+
References
135+
----------
136+
.. [1] Penzl, T., "Numerical solution of generalized Lyapunov
137+
equations", Advances in Computational Mathematics, 8:33-48, 1998.
138+
106139
"""
107140
# Decide what method to use
108141
method = _slycot_or_scipy(method)
@@ -162,8 +195,23 @@ def lyap(A, Q, C=None, E=None, method=None):
162195
_check_shape(E, n, n, square=True, name="E")
163196

164197
if method == 'scipy':
165-
raise ControlArgument(
166-
"method='scipy' not valid for generalized Lyapunov equation")
198+
# Transform to a standard Lyapunov equation by multiplying
199+
# from the left by inv(E) and from the right by inv(E).T:
200+
#
201+
# (E^-1 A) X + X (E^-1 A)^T + E^-1 Q E^-T = 0
202+
#
203+
# This requires E to be nonsingular; the SLICOT routine
204+
# SG03AD used by method='slycot' (based on the generalized
205+
# Schur method of Penzl (1998)) also handles singular E.
206+
try:
207+
At = solve(E, A)
208+
Qt = solve(E, solve(E, Q).T).T
209+
except np.linalg.LinAlgError:
210+
raise ControlArgument(
211+
"method='scipy' requires E to be nonsingular; "
212+
"use method='slycot' (SLICOT sg03ad) for singular E")
213+
_warn_ill_conditioned_E(E)
214+
return sp.linalg.solve_continuous_lyapunov(At, -Qt)
167215

168216
# Make sure we have access to the write Slycot routine
169217
try:
@@ -229,6 +277,32 @@ def dlyap(A, Q, C=None, E=None, method=None):
229277
X : 2D array (or matrix)
230278
Solution to the Lyapunov or Sylvester equation.
231279
280+
Notes
281+
-----
282+
For the generalized Lyapunov equation, method='slycot' uses the
283+
SLICOT routine SG03AD, based on the generalized Schur method of
284+
Penzl [1]_, which also handles singular E. With method='scipy', the
285+
equation is transformed to a standard Lyapunov equation by inverting
286+
E, which requires E to be nonsingular and loses accuracy when E is
287+
ill-conditioned (a UserWarning is then issued); method='slycot' does
288+
not invert E and is preferable in that case.
289+
290+
For the Sylvester equation, method='slycot' uses the
291+
Hessenberg-Schur method of the SLICOT routine SB04QD [2]_ and
292+
method='scipy' uses the Bartels-Stewart method [3]_; both reduce the
293+
coefficient matrices to (Hessenberg-)Schur form and solve the result
294+
by back-substitution, with O(n^3 + m^3) cost.
295+
296+
References
297+
----------
298+
.. [1] Penzl, T., "Numerical solution of generalized Lyapunov
299+
equations", Advances in Computational Mathematics, 8:33-48, 1998.
300+
.. [2] Golub, G.H., Nash, S., and Van Loan, C., "A Hessenberg-Schur
301+
method for the problem AX + XB = C", IEEE Trans. Automatic
302+
Control, AC-24, pp. 909-913, 1979.
303+
.. [3] Bartels, R.H. and Stewart, G.W., "Solution of the matrix
304+
equation AX + XB = C", Comm. ACM, 15(9), pp. 820-826, 1972.
305+
232306
"""
233307
# Decide what method to use
234308
method = _slycot_or_scipy(method)
@@ -279,8 +353,40 @@ def dlyap(A, Q, C=None, E=None, method=None):
279353
_check_shape(C, n, m, name="C")
280354

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

285391
# Solve the Sylvester equation by calling Slycot function sb04qd
286392
X = sb04qd(n, m, -A, Q.T, C)
@@ -292,8 +398,23 @@ def dlyap(A, Q, C=None, E=None, method=None):
292398
_check_shape(E, n, n, square=True, name="E")
293399

294400
if method == 'scipy':
295-
raise ControlArgument(
296-
"method='scipy' not valid for generalized Lyapunov equation")
401+
# Transform to a standard Lyapunov equation by multiplying
402+
# from the left by inv(E) and from the right by inv(E).T:
403+
#
404+
# (E^-1 A) X (E^-1 A)^T - X + E^-1 Q E^-T = 0
405+
#
406+
# This requires E to be nonsingular; the SLICOT routine
407+
# SG03AD used by method='slycot' (based on the generalized
408+
# Schur method of Penzl (1998)) also handles singular E.
409+
try:
410+
At = solve(E, A)
411+
Qt = solve(E, solve(E, Q).T).T
412+
except np.linalg.LinAlgError:
413+
raise ControlArgument(
414+
"method='scipy' requires E to be nonsingular; "
415+
"use method='slycot' (SLICOT sg03ad) for singular E")
416+
_warn_ill_conditioned_E(E)
417+
return sp.linalg.solve_discrete_lyapunov(At, Qt)
297418

298419
# Solve the generalized Lyapunov equation by calling Slycot
299420
# function sg03ad

‎control/tests/mateqn_test.py‎

Lines changed: 56 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -90,19 +90,22 @@ def test_lyap_sylvester(self, method):
9090
X_scipy = lyap(A, B, C, method='scipy')
9191
assert_array_almost_equal(X_scipy, X)
9292

93-
@pytest.mark.slycot
94-
def test_lyap_g(self):
93+
@pytest.mark.parametrize('method',
94+
['scipy',
95+
pytest.param('slycot', marks=pytest.mark.slycot)])
96+
def test_lyap_g(self, method):
9597
A = array([[-1, 2], [-3, -4]])
9698
Q = array([[3, 1], [1, 1]])
9799
E = array([[1, 2], [2, 1]])
98-
X = lyap(A, Q, None, E)
100+
X = lyap(A, Q, None, E, method=method)
99101
# print("The solution obtained is ", X)
100102
assert_array_almost_equal(A @ X @ E.T + E @ X @ A.T + Q,
101103
zeros((2,2)))
102104

103-
# Make sure that trying to solve with SciPy generates an error
104-
with pytest.raises(ControlArgument, match="'scipy' not valid"):
105-
X = lyap(A, Q, None, E, method='scipy')
105+
# Compare methods
106+
if method == 'slycot':
107+
X_scipy = lyap(A, Q, None, E, method='scipy')
108+
assert_array_almost_equal(X_scipy, X)
106109

107110
@pytest.mark.parametrize('method',
108111
['scipy',
@@ -125,39 +128,73 @@ def test_dlyap(self, method):
125128
X_scipy = dlyap(A,Q, method='scipy')
126129
assert_array_almost_equal(X_scipy, X)
127130

128-
@pytest.mark.slycot
129-
def test_dlyap_g(self):
131+
@pytest.mark.parametrize('method',
132+
['scipy',
133+
pytest.param('slycot', marks=pytest.mark.slycot)])
134+
def test_dlyap_g(self, method):
130135
A = array([[-0.6, 0],[-0.1, -0.4]])
131136
Q = array([[3, 1],[1, 1]])
132137
E = array([[1, 1],[2, 1]])
133-
X = dlyap(A, Q, None, E)
138+
X = dlyap(A, Q, None, E, method=method)
134139
# print("The solution obtained is ", X)
135140
assert_array_almost_equal(A @ X @ A.T - E @ X @ E.T + Q,
136141
zeros((2,2)))
137142

138-
# Make sure that trying to solve with SciPy generates an error
139-
with pytest.raises(ControlArgument, match="'scipy' not valid"):
140-
X = dlyap(A, Q, None, E, method='scipy')
143+
# Compare methods
144+
if method == 'slycot':
145+
X_scipy = dlyap(A, Q, None, E, method='scipy')
146+
assert_array_almost_equal(X_scipy, X)
141147

142-
@pytest.mark.slycot
143-
def test_dlyap_sylvester(self):
148+
@pytest.mark.parametrize('method',
149+
['scipy',
150+
pytest.param('slycot', marks=pytest.mark.slycot)])
151+
def test_dlyap_sylvester(self, method):
144152
A = 5
145153
B = array([[4, 3], [4, 3]])
146154
C = array([2, 1])
147-
X = dlyap(A,B,C)
155+
X = dlyap(A, B, C, method=method)
148156
# print("The solution obtained is ", X)
149157
assert_array_almost_equal(A * X @ B.T - X + C, zeros((1,2)))
150158

151159
A = array([[2, 1], [1, 2]])
152160
B = array([[1, 2], [0.5, 0.1]])
153161
C = array([[1, 0], [0, 1]])
154-
X = dlyap(A, B, C)
162+
X = dlyap(A, B, C, method=method)
155163
# print("The solution obtained is ", X)
156164
assert_array_almost_equal(A @ X @ B.T - X + C, zeros((2,2)))
157165

158-
# Make sure that trying to solve with SciPy generates an error
159-
with pytest.raises(ControlArgument, match="'scipy' not valid"):
160-
X = dlyap(A, B, C, method='scipy')
166+
# Compare methods
167+
if method == 'slycot':
168+
X_scipy = dlyap(A, B, C, method='scipy')
169+
assert_array_almost_equal(X_scipy, X)
170+
171+
@pytest.mark.parametrize("cdlyap", [lyap, dlyap])
172+
def test_lyap_g_singular_E_scipy(self, cdlyap):
173+
"""Generalized Lyapunov with singular E raises on scipy path"""
174+
A = array([[-1, 2], [-3, -4]])
175+
Q = array([[3, 1], [1, 1]])
176+
E = array([[1, 2], [2, 4]]) # singular
177+
with pytest.raises(ControlArgument, match="E to be nonsingular"):
178+
cdlyap(A, Q, None, E, method='scipy')
179+
180+
@pytest.mark.parametrize("cdlyap", [lyap, dlyap])
181+
def test_lyap_g_ill_conditioned_E_scipy(self, cdlyap):
182+
"""Generalized Lyapunov with ill-conditioned E warns on scipy path"""
183+
A = array([[-1, 2], [-3, -4]])
184+
Q = array([[3, 1], [1, 1]])
185+
E = array([[1, 1], [1, 1 + 1e-12]]) # nonsingular, cond ~ 4e12
186+
with pytest.warns(UserWarning, match="ill-conditioned"):
187+
cdlyap(A, Q, None, E, method='scipy')
188+
189+
def test_dlyap_sylvester_singular_scipy(self):
190+
"""Discrete Sylvester with an eigenvalue product of 1 is singular"""
191+
# eig(A) = {2, 3}, eig(Q) = {0.5, 0.7}; the product 2 * 0.5 = 1
192+
# makes the discrete-time Sylvester operator singular
193+
A = array([[2., 0.], [0., 3.]])
194+
Q = array([[0.5, 0.], [0., 0.7]])
195+
C = array([[1., 0.], [0., 1.]])
196+
with pytest.raises(ControlArgument, match="singular"):
197+
dlyap(A, Q, C, method='scipy')
161198

162199
@pytest.mark.parametrize('method',
163200
['scipy',
@@ -274,7 +311,6 @@ def test_dare(self, method):
274311
lam = eigvals(A - B @ G)
275312
assert_array_less(abs(lam), 1.0)
276313

277-
@pytest.mark.slycot
278314
def test_dare_compare(self):
279315
A = np.array([[-0.6, 0], [-0.1, -0.4]])
280316
Q = np.array([[2, 1], [1, 0]])

0 commit comments

Comments
 (0)