Skip to content

fix: Canon modal Schur fallback for ill-conditioned systems#24

Merged
jamestjsp merged 3 commits intomainfrom
fix/canon-modal-schur-fallback
Apr 4, 2026
Merged

fix: Canon modal Schur fallback for ill-conditioned systems#24
jamestjsp merged 3 commits intomainfrom
fix/canon-modal-schur-fallback

Conversation

@jamestjsp
Copy link
Copy Markdown
Owner

@jamestjsp jamestjsp commented Apr 4, 2026 •
edited
Loading

Summary

Hardening fixes from auditing python-control issues and merged PRs for numerical edge cases.

Commit 1 — Canon modal Schur fallback:

Commit 2 — Input validation hardening:

  • StateSpace() rejects improper TFs (num deg > den deg) — python-control#235
  • Kalman/Lqg/H2Syn/HinfSyn reject descriptor systems (E≠nil) — python-control#36
  • Care/Dare validate Q is PSD via shifted Cholesky — python-control#252

Commit 3 — Numerical hardening from PR audit:

  • Canon: replace explicit T⁻¹ with lu.SolveTo() — avoids forming inverse (pc PR#101)
  • c2d matched: fix sign loss in gain fallback — cmplx.Abs ratio → real(contVal/discVal) (pc PR#951)
  • IsStabilizable: add relative tolerance to continuous-time check — roundoff Re=1e-16 no longer falsely unstable (pc PR#345)

Test plan

  • Full test suite passes (all 3 commits)
  • New tests: improper TF rejected, proper/biproper accepted
  • New tests: negative-definite/indefinite Q rejected, PSD Q accepted, non-PD R rejected
  • New tests: descriptor systems rejected in Kalman/Lqg
  • Existing canon tests pass with Solve-based implementation
  • Existing c2d matched tests pass with sign-preserving gain

🤖 Generated with Claude Code

jamestjsp and others added 3 commits April 4, 2026 08:02
@jamestjsp @claude
fix: Canon modal falls back to Schur for ill-conditioned eigenvectors
982472e 
canonModal used eigendecomposition to build the transformation matrix T.
When eigenvalues are nearly repeated (e.g. triple pole at ≈-0.01 from
1e6*s³+30000*s²+300*s+1), eigenvectors become nearly parallel, giving
κ(T)≈8e13. Inverting T amplifies B entries to ~1e14, destroying the
frequency response.

Now checks κ(T)>1e12 (matching eigdecomp.go threshold) and falls back
to ordered real Schur form (DGEES+DTREXC), which is numerically stable
for near-repeated eigenvalues.

Reported via python-control/python-control#1077 — the MIMO system there
triggers the ill-conditioning.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
@jamestjsp @claude
fix: harden input validation (improper TF, descriptor Riccati, Q PSD)
44d7971 
Three defensive checks inspired by python-control edge cases:

1. StateSpace() rejects improper transfer functions (num deg > den deg)
   before conversion, preventing silent wrong results. (pc#235)

2. Kalman, Lqg, H2Syn, HinfSyn reject descriptor systems (E != nil)
   since standard CARE/DARE ignore E, producing wrong results. (pc#36)

3. Care/Dare validate Q is positive semi-definite via Cholesky on
   Q + ε‖Q‖nI. Negative/indefinite Q was silently accepted. (pc#252)
   R positive-definiteness was already caught by Cholesky.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
@jamestjsp @claude
fix: numerical hardening from python-control PR audit
2392746 
Three fixes from auditing merged python-control PRs:

1. canon.go: Replace explicit Tinv computation with lu.SolveTo() in
   both canonModalEig and canonCompanion. Avoids forming the inverse
   matrix, reducing numerical error amplification. (pc PR#101)

2. convert.go: matchedGainFallback used cmplx.Abs ratio which discards
   sign, producing wrong gain for systems with negative DC gain. Now
   uses real(contVal/discVal) matching the non-fallback path. (pc PR#951)

3. ctrbobsv.go: IsStabilizable continuous-time check used zero tolerance
   (real(v) >= 0), so roundoff eigenvalue at Re=1e-16 was falsely
   classified as unstable. Now uses relative tolerance consistent with
   the discrete-time check. (pc PR#345)

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
@jamestjsp
Copy link
Copy Markdown
Owner Author

jamestjsp commented Apr 4, 2026

Benchmark Results (benchstat, 5 runs, Apple M1 Pro)

Execution time (sec/op)

Benchmark Before After Δ Notes
Canon_Modal_N2 2.63µs 2.49µs -5.1% Solve faster than explicit Tinv
Canon_Modal_N10 23.8µs 21.6µs -9.3%
Canon_Modal_N50 755µs 666µs -11.7%
Canon_Modal_N100 4.62ms 8.35ms +80.7%¹ ¹Was returning error (κ=5.6e17); now succeeds via Schur
Canon_Comp_N2 3.10µs 2.95µs -4.8%
Canon_Comp_N10 24.2µs 22.1µs -8.9%
Care_N10 73.0µs 73.5µs +0.8% isPSD check cost
Care_N100 31.2ms 31.3ms +0.6%
Geomean (all) +0.09% Flat

Memory (B/op)

Benchmark Before After Δ Notes
Canon_Modal_N2–N50 -11 to -16% No Tinv allocation needed
Care/Dare_N10 914B 1810B +98% isPSD temp n×n array
Care/Dare_N100 75KB 155KB +106% isPSD temp n×n array

Analysis

  • Canon N2–N50: 5–12% faster, 11–16% less memory from replacing explicit T⁻¹ with lu.SolveTo()
  • Canon N100 (+80%): Not a regression — old code returned ErrSingularTransform for this system (κ(T) = 5.6×10¹⁷). New code succeeds via Schur fallback. The "slowdown" is real work vs fast failure.
  • Care/Dare (+0.6% time, +2× B/op): Cost of isPSD(Q) validation — one temporary n×n Cholesky per call
  • Everything else: within noise

@jamestjsp jamestjsp merged commit 11d7737 into main Apr 4, 2026
2 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

No reviews

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

1 participant

@jamestjsp