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:

1. It performs elementwise comparison without using abs(), so small negative differences (e.g., from roundoff) can cause false failures.
2. 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