Repository navigation
better symmetry check test #1174
Description
Activity
Hi I happen to be the person who wrote that piece of code in SciPy.
That symmetricity test is also not a good one and should be replaced. You can use issymmetric and ishermitian from SciPy yourself. I wrote them in C and they are considerably faster, moreover they have quick exit implemented; they quit at the first instance of mismatch and don't go all the way to the end regardless (unlike
arr.all()).If symmetricity is required just do
(M + M.T)*0.5where you need it. It will rarely have considerable performance impact compared to other operations in control codebases.Thanks for the update.
Sorry I should have also added the context where this issue arises in the control library:Since
dlqetakesA, G, C, Qas inputs, and computes the process noise covariance internally asQN = G @ Q @ G.T, there's no way for the user to ensure thatQNis symmetric to machine precision whenG ≠ I.While I could compute
QNexternally and passG = Ito dlqe, that defeats the purpose for general users who rely on G for modeling how noise enters the system.Hence the request.
@dakeprithvi Can you provide the example that is triggering the error? I'm assuming the issue is that
Qis not (sufficiently) symmetric?A couple of thoughts on things we might do:
- We should almost certainly replace
_issymmetricwithscipy.issymmetric, since there is no reason for the duplication. - We could add a way to allow
rtolandatolto be passed through toscipy.issymmetric, so that the user can control the behavior better. There are several other examples where we pass down options to scipy functions. - We might also include an option to symmetrize either
QorQN(via(M + M.T)*0.5, as @ilayn suggests).
Whoever picks up this issue should look through the code and see what makes the most sense.
As a quick fix in the meantime, you can presumably use the symmetrization fix as a workaround?
- We should almost certainly replace
Thanks for the quick reply!
Here's a test case (for me the code runs 5/10 times with success)import numpy as np import control # Dimensions n1, n2, p, g = 3, 2, 3, 3 n = n1 + n2 # Random system setup rng = np.random.default_rng() A = np.block([ [rng.random((n1, n1)), rng.random((n1, n2))], [np.zeros((n2, n1)), 2*np.eye(n2) + rng.random((n2, n2))] ]) G = np.vstack([rng.random((n1, g)), np.zeros((n2, g))]) C = np.hstack([rng.random((p, n1)), rng.random((p, n2))]) # Noise covariances Q = rng.standard_normal((g, g)); Q = Q @ Q.T R = rng.standard_normal((p, p)); R = R @ R.T Q = (Q + Q.T) / 2 # Doesn't help, what matters is QN = G @ Q @ G.T R = (R + R.T) / 2 # Kalman gain L, P, E = control.dlqe(A, G, C, Q, R)Very helpful! Sounds like we need to put in a
QN = (QN + QN.T) / 2in the code.However, I’d recommend adopting the improved symmetry test (as @ilayn suggested), since forcing
QN = (QN + QN.T)/2will silently make even incorrect user-supplied Q appear valid.Leaving the issue open for now as the change hasn’t been implemented yet. Devs can close it once done.
Thanks.I definitely think the matrix needs to be made symmetric inside
dlqe. I believe I ran into this issue withdlqrwith a colleague (and we got distracted before we could raise an issue 🙃)@SepehrMoalemi do you still have the code that reproduces this?
I ran into this issue a year ago using
lqe. I ended up just doing the following:# Compute Estimator using LQE L, _, _ = control.lqe(A, B, C, Qe, Re) # Note, QN can be illconditioned and should be symmetric QN = B @ Qe @ B.T QN = (QN + QN.T)/2 _, _, LT = control.care(A, C.T, QN, Re, B_s="C", Q_s="QN", R_s="RN", S_s="NN")
I think the current implementation shortcomings are agreed here. So we should do something about that but even if that is fixed, your test code still might generate bad results.
Note that, the symmetricity and other issues always need to start with the limitations of the algorithm such as observability and the weight matrix semi-definiteness. So when you start with large random matrices they tend to be very close to singularity (you can actually show this theoretically, Wigner, Wishart, real Ginibre ensembles and whatnot). This is not always numerically guaranteed to have sane answers especially in Riccati based algorithms where you solve large matrix equations.
Hence I would suggest you do tests with well-conditioned problems to start with. These are not numerically robust algorithms hence you might encounter strange results even after this symmetricity check problem is fixed.
Reacted by Steven Dahdah- added a commit that references this issue
on Sep 29, 2026 fixed in #1248
Hi,
While porting our MATLAB code to Python using the control library, we noticed that the current symmetry check:
has some limitations:
Would it be possible to adopt a more numerically robust check, similar to what SciPy uses?
This approach allows for small, scale-aware numerical asymmetries and works well for both real symmetric and complex Hermitian matrices.
Prithvi