better symmetry check test

Hi,

While porting our MATLAB code to Python using the control library, we noticed that the current symmetry check:

def _is_symmetric(M):
    M = np.atleast_2d(M)
    if isinstance(M[0, 0], inexact):
        eps = finfo(M.dtype).eps
        return ((M - M.T) < eps).all()
    else:
        return (M == M.T).all()

has some limitations:

It performs elementwise comparison without using abs(), so small negative differences (e.g., from roundoff) can cause false failures.
It uses a fixed eps, which is stricter than MATLAB's behavior and doesn’t scale with matrix magnitude.

Would it be possible to adopt a more numerically robust check, similar to what SciPy uses?

for ind, mat in enumerate((q, r)):
    if norm(mat - mat.conj().T, 1) > np.spacing(norm(mat, 1)) * 100:
        raise ValueError(f"Matrix {'qr'[ind]} should be symmetric/hermitian.")

This approach allows for small, scale-aware numerical asymmetries and works well for both real symmetric and complex Hermitian matrices.

Prithvi

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.5 where 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 dlqe takes A, G, C, Q as inputs, and computes the process noise covariance internally as QN = G @ Q @ G.T, there's no way for the user to ensure that QN is symmetric to machine precision when G ≠ I.

While I could compute QN externally and pass G = I to dlqe, that defeats the purpose for general users who rely on G for modeling how noise enters the system.

Hence the request.

@dakeprithvi

@dakeprithvi Can you provide the example that is triggering the error? I'm assuming the issue is that Q is not (sufficiently) symmetric?

A couple of thoughts on things we might do:

We should almost certainly replace _issymmetric with scipy.issymmetric, since there is no reason for the duplication.
We could add a way to allow rtol and atol to be passed through to scipy.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 Q or QN (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?

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) / 2 in the code.

@ilayn

However, I’d recommend adopting the improved symmetry test (as @ilayn suggested), since forcing QN = (QN + QN.T)/2 will 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.

@SepehrMoalemi

I definitely think the matrix needs to be made symmetric inside dlqe. I believe I ran into this issue with dlqr with 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.

use symmetry check from SciPy

fixed in #1248

Navigation Menu

Repository navigation

better symmetry check test #1174

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

better symmetry check test #1174

Description

Activity

ilayn commented on Jul 23, 2025

dakeprithvi commented on Jul 23, 2025

murrayrm commented on Jul 23, 2025

dakeprithvi commented on Jul 23, 2025

murrayrm commented on Jul 23, 2025

dakeprithvi commented on Jul 23, 2025

dakeprithvi commented on Jul 23, 2025

sdahdah commented on Jul 23, 2025

SepehrMoalemi commented on Jul 23, 2025

ilayn commented on Jul 23, 2025

slivingston commented on Sep 29, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions