diff --git a/meegkit/utils/covariances.py b/meegkit/utils/covariances.py
index bf31093a..2f05bae5 100644
--- a/meegkit/utils/covariances.py
+++ b/meegkit/utils/covariances.py
@@ -510,5 +510,15 @@ def ehess(X, U):
     problem = Problem(manifold=manifold, cost=cost, euclidean_gradient=egrad,
                       euclidean_hessian=ehess)
     Xsol = solver.run(problem, initial_point=U0)
-
-    return S0, Xsol.point
+    X = Xsol.point
+
+    # The Grassmann manifold only fixes the *subspace* spanned by X, so X carries
+    # an arbitrary internal rotation.  Diagonalise the k×k projected matrix to
+    # obtain proper Ritz pairs: after rotation, X.T @ L @ X is diagonal and each
+    # column of X is a genuine approximate eigenvector of L.
+    T = X.T @ L @ X
+    T = (T + T.T) / 2  # symmetrise against numerical noise
+    S, R = linalg.eigh(T)
+    X = X @ R
+
+    return np.real(S), np.real(X)
diff --git a/tests/test_cov.py b/tests/test_cov.py
index 0ca46108..a7a25ccc 100644
--- a/tests/test_cov.py
+++ b/tests/test_cov.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pytest
 from numpy.testing import assert_almost_equal
 
 from meegkit.utils import convmtx, tscov, tsxcov
@@ -99,6 +100,32 @@ def test_convmtx():
                   ])
     )
 
+
+def test_nonlinear_eigenspace_consistent_eigpairs():
+    """Each returned eigenvalue must match the Rayleigh quotient of its column."""
+    pytest.importorskip("pymanopt")
+    from meegkit.utils.covariances import nonlinear_eigenspace
+
+    A = rng.standard_normal((6, 6))
+    L = A.T @ A + np.eye(6) * 1e-6
+
+    S, X = nonlinear_eigenspace(L, 6)
+    # X.T @ L @ X should be (near) diagonal; check the diagonal matches S
+    # element-wise so a permutation between S and X columns cannot slip through.
+    rayleigh = np.real(np.diag(X.T @ L @ X))
+
+    np.testing.assert_allclose(
+        S,
+        rayleigh,
+        rtol=1e-3,
+        atol=1e-6,
+    )
+
+    # Off-diagonal entries of the projected matrix should be negligible.
+    proj = X.T @ L @ X
+    off_diag = proj - np.diag(np.diag(proj))
+    assert np.max(np.abs(off_diag)) < 1e-6 * np.max(np.abs(np.diag(proj)))
+
 if __name__ == "__main__":
     import pytest
     pytest.main([__file__])