diff --git a/pygpcca/_sorted_schur.py b/pygpcca/_sorted_schur.py
index b45e833..2c6eb5f 100644
--- a/pygpcca/_sorted_schur.py
+++ b/pygpcca/_sorted_schur.py
@@ -106,7 +106,7 @@ def _check_conj_split(eigenvalues: ArrayLike) -> bool:
     """
     last_eigenvalue, second_last_eigenvalue = eigenvalues[-1], eigenvalues[-2]
     splits_block = False
-    if last_eigenvalue.imag > EPS:
+    if abs(last_eigenvalue.imag) > EPS:
         splits_block = not np.isclose(last_eigenvalue, np.conj(second_last_eigenvalue))
 
     return splits_block
@@ -221,7 +221,7 @@ def sorted_krylov_schur(
     eigenvalues_error
         Array of shape `(k,)` containing the error, based on the residual
         norm, of the `i`th eigenpair at index `i`.
-    """  # noqa: D205, D400
+    """
     # We like to thank A. Sikorski and M. Weber for pointing us to SLEPc for partial Schur decompositions of
     # sparse matrices.
     # Further parts of sorted_krylov_schur were developed based on the function krylov_schur
diff --git a/tests/conftest.py b/tests/conftest.py
index 45f8464..3204d3b 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -39,6 +39,25 @@ def assert_allclose(actual, desired, rtol=1.0e-5, atol=1.0e-8, err_msg="", verbo
     return assert_allclose_np(actual, desired, rtol=rtol, atol=atol, err_msg=err_msg, verbose=True)
 
 
+def normalize_conj_pairs(eigenvalues: np.ndarray) -> np.ndarray:
+    """Normalize conjugate pair ordering so positive imaginary part always comes first.
+
+    The ordering of eigenvalues within complex conjugate pairs from a Schur
+    decomposition is not guaranteed by LAPACK and may vary across scipy/numpy
+    versions. This function ensures a canonical ordering for testing purposes.
+    """
+    result = eigenvalues.copy()
+    i = 0
+    while i < len(result) - 1:
+        if abs(result[i].imag) > 1e-10 and np.isclose(result[i], np.conj(result[i + 1])):
+            if result[i].imag < 0:
+                result[i], result[i + 1] = result[i + 1], result[i]
+            i += 2
+        else:
+            i += 1
+    return result
+
+
 def get_known_input(Tc: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
     assert not np.allclose(Tc, 0.0), "Tc doesn't seem to be a count matrix."
     assert Tc.dtype == np.float64, "Expected double precision"
diff --git a/tests/test_gpcca.py b/tests/test_gpcca.py
index 20dc80b..36d0402 100644
--- a/tests/test_gpcca.py
+++ b/tests/test_gpcca.py
@@ -48,7 +48,13 @@
     gpcca_coarsegrain,
     _initialize_rot_matrix,
 )
-from tests.conftest import mu, assert_allclose, get_known_input, skip_if_no_petsc_slepc
+from tests.conftest import (
+    mu,
+    assert_allclose,
+    get_known_input,
+    normalize_conj_pairs,
+    skip_if_no_petsc_slepc,
+)
 from pygpcca._sort_real_schur import sort_real_schur
 
 eps = np.finfo(np.float64).eps * 1e10
@@ -845,8 +851,14 @@ def test_P_2_LM(
         assert_allclose(g.crispness_values, crispness_values_P_2_LM)
         assert_allclose(g.optimal_crispness, optimal_crispness_P_2_LM)
         assert_allclose(n_m, n_m_P_2_LM)
-        assert_allclose(g.top_eigenvalues, top_eigenvalues_P_2_LM)
-        assert_allclose(g.dominant_eigenvalues, top_eigenvalues_P_2_LM[:n_m])
+        assert_allclose(
+            normalize_conj_pairs(g.top_eigenvalues),
+            normalize_conj_pairs(top_eigenvalues_P_2_LM),
+        )
+        assert_allclose(
+            normalize_conj_pairs(g.dominant_eigenvalues),
+            normalize_conj_pairs(top_eigenvalues_P_2_LM[:n_m]),
+        )
 
     def test_split_warning_LM(self, P_2: np.ndarray):
         g = GPCCA(P_2, eta=None, z="LM")
@@ -930,8 +942,14 @@ def test_P_2_LR(
         assert_allclose(g.crispness_values, crispness_values_P_2_LR)
         assert_allclose(g.optimal_crispness, optimal_crispness_P_2_LR)
         assert_allclose(n_m, n_m_P_2_LR)
-        assert_allclose(g.top_eigenvalues, top_eigenvalues_P_2_LR)
-        assert_allclose(g.dominant_eigenvalues, top_eigenvalues_P_2_LR[:n_m])
+        assert_allclose(
+            normalize_conj_pairs(g.top_eigenvalues),
+            normalize_conj_pairs(top_eigenvalues_P_2_LR),
+        )
+        assert_allclose(
+            normalize_conj_pairs(g.dominant_eigenvalues),
+            normalize_conj_pairs(top_eigenvalues_P_2_LR[:n_m]),
+        )
 
     def test_split_warning_LR(self, P_2: np.ndarray):
         g = GPCCA(P_2, eta=None, z="LR")