xLAEIN: apply the growth test to the growth, not to |x| alone - #1400
Open
ACSimon33 wants to merge 2 commits into
Open
xLAEIN: apply the growth test to the growth, not to |x| alone#1400ACSimon33 wants to merge 2 commits into
ACSimon33 wants to merge 2 commits into
Conversation
The residual of the vector a solve returns is SCALE*|v|/|x|, so the acceptance test belongs on that ratio. Testing |x| against GROWTO alone made it depend on the starting vector: GROWTO suits the restart vectors, whose 1-norm is of order SQRT(N)*EPS3, and was SQRT(N) too lenient for the initial vector, whose 1-norm is N*EPS3 - and that is where almost every vector is accepted. One NEP matrix reached 27.2 against a test threshold of 20; the worst ratio is now 4.2. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
❌ 156 Tests Failed:
View the top 3 failed test(s) by shortest run time
View the full list of 76 ❄️ flaky test(s)
To view more test analytics, go to the Test Analytics Dashboard |
Each starting vector gets exactly one solve, so there is no iterate beyond the first; 0 and 1 say that where i and i+1 do not. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
martin-frbg
approved these changes
Sep 10, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
xLAEINdecides when one step of inverse iteration has produced an acceptable eigenvector. The residual of that step is governed by how much the solution grew relative to the starting vector, but the test compared the solution norm againstGROWTOon its own. That silently assumes every starting vector has the same constant norm, which is not the case.The criterion
One step of inverse iteration with shift
sigmaisThe residual of the vector it produces is an identity, not a bound:
The reciprocal of the right-hand side is the norm growth, and the vector is accepted once that growth is large enough that (1) delivers a residual at the level of the accuracy of the eigenvalue:
Then
|| r || <= (1/C) n eps ||A|| || v^(i+1) ||, the standard requirement on an approximate eigenpair. The essential point is that the criterion is a statement about the ratio: the norm of the starting vector belongs in it.This is Dhillon's section 3, issue IV (PDF), which attributes the criterion to Wilkinson's The Algebraic Eigenvalue Problem, p. 324.
What the old code did
Here
VNORMis the 1-norm of the solution andSCALEisxLATRS's scale factor, so dividing the test through bySCALEtimes the starting vector's 1-norm puts it in the same form as (2):The right-hand side is what (2) asks for only if
So the starting vector's norm was not so much ignored as assumed constant, and assumed equal to
ROOTN*EPS3in particular. That is close to the norm of the restarting vectors the routine builds after a failed attempt, so the assumption holds to within a factor of two for those. It does not hold for the vector the iteration actually starts from:NOINIT)EPS3N*EPS3ROOTNEPS3, thenEPS3/(ROOTN+1), one entry- EPS3*ROOTN~2*ROOTN*EPS3~2INITV='U')ROOTN*EPS3ROOTN*EPS3 .. N*EPS31 .. ROOTNThe consequence follows from (1). Accepting a growth
ROOTNtimes smaller than required admits a residualROOTNtimes larger, so the bound degrades fromO( n eps ||H|| )toO( n**1.5 eps ||H|| ). NEP test 11 measures exactly that quantity divided bynagainst a threshold of 20, so the old bound permits a ratio of order10*ROOTNand exceeds the threshold as soon asN > 4. The initial vector is where 28,451 of the 28,452 eigenvectors in the NEP test suite are accepted, so it is the case that usually decides the outcome, and one matrix crossed the threshold on aarch64 gfortran:Instrumenting the routine to print the residual of every accepted vector shows that case as
n=10, try 1, residual 27.9in units ofN*EPS3*|x|, while the very next starting vector gives0.0209- a residual 1300 times smaller for the same eigenvalue. So this was never an ill-conditioned eigenvector that no starting vector could resolve. A far better one was one restart away, and the test stopped short because it wasn't strict enough.The fix
xLAEINtakes exactly one solve from each starting vector - on insufficient growth it chooses a new starting vector rather than iterating on the result - soiis always 0 here. Mapping the notation of (2) onto this routine:A - sigma IBEPS3tau^(0)SCALExLATRSto avoid overflowv^(0)VRon entry, the starting vectorV0NORMv^(1)VRafter the solveV1NORMepstimes the norm of AEPS3xHSEINtoULPtimes the block normnNCTENTHso (2) reads
GROWTObecomes the right-hand side of (3), and the test is then (3) verbatim:GROWTOchanges value, from0.1/SQRT(N)to0.1/(N*EPS3), because it changes meaning: it is no longer a bound on||x||but the growth ratio that (3) requires. Nothing else in the routine uses it.The effective threshold each starting vector now faces, relative to before:
TENTH*SCALE/ROOTNROOTNtimes stricterTENTH*SCALE/ROOTN~2times stricterTENTH*SCALE/ROOTN1toROOTNtimes stricterEffect
Worst
xneptest-11 ratio,|HX - XW| / (|H| |X| ulp)on aarch64 now:Because only the accept/reject decision changes and no iterate is touched, vectors that already passed come back bit-identical - visible above in the D and Z columns.
Notes
The
0.1/SQRT(N)threshold is inherited verbatim from EISPACK'sinvit, which hasgrowto = 0.1d0 / ukroot,ukroot = dsqrt(uk)under the commentgrowto is the criterion for the growth, and which implements [3]. The same missing normalisation is therefore present there. Every commit toSRC/[scdz]laein.fon master is cosmetic, so the criterion is untouched since that import.xSTEINavoids the problem from the other direction:SRC/dstein.frescales the right-hand side before every solve so its norm is a known constant, which is what makes its absolute test onNRMsound.References