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

Open
kwlee2025cpp wants to merge 1 commit into
base: main
Choose a base branch
from
+185 −28
Open

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
1 commit
Select commit
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
137 changes: 129 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,33 @@ 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 the inv(E) congruence transform will lose 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). SLICOT sg03ad
(method='slycot') avoids inverting E and is preferable in that case.
"""
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 "
"method='scipy' generalized Lyapunov solution may have reduced "
"accuracy. Use method='slycot' (SLICOT sg03ad) for a more "
"robust solution.", UserWarning, stacklevel=3)

#
# Lyapunov equation solvers lyap and dlyap
#
Expand Down Expand Up @@ -103,6 +121,21 @@ 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 also handles singular 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); method='slycot' does
not invert E and is preferable in that case.

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 +195,23 @@ 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; the SLICOT routine
# SG03AD used by method='slycot' (based on the generalized
# Schur method of Penzl (1998)) also handles singular E.
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; "
"use method='slycot' (SLICOT sg03ad) for singular E")
_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 +277,32 @@ 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 also handles singular 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); method='slycot' does
not invert E and is preferable in that case.

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 +353,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)
Expand All @@ -292,8 +398,23 @@ 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; the SLICOT routine
# SG03AD used by method='slycot' (based on the generalized
# Schur method of Penzl (1998)) also handles singular E.
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; "
"use method='slycot' (SLICOT sg03ad) for singular E")
_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
76 changes: 56 additions & 20 deletions control/tests/mateqn_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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',
Expand All @@ -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',
Expand Down Expand Up @@ -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]])
Expand Down
Loading