Skip to content

Improve numerical symmetry check #1248

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
CNZHM666 wants to merge 4 commits into
base: main
Choose a base branch
from
+27 −7
Open

Improve numerical symmetry check #1248

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
29fd88c
Improve numerical symmetry check
CNZHM666 Sep 13, 2026
fd72ea8
Fix whitespace in symmetry tests
CNZHM666 Sep 15, 2026
0dfa703
Fix remaining style issues
CNZHM666 Sep 15, 2026
7ac6daa
Use SciPy symmetry check method
CNZHM666 Sep 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 5 additions & 6 deletions control/mateqn.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@

import numpy as np
import scipy as sp
from numpy import eye, finfo, inexact
from numpy import eye, finfo
from scipy.linalg import eigvals, solve

from .exception import ControlArgument, ControlDimension, ControlSlycot, \
Expand Down Expand Up @@ -787,8 +787,7 @@ def _check_shape(M, n, m, square=False, symmetric=False, name="??"):
# Utility function to check if a matrix is symmetric
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()
return (
sp.linalg.norm(M - M.conj().T, 1)
<= np.spacing(sp.linalg.norm(M, 1)) * 100
)
Comment on lines +790 to +793

@slivingston slivingston Sep 15, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Where is the "SciPy symmetry check method" in this commit?

A few quick questions:

  1. Can you read the comments that I linked to previously, which describe the desired changes, without using an AI bot?
  2. Can you try to write the change without using an AI bot?

Using AI tools is OK, but I want to make sure that you understand the proposed solution (and thus, can write/understand the code).

@CNZHM666 CNZHM666 Sep 16, 2026

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. I read the linked comments and SciPy documentation myself.

I misunderstood the requested change earlier. I now understand that I should use SciPy’s built-in symmetry/Hermitian check rather than reimplementing the norm/spacing test. I will rewrite the change myself.

Should the implementation use SciPy’s default exact comparison, or should atol/rtol also be passed through?

@slivingston slivingston Sep 16, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@CNZHM666 I recommend that atol/rtol are passed through. Let me highlight text from one of the comments that I linked:

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.

@CNZHM666 CNZHM666 Sep 16, 2026

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I traced the call chain. lqe/dlqe already accept keyword arguments and call care/dare, while care, dare, lyap, and dlyap currently do not expose symmetry tolerances. Do you want atol/rtol to be added to all of these public matrix-equation functions, or only passed through the LQE path for this issue?

@slivingston slivingston Sep 16, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's try to use kwargs or a similar name to support passing through the parameter to scipy. (Probably better not to add explict atol, rtol parameters.) As in the above comment, "There are several other examples where we pass down options to scipy functions." Find those examples, and try to follow the pattern in this PR.

23 changes: 22 additions & 1 deletion control/tests/mateqn_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@
import pytest
from scipy.linalg import eigvals, solve

from control.mateqn import lyap, dlyap, care, dare
from control.mateqn import lyap, dlyap, care, dare, _is_symmetric
from control.exception import ControlArgument, ControlDimension


Expand Down Expand Up @@ -466,3 +466,24 @@ def test_raise(self):
cdare(A, B, Qfs, R, S, E)
with pytest.raises(ControlArgument):
cdare(A, B, Q, Rfs, S, E)

def test_is_symmetric_scale_aware(self):

@slivingston slivingston Sep 15, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you improve the whitespace of your changes?

  • Always at least 1 blank line before the start of a function def.
  • The blank line between each M definition and respective assert statement below does not improve clarity.

@CNZHM666 CNZHM666 Sep 15, 2026

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for pointing that out. I have cleaned up the whitespace and rerun the relevant tests.

M = np.array([
[1e8, 1e8],
[1e8 + 1e-8, 1e8]
])
assert _is_symmetric(M)

def test_is_symmetric_rejects_asymmetric(self):
M = np.array([
[1., 2.],
[5., 1.]
])
assert not _is_symmetric(M)

def test_is_symmetric_complex_hermitian(self):
M = np.array([
[1., 2. + 1.j],
[2. - 1.j, 3.]
])
assert _is_symmetric(M)