Add scipy fallbacks for generalized Lyapunov and discrete Sylvester equations by kwlee2025cpp · Pull Request #1234 · python-control/python-control

kwlee2025cpp · 2026-06-15T10:54:07Z

Summary

lyap and dlyap currently raise ControlArgument("method='scipy' not valid for ...") for three problems, so solving them requires the optional
slycot package:

the generalized continuous Lyapunov equation A X Eᵀ + E X Aᵀ + Q = 0,
the generalized discrete Lyapunov equation A X Aᵀ − E X Eᵀ + Q = 0,
the discrete-time Sylvester equation A X Qᵀ − X + C = 0.

This PR implements pure scipy/numpy fallbacks for all three, so they work
when slycot is not installed. When slycot is installed it remains the
default (method='slycot'), and the new tests cross-validate the scipy
results against it.

What it adds

Generalized Lyapunov (continuous and discrete). The equation is reduced
to a standard Lyapunov equation by a congruence transform with inv(E)
(Ã = E⁻¹A, Q̃ = E⁻¹Q E⁻ᵀ, leaving X unchanged), then solved with
scipy.linalg.solve_{continuous,discrete}_lyapunov. This requires E to be
nonsingular. SLICOT SG03AD (method='slycot'), based on the generalized
Schur method of Penzl, additionally handles singular E; that case raises a
clear ControlArgument pointing the user to method='slycot'.

Discrete-time Sylvester. Solved by the Bartels–Stewart method: complex
Schur factorizations of A and Qᵀ, then a column-by-column back
substitution. This matches the O(n³ + m³) complexity of the
Hessenberg–Schur algorithm in SLICOT SB04QD (method='slycot').
Singularity (a pair of eigenvalues with λᵢ·λⱼ ≈ 1) is detected and raises
a clear ControlArgument.

Ill-conditioned E. Because the generalized-Lyapunov fallback inverts
E, accuracy degrades when E is poorly conditioned. A UserWarning is
now issued when cond(E) > 1/√eps, recommending method='slycot'. (This
also surfaces a case in the continuous path that was previously silently
inaccurate.)

Fixes included

Import-alias bug in the optional-slycot guards: the except ImportError branches bound the sentinels under the wrong names
(sb0qmd for sb04qd, sb04ad for sg03ad), so on a slycot-less
install the intended None fallbacks were never set under the names the
code references.
Stale dlyap docstring note (implied a Kronecker-product method);
corrected to Bartels–Stewart, with Golub–Nash–Van Loan and
Bartels–Stewart references added to the Notes/References.

Testing

scipy↔slycot cross-validation added to test_lyap_g, test_dlyap_g, and
test_dlyap_sylvester (scipy result matches slycot to ~1e-9…1e-16).
New test_lyap_g_ill_conditioned_E_scipy (parametrized over lyap/dlyap)
asserts the ill-conditioned-E UserWarning.
New test_dlyap_sylvester_singular_scipy asserts the singular-pencil
ControlArgument.
Green both with and without slycot: 19 passed, 12 skipped (no slycot) /
31 passed (slycot).

Only control/mateqn.py and control/tests/mateqn_test.py are touched.

Scope

This is an intentionally small, self-contained slice. It is the first of a
short series factoring a larger "scipy fallbacks for slycot-gated routines"
effort into independently reviewable PRs; the matrix-equation solvers go
first because they are foundational and easy to validate against SLICOT.
gram is the planned next slice.

Alternatives considered

The same inv(E) reduction is used by pyMOR's dense scipy Lyapunov
backend. Factoring the pencil (A, E) with scipy.linalg.qz instead of
inverting E (as harold does) is more robust for ill-conditioned E
and is a natural follow-up. Truly singular E (descriptor systems with
infinite generalized eigenvalues) is out of scope here and continues to
require method='slycot'; a pure-scipy path for that case would need
spectral-projector deflation of the infinite spectrum.

…quations 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>

slivingston · 2026-06-16T07:12:06Z

@kwlee2025cpp Thanks for proposing these changes. Can you answer a few questions to help with review and consideration for merging?

Is there a prior issue or discussion that provides context and motivation for this addition?
If the feature is already provided when slycot is installed, is it worth the mantenance effort of an alternative implementation that does not require it? For example, how does the performance compare between this PR (the "fallbacks") and those requiring slycot?
Are you using this in your own work and now want to share the result as open source?

I am a little concerned that this PR is just a Claude-generated thing that no one is asking for. I am also concerned we are about to get many more Claude-generated PRs as your description suggestions ("...It is the first of a
short series factoring a larger 'scipy fallbacks for slycot-gated routines' effort into independently reviewable PRs..."). I am happy to be mistaken about that, so please let me know.

bnavigator · 2026-06-16T09:00:18Z

Approving the workflows because the code on first glance looks good and obviously not malicious. No assessment about scientific merit, I fully support @slivingston's inquiry.

coveralls · 2026-06-16T09:17:26Z

coverage: 94.583% (-0.1%) from 94.73% — kangwonlee:scipy-matrix-equation-fallbacks into python-control:main

The generalized-Lyapunov Notes, the scipy-path comments, and the nonsingular-E error messages stated that SLICOT SG03AD "also handles singular E". It does not: SG03AD factors the matrix pencil without inverting E (its advantage for a nonsingular but ill-conditioned E), but a truly singular (descriptor) E returns a degenerate-pair warning and is out of scope for both the scipy and slycot paths. Correct all six spots (continuous lyap and discrete dlyap) to say so. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

kwlee2025cpp · 2026-06-26T11:42:34Z

Thanks @slivingston, and @bnavigator — these are the right questions, and I would
rather answer them squarely than let this sit as an unmotivated patch.

1. Prior issue / motivation. It is the project's own direction. In #27 ("Move
away from SLICOT?") @murrayrm's 2014 plan was to "replace functions that require
SLICOT with python versions … and if there is some advanced function that really
requires SLICOT, we just leave that in as something that requires an optional
package." That issue is closed-completed, but the pattern is in the tree today —
lqr / dlqr / dare already carry this method= slycot-or-scipy dispatch —
and the coverage is still openly requested (#36 asks the same for dare; #1057 raises the
GPL/BSD question). So this extends an accepted line, not a new one.

One more disclosure, a personal thread: in 2017 (as @beachdweller) I was on #27
thinking about translating the SLICOT Fortran myself — the same year #133 added a
no-Slycot test run. That effort and #27's pointed the same way but never met; I
wish I had carried it through then instead of going quiet — belated apologies. The
viable road was not translation but these scipy fallbacks.

2. Worth it? — let me be straight, because a survey moved under me.

Installation is no longer the argument. slycot 0.7.0 now ships pip
wheels across platforms — manylinux/musllinux on x86_64 and aarch64, macOS,
Windows. "Hard to install" is essentially solved; I will not pretend otherwise.
The honest core is completeness. Without slycot, generalized-E Lyapunov and
discrete Sylvester currently raise — so for the many who never install the
optional package, these basic operations do not exist. The scipy path makes them
work in a stock pip install control: the function is there instead of erroring.
That is independent of wheels and of license.
License is a bounded secondary. slycot's wrapper is GPLv2 (the underlying
SLICOT has been BSD-3-Clause since 2020; only the wrapper relicensing is stuck,
on all-author consent). A scipy path gives BSD-clean coverage for GPL-averse
distributors — but it is distribution-only and does not change python-control's
own BSD license, so I will not overstate it.
Precedent: @murrayrm merged the same move in scipy-based implementation of ss2tf #1108 (2025), "allow MIMO ss2tf
without requiring slycot."
Performance: compiled SLICOT should win at large n; the scipy path is
LAPACK-backed so it is not slow, but I am not claiming it beats SLICOT — happy to
benchmark.
One correction I owe you: a docstring note here overclaims that SLICOT SG03AD
"also handles singular E." It does not — a truly singular E is out of scope for
both paths — and I have pushed a commit fixing the wording.

3. Am I using it. Not in a production project that calls these specific
functions today, and I will retract part of my own pitch: I would have cited
Fortran-free robotics-grader containers (RoboRacer / F1TENTH on ROS 2), but 0.7.0's
wheels just narrowed that to a niche. The honest drivers are the project's
direction (#27 / #36) and completeness — and, further out, a descriptor-system
(singular-E) case that slycot itself cannot do, which is where the real new value
lies.

On the "Claude-generated" concern. Claude assisted — implementation from one
model, review from another, and the commit trailers say so; I am not hiding it. But
the motivation, the mathematics, and the direction are mine: I have been through
every line — the E ≠ I congruence, the Bartels–Stewart back-substitution and its
λ_i·μ_j ≠ 1 solvability condition, the stability argument — and can defend each. I
stand behind every line. If a reviewer finds one I cannot account for, that is a
bug and I want to hear it.

On "many more PRs." It is a finite, bounded effort, sliced small to make review
easier, not to flood you — I am glad to pace it, hold the rest until this one is
judged, or stop. And knowing the honest and modest case for it, if you would still
rather not carry it, I am genuinely happy to withdraw — I do not want to add
maintenance you did not ask for. No hard feelings either way.

slivingston

@kwlee2025cpp Thanks for your reply. Some of it is unclear, but I think the only part critical for review is

Performance: compiled SLICOT should win at large n; the scipy path is
LAPACK-backed so it is not slow, but I am not claiming it beats SLICOT — happy to
benchmark.

Yes, please benchmark. This will also be useful for the other PRs that you plan to make/submit.

kwlee2025cpp · 2026-06-27T12:14:29Z

Thanks @slivingston — benchmark (scipy vs slycot across n, with residuals) coming shortly.

In the meantime, so I fix the unclear parts rather than guess at them: which sections read as unclear to you? I will tighten or rewrite them — I know the reply ran long.

The generalized-Lyapunov Notes and the _warn_ill_conditioned_E helper (its docstring and warning message) claimed method='slycot' is "preferable" / "more robust" than method='scipy' for ill-conditioned E. Benchmarking (scipy vs slycot, n=100, cond(E) swept 1 to 1e12 over 5 seeds) shows both methods lose accuracy at the same rate -- slycot is in fact 1.3-2x *less* accurate for cond(E) >= 1e4, never better. The generalized Lyapunov problem is itself ill-conditioned (about cond(E)^2) when E is, so no method beats that floor. Reword to say both degrade and slycot is not measurably more accurate, and drop the false "use method='slycot' for a more robust solution" line from the warning. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

kwlee2025cpp · 2026-06-29T15:01:56Z

Benchmark — scipy vs slycot. Apple M-series, macOS, Python 3.13, OpenBLAS, slycot 0.7.0. Each problem is built from a known solution; time is the median over seeds (5 for n ≤ 200, 3 at 400, 1 at 800), ≥ 2–3 runs each. Ratio < 1 ⇒ scipy faster.

generalized continuous Lyapunov — `lyap(A, Q, E=E)`

n	scipy (ms)	slycot (ms)	scipy/slycot
10	0.09	0.07	1.27
20	0.19	0.26	0.75
50	1.0	4.6	0.22
100	5.8	22.8	0.26
200	20.9	161	0.13
400	90.7	1003	0.09
800	797	6780	0.12

generalized discrete Lyapunov — `dlyap(A, Q, E=E)`

n	scipy (ms)	slycot (ms)	scipy/slycot
10	0.15	0.09	1.72
20	0.24	0.32	0.75
50	1.6	4.5	0.36
100	6.7	24.8	0.27
200	23.5	181	0.13
400	100	1109	0.09
800	908	7215	0.13

discrete Sylvester — `dlyap(A, Q, C)`

n	scipy (ms)	slycot (ms)	scipy/slycot
10	0.27	0.03	9.78
20	0.68	0.10	6.68
50	8.6	1.7	5.14
100	45.2	10.2	4.42
200	116	56.9	2.05
400	598	352	1.70
800	5187	2728	1.90

Residuals are at machine precision (~1e-15…1e-14) for both methods on all three routines; scipy and slycot agree to ~1e-14 (≤ 1e-9 at n = 800, reflecting the solution's sensitivity at that size).

Headline — I was wrong that SLICOT wins at large n; it depends on the routine:

Generalized Lyapunov (continuous and discrete): scipy is faster, ~8–11× at n ≥ 400. The scipy path reduces via inv(E) to a standard Lyapunov equation (LAPACK Bartels–Stewart); SG03AD solves the generalized equation via generalized Schur, which costs more.
Discrete Sylvester: slycot is faster, ~1.7–10× (widest at small n). My Bartels–Stewart Schur-factors both A and Qᵀ plus a Python column loop; SB04QD's Hessenberg–Schur in Fortran is leaner.

Accuracy vs cond(E) — lyap, n = 100, median over 5 seeds, err = ‖X̂ − X_true‖ / ‖X_true‖:

cond(E)	scipy err	slycot err	slycot / scipy
1e0	7e-15	7e-15	0.99
1e2	1.1e-13	8.5e-14	0.80
1e4	2.6e-10	3.4e-10	1.31
1e6	1.2e-6	1.7e-6	1.61
1e8	7.6e-3	1.5e-2	1.95
1e10	67	110	1.71
1e12	5.2e5	9.2e5	1.82

This one corrected me on my own docstring. I had written that method='slycot' is "preferable" for ill-conditioned E. It isn't: both methods degrade together — the generalized Lyapunov problem is itself ill-conditioned (about cond(E)²) when E is — and slycot is in fact 1.3–2× less accurate for cond(E) ≥ 1e4, never better. I've pushed a commit fixing the docstrings and the runtime warning accordingly.

Net: the fallback is correct everywhere and as accurate as slycot across the whole conditioning range; for generalized Lyapunov it's also a speedup, while for discrete Sylvester slycot remains the faster path when it's installed. (Happy to add the benchmark script to the PR or a gist if useful.)

slivingston · 2026-06-30T06:23:40Z

@kwlee2025cpp Thanks for comparing performance. In this PR, can you add the code for replicating the experiment? I think the best path is in benchmarks/. Give the file a name like scipy_fallback_bench.py. In the future (after this PR is merged), we can expand the file for comparing other fallback routines with SLICOT.

Add benchmarks/scipy_fallback_bench.py (asv) comparing the pure-scipy fallbacks with slycot for generalized Lyapunov (continuous and discrete) and discrete Sylvester, plus an accuracy-vs-cond(E) track for the generalized Lyapunov solution. slycot parameterizations are skipped when slycot is not installed. Requested in the PR python-control#1234 review. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

kwlee2025cpp · 2026-07-02T21:45:32Z

Done — added benchmarks/scipy_fallback_bench.py.

MatrixEquationTiming times each (routine, n, method) for the three routines that gained method= (generalized continuous/discrete Lyapunov and discrete Sylvester). setup() builds each problem from a known solution and asserts the chosen method actually solves it (relerr < 1e-6) before timing.
GenLyapAccuracy is a track_ benchmark recording the generalized-Lyapunov relative error vs cond(E) — the evidence behind the docstring/warning correction in this PR (both paths degrade about cond(E)**2; slycot is not more accurate for ill-conditioned E).
method='slycot' cases raise NotImplementedError in setup when slycot isn't installed, so the suite runs with or without slycot.

Run with:

PYTHONPATH=`pwd` asv run --python=python --bench scipy_fallback

or, since these are plain classes, call the setup/time_/track_ methods directly to reproduce the numbers without asv.

One note on reproducibility: the tables I posted earlier were medians over several seeds, whereas the committed benchmark fixes a single seed per problem so runs are deterministic and comparable across commits (asv's repeated sampling averages the timing). Same generators, same ratios.

I ran the suite locally with asv (existing-environment mode) and it reproduces both the timing and the accuracy-vs-cond(E) tables. Happy to adjust the size grid or fold in more routines as the fallbacks grow.

slivingston · 2026-07-04T04:50:37Z

Thanks @kwlee2025cpp I will review it tomorrow.

kwlee2025cpp force-pushed the scipy-matrix-equation-fallbacks branch from 2c68e3e to c1ab92c Compare June 15, 2026 15:38

slivingston requested changes Jun 27, 2026

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

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

Add scipy fallbacks for generalized Lyapunov and discrete Sylvester equations#1234
kwlee2025cpp wants to merge 4 commits into
python-control:mainfrom
kangwonlee:scipy-matrix-equation-fallbacks

kwlee2025cpp commented Jun 15, 2026

Uh oh!

slivingston commented Jun 16, 2026

Uh oh!

bnavigator commented Jun 16, 2026

Uh oh!

coveralls commented Jun 16, 2026 •

edited

Loading

Uh oh!

kwlee2025cpp commented Jun 26, 2026

Uh oh!

slivingston left a comment

Uh oh!

kwlee2025cpp commented Jun 27, 2026

Uh oh!

kwlee2025cpp commented Jun 29, 2026

Uh oh!

slivingston commented Jun 30, 2026

Uh oh!

kwlee2025cpp commented Jul 2, 2026

Uh oh!

slivingston commented Jul 4, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Uh oh!

Conversation

kwlee2025cpp commented Jun 15, 2026

Summary

What it adds

Fixes included

Testing

Scope

Alternatives considered

Uh oh!

slivingston commented Jun 16, 2026

Uh oh!

bnavigator commented Jun 16, 2026

Uh oh!

coveralls commented Jun 16, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

kwlee2025cpp commented Jun 26, 2026

Uh oh!

slivingston left a comment

Choose a reason for hiding this comment

Uh oh!

kwlee2025cpp commented Jun 27, 2026

Uh oh!

kwlee2025cpp commented Jun 29, 2026

generalized continuous Lyapunov — lyap(A, Q, E=E)

generalized discrete Lyapunov — dlyap(A, Q, E=E)

discrete Sylvester — dlyap(A, Q, C)

Uh oh!

slivingston commented Jun 30, 2026

Uh oh!

kwlee2025cpp commented Jul 2, 2026

Uh oh!

slivingston commented Jul 4, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

coveralls commented Jun 16, 2026 •

edited

Loading

generalized continuous Lyapunov — `lyap(A, Q, E=E)`

generalized discrete Lyapunov — `dlyap(A, Q, E=E)`

discrete Sylvester — `dlyap(A, Q, C)`