[dongxuelian2]

Part of #1081: the matrix relative_gain_array calculation and its numerical tolerances.

The exact det(A) != 0 branch can switch between an inverse and a generalized inverse after a change of units. For example, the third column below is the first plus twice the second:

using ControlSystemsBase
A = [1.0 8 17; -4 6 8; -8 1 -6]
relative_gain_array(A)
relative_gain_array(1e-3A)

On master 9baa912a, the largest absolute entries are approximately 0.843 and 6.92e15. Both are 0.843 after this change, agreeing with an independent 256-bit UC reference. The Float32 generalized path also promotes its result to Float64 and can interpret Float32 normalization noise as an independent direction: scaling [1 2 3;2 4 6;3 6 9] by 1e-12 gives entries around 524288 instead of 1/9.

There is a separate range failure even with a finite result: [1e300 1e-300;1e300 2e-300] gives [1 -1;0 2] on master. Its correct RGA is [2 -1;-1 2], which the new calculation recovers.

Simply always using geometric normalization/SVD is insufficient. [1 1e-100 .1;1e-100 1 .1;.1 .1 1] has condition number 1.33, but geometric normalization makes it extremely ill-conditioned. Truncating its canonical SVD can discard an independent direction after a channel-unit change. Conversely, a triangular matrix with unit diagonal can have an enormous inverse norm while its RGA remains exactly the identity. Both are regression tests.

Numerical choices

• Keep an LU fast path, reusing its factorization for the inverse and infinity-norm condition estimate. κ∞*sqrt(eps(T)) < 1 is a conservative screen for using that result, not a singularity threshold.
• Before truncating a square BLAS input, try binary equilibration based on a maximum-weight matching of element exponents, using the existing Hungarian dependency. Bounded integer difference-constraint relaxation produces scales with matched maximum components in [1,2) and all components below 2, without forming an intermediate row-scaled matrix. This follows the idea of matching-based scaling; it does not add MC64 or another dependency.
• Keep the balanced inverse if its condition screen passes or its LU pivots exceed the factor-error scale γn*norm(L,Inf)*norm(U,Inf), with γn=n*eps(T)/(1-n*eps(T)). This is motivated by the LU backward-error bound, not a forward-error or exact-rank certificate. An exact complex rank-four regression has a spurious 5.65e-6 pivot above n*eps*maximum(abs,U) but below the factor scale. The independent triangular diagonal is retained. Single row/column max-norm balancing in either order was considered but still missed six randomized unit-change cases.
• Otherwise use Uhlmann's geometric normalization and the usual pinv cutoff, min(m,n)*eps(T)*σmax. Raw conditioning caused by channel units alone must not discard independent directions. Neither the LU nor SVD threshold is a certificate of mathematical singularity.
• If the normalized square matrix has full numerical rank, invert a copy equilibrated only by powers of two. These scales avoid extra rounding while the balanced values remain representable without underflow. In particular, the tests retain the large, correct RGA for [1 1;1 1+δ], with δ=2^-12 in Float32 and δ=2^-40 in Float64. High conditioning alone does not select a truncated result.
• For the UC result, if S=Dl*A*Dr, then Aᵁ=Dr*pinv(S)*Dl and the factors cancel: RGA(A)=S .* transpose(pinv(S)). This avoids recovering inverse scales. The original entries are scaled using split binary exponents and bounded factors; sign(A).*exp(L) can erase small differences, and separately forming exp(u)/exp(v) can exceed the exponent range.
• Keep real work arrays at the input floating precision. Use dimensionless log corrections, a precision-dependent default tol, a roundoff floor relative to the residual logs, and an explicit maxiter=10000 cap. Nonconvergence, nonfinite input and unrepresentable normalized support produce errors.

The existing tol keyword still controls normalization, not numerical rank. The matrix method gains optional maxiter. Singular/rectangular Float32 output now remains Float32. Nonsingular non-BLAS square types retain their existing LU path, including exact rational inverses; generalized BigFloat cases still require an appropriate generic SVD implementation. No dependency is added.

The UC construction preserves diagonal-unit changes; a Moore–Penrose inverse applied directly to the original matrix generally does not. References: Uhlmann's RGA paper and Appendix A, general UC construction, Theorem V.4, and Julia's pinv cutoff.

Validation

Windows x86_64, AMD Ryzen 7 7840H, one OpenBLAS ILP64 thread:

• The unmodified master method fails all three standalone regressions: uniform scaling, the finite but incorrect extreme-range inverse, and Float32 output precision.
• Julia 1.10.12: 216 new assertions cover deterministic failures, independent diagonal/scalar scaling, square/tall/wide/rank-deficient matrices, disconnected and zero support, complex phases and transpose semantics, exponent limits, resolvable weak directions, exact rational inverses, BigFloat's existing support boundary, type inference, explicit convergence failures, and the geometric-conditioning/triangular/complex-rank examples above.
• Fixed-seed independent experiments: 2,231 assertions passed, including 256-bit UC references constructed by closed-form log centering and known rank factorizations, rather than the production iteration/SVD. Random invertible inputs are checked against 256-bit inverses of the actual rounded input values. Seeds 10812026, 11081, 21081; these include 800 checks of well-conditioned matrices with tiny interactions and two-sided units, and 400 complex rank/phase checks.
• Final complete ControlSystemsBase suite, Julia 1.13.2: 7,356 passed, 6 existing broken, zero failures; 5 Aqua checks passed.
• Complete ControlSystems suite, Julia 1.10.12: 1,156 passed.
• git diff --check passes. Existing test tolerances are unchanged.

The new regression file is included by test_analysis.jl. It can also run separately with using ControlSystemsBase, LinearAlgebra, Random, Test; import SparseArrays: sparse; include("lib/ControlSystemsBase/test/test_relative_gain_array.jl").

Performance

BenchmarkTools 1.8.0, Julia 1.10.12, warmed functions, 1,000 samples, 10 evaluations/sample for ordinary matrices, one BLAS thread, no other Julia jobs running. Inputs are randn(rng,T,n,n)+n*I, seed 1081. Median microseconds for Float64:

Size	Master	This change	Always UC/SVD candidate
2×2	0.610	0.480	3.620
4×4	0.650	0.570	4.710
8×8	1.625	1.040	12.310
16×16	4.235	3.450	40.805
64×64	64.465	54.270	525.030

Float32 sizes 2/4/8/16/64 likewise show no measured regression; allocated bytes/counts are identical to master for every ordinary case. Always taking the geometric UC path was considered, but adds iteration and SVD to every frequency slice and has the conditioning problem above. The rank-two 3×3 counterexample now costs 5.0 → 7.4 µs (10,464 → 8,480 bytes): the extra square-inverse safeguard has a real cost. A random 4×6 input improves from 9.8 to 5.4 µs (16,520 → 7,328 bytes), each measured with one evaluation/sample. These are local microbenchmarks, not universal speedup guarantees.

An ordinary square call remains O(n³), and storage is O(mn). The bounded generalized path costs O(Kmn + mn·min(m,n)); square fallback additionally pays for matching and O(n³) relaxation/factorization. A conservative bound for the existing matrix-form Munkres implementation is O(n⁴). Ordinary calls bypass matching and iteration.

Numerical limits

This supplies a numerical rank model, not a certificate of mathematical singularity. In 50 fixed-seed 4×4 cases per requested condition number, median relative errors against a 256-bit inverse of the rounded input were:

Requested κ	Float32	Float64
10⁴	3.25e-5	5.57e-14
10⁸	approximately 1 in most cases	3.43e-10
10¹²	approximately 1 in most cases	4.50e-6
10¹⁶	approximately 1	approximately 1 in most cases

Generic Float32 inputs cannot reliably encode those smallest directions at the larger requested condition numbers. At Float64 10¹⁶, 40/50 cases discarded the small direction; others retained a sensitive inverse. A unit change can still cross the rank cutoff in this ambiguous region: maximum normalized scaling jumps reached 1.38 (Float32) and 2.50e14 (Float64), with a possibly projected result in the denominator. A projected UC result is different from the classical inverse of a strictly invertible rounded input. Resolved pivots also do not guarantee an accurate sensitive inverse. These failures of universal invariance/accuracy are documented rather than counted as passing cases.

Sparse support can converge slowly. A 32×33 bidiagonal-support matrix with values 2 and 3 reaches the default cap in Float64 and succeeds with maxiter=50000; both behaviors are tested. A 64×65 chain can exceed even that larger cap. The path's small projection spectral gap also means local eps-sized corrections do not imply eps-sized global error: measured 32-row UC errors were 5.20e-5 (Float32) and 1.02e-13 (Float64). Normalized entries themselves must fit in the working exponent range. This PR does not claim universal scaling invariance at rank/range boundaries or add generic high-precision SVD support.

Remaining generalized calculations can be sensitive after canonical scaling even when the raw input looks well-conditioned; they still apply the stated numerical rank model. Non-BLAS square types retain their existing LU behavior, without the new BLAS rank safeguards.

AI-assisted implementation, numerical experiments, and review with Codex. Independent references here refer to separate mathematical implementations, not a human review.

[codecov]

Codecov Report

❌ Patch coverage is 98.06452% with 3 lines in your changes missing coverage. Please review.
✅ Project coverage is 92.39%. Comparing base (9baa912) to head (e7a0b92).
⚠️ Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
lib/ControlSystemsBase/src/analysis.jl	98.06%	3 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1092      +/-   ##
==========================================
+ Coverage   92.26%   92.39%   +0.12%     
==========================================
  Files          42       42              
  Lines        6170     6299     +129     
==========================================
+ Hits         5693     5820     +127     
- Misses        477      479       +2     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.