@@ -200,7 +200,7 @@ def lyap(A, Q, C=None, E=None, method=None, symmetric_kwargs=None):
     # Solve the generalized Lyapunov equation
     elif C is None and E is not None:
         # Check to make sure input matrices are the right shape and type
-        _check_shape(Q, n, n, square=True, symmetric=True, name="Q")
+        _check_shape(Q, n, n, square=True, symmetric=True, name="Q", symmetric_kwargs=symmetric_kwargs)
         _check_shape(E, n, n, square=True, name="E")
 
         if method == 'scipy':
@@ -412,7 +412,7 @@ def dlyap(A, Q, C=None, E=None, method=None, symmetric_kwargs=None):
     # Solve the generalized Lyapunov equation
     elif C is None and E is not None:
         # Check to make sure input matrices are the right shape and type
-        _check_shape(Q, n, n, square=True, symmetric=True, name="Q")
+        _check_shape(Q, n, n, square=True, symmetric=True, name="Q", symmetric_kwargs=symmetric_kwargs)
         _check_shape(E, n, n, square=True, name="E")
 
         if method == 'scipy':
@@ -803,4 +803,4 @@ def _is_symmetric(M, symmetric_kwargs=None):
 
     if np.iscomplexobj(M):
         return sp.linalg.ishermitian(M, **symmetric_kwargs)
-    return sp.linalg.issymmetric(M, **symmetric_kwargs)
+    return sp.linalg.issymmetric(M, **symmetric_kwargs)
@@ -318,7 +318,7 @@ def create_estimator_iosystem(
         control_indices=None, disturbance_indices=None,
         estimate_labels='xhat[{i}]', covariance_labels='P[{i},{j}]',
         measurement_labels=None, control_labels=None,
-        inputs=None, outputs=None, states=None, **kwargs):
+        inputs=None, outputs=None, states=None, symmetric_kwargs=None, **kwargs):
     r"""Create an I/O system implementing a linear quadratic estimator.
 
     This function creates an input/output system that implements a
@@ -416,6 +416,9 @@ def create_estimator_iosystem(
     name : string, optional
         System name (used for specifying signals). If unspecified, a generic
         name 'sys[id]' is generated with a unique integer id.
+    symmetric_kwargs : dict, optional
+        Keyword arguments passed to `scipy.linalg.issymmetric` or
+        `scipy.linalg.ishermitian`.
 
     Notes
     -----
@@ -498,7 +501,7 @@ def create_estimator_iosystem(
     if P0 is None:
         # Initialize P0 to the steady state value
         _, P0, _ = lqe(A, G, C, QN, RN)
-    P0 = _check_shape(P0, sys.nstates, sys.nstates, symmetric=True, name='P0')
+    P0 = _check_shape(P0, sys.nstates, sys.nstates, symmetric=True, name='P0', symmetric_kwargs=symmetric_kwargs)
 
     # Figure out the labels to use
     estimate_labels = _process_labels(