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

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

1. Is there a prior issue or discussion that provides context and motivation for this addition?
2. 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?
3. 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.

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.

coverage: 94.757% (+0.01%) from 94.743% — kangwonlee:scipy-matrix-equation-fallbacks into python-control:main

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.

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.

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)

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

discrete Sylvester — dlyap(A, Q, C)

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‖:

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.)

@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.

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.

Thanks @kwlee2025cpp I will review it tomorrow.

Hi @slivingston — no rush at all, just a friendly check-in on this one.

The benchmark is in as benchmarks/scipy_fallback_bench.py and reproduces both the timing and the accuracy-vs-cond(E) tables.

In case the red CI is what's giving pause: those failures are pre-existing and unrelated to this PR — test_pole_zero_subplots and test_nyquist_basic fail only on Py3.12 from newer matplotlib/numpy (main was last green on 2026-06-09 and hasn't re-run since). The mateqn changes here pass across the matrix, with and without slycot.

Happy to rebase or adjust whenever you have a chance — thanks for taking a look.

Hi @slivingston — no rush at all, just a friendly check-in on this one.

Thanks for checking in. Sorry I have been slow to reply. I will be done soon.

@kwlee2025cpp Almost done with review. Quick question: Can you tell me the relationship between your account (@kwlee2025cpp ) and @kangwonlee ?

I ask because the accounts have similar names and one has forks of the other's repositories.

Both accounts are mine.

@kangwonlee is my long-standing personal account, and it owns the fork this PR is branched from. @kwlee2025cpp is the account I use for work where I am getting help from an LLM (Claude). I set it up to be upfront about that, rather than let such work blend into my personal account's history. The commits here are authored under my own name and email, and the pushes use @kwlee2025cpp's SSH key, which has write access to that fork — that is why one account appears to hold forks of the other's repositories.

That is the same disclosure I made earlier in this thread: I use Claude's help, I verify the mathematics myself, and I stand behind every line.

Happy to answer anything further. I would rather not re-open this PR from the other account, since that would lose your review thread, but I can consolidate for future contributions if you prefer.

FYI: CI failures have been addressed in #1243. Once that has been merged into main, you should be able to rebase and clear the errors.

lgtm, after a rebase

Rebased onto main, which now includes #1243 — thanks @murrayrm for the CI fix and @sawyerbfuller for the review.

The four commits replay unchanged; the patch is byte-identical to the pre-rebase one. Locally on numpy 2.4.4 / scipy 1.17.1 / matplotlib 3.10.8, mateqn_test.py (31 passed), kwargs_test.py (66 passed), and both tests that were red before #1243 — ctrlplot_test.py::test_pole_zero_subplots and nyquist_test.py::test_nyquist_basic — all pass.

The workflow runs on the new head are showing as action_required, so they may need an approval click before CI reports.

@slivingston Do you have any pending comments on this PR? If not, I'll merge as soon as CI checks clear.

Both answered, briefly:

1. Yes — @beachdweller is an earlier account of mine. That 2017 comment on Move away from SLICOT? #27 was me, back when I was trying to get SLICOT working through f2py/f2c and did not get far. It predates the @kangwonlee / @kwlee2025cpp split and has nothing to do with it.

2. "robotics-grader" was my own shorthand, sorry — it is not a package. I teach RoboRacer / F1TENTH labs and grade them with autograder containers running ROS 2 on ARM (Jetson-class) hardware, so a Fortran-free dependency chain used to matter for the image build. slycot 0.7.0's wheels largely removed that motivation, which is why I called it a niche rather than a selling point. To be precise: those labs use ct.dlqr, so they motivate the general goal rather than this PR's specific functions.

Thank you for the review and for merging. The feedback on warning messages and code comments is welcome — happy to prepare that as a small follow-up whenever suits you.