File tree Expand file tree Collapse file tree
Expand file tree Collapse file tree Original file line number Diff line number Diff line change @@ -505,8 +505,13 @@ def ehess(X, U):
505505 Xsol = solver .run (problem , initial_point = U0 )
506506 X = Xsol .point
507507
508- # Keep eigpairs consistent: compute each optimized vector's Rayleigh
509- # quotient against L instead of returning the random-init eigenvalues S0.
510- S = np .real (np .diag (X .T @ L @ X ))
511-
512- return S , X
508+ # The Grassmann manifold only fixes the *subspace* spanned by X, so X carries
509+ # an arbitrary internal rotation. Diagonalise the k×k projected matrix to
510+ # obtain proper Ritz pairs: after rotation, X.T @ L @ X is diagonal and each
511+ # column of X is a genuine approximate eigenvector of L.
512+ T = X .T @ L @ X
513+ T = (T + T .T ) / 2 # symmetrise against numerical noise
514+ S , R = linalg .eigh (T )
515+ X = X @ R
516+
517+ return np .real (S ), np .real (X )
Original file line number Diff line number Diff line change @@ -91,23 +91,30 @@ def test_convmtx():
9191
9292
9393def test_nonlinear_eigenspace_consistent_eigpairs ():
94- """Returned eigenvalues must match returned eigenvectors ."""
94+ """Each returned eigenvalue must match the Rayleigh quotient of its column ."""
9595 pytest .importorskip ("pymanopt" )
9696 from meegkit .utils .covariances import nonlinear_eigenspace
9797
9898 A = rng .standard_normal ((6 , 6 ))
9999 L = A .T @ A + np .eye (6 ) * 1e-6
100100
101101 S , X = nonlinear_eigenspace (L , 6 )
102+ # X.T @ L @ X should be (near) diagonal; check the diagonal matches S
103+ # element-wise so a permutation between S and X columns cannot slip through.
102104 rayleigh = np .real (np .diag (X .T @ L @ X ))
103105
104106 np .testing .assert_allclose (
105- np . sort ( np . real ( S )) ,
106- np . sort ( rayleigh ) ,
107+ S ,
108+ rayleigh ,
107109 rtol = 1e-3 ,
108110 atol = 1e-6 ,
109111 )
110112
113+ # Off-diagonal entries of the projected matrix should be negligible.
114+ proj = X .T @ L @ X
115+ off_diag = proj - np .diag (np .diag (proj ))
116+ assert np .max (np .abs (off_diag )) < 1e-6 * np .max (np .abs (np .diag (proj )))
117+
111118if __name__ == "__main__" :
112119 import pytest
113120 pytest .main ([__file__ ])
You can’t perform that action at this time.
0 commit comments