REVIEW 3 major objections 4 minor 30 references
Rayleigh Quotient Iteration, cubic convergence, and second covariant derivative
T0 review · 3 major / 4 minor · reviewed 2026-08-14 · deepseek-v4-flash
Pith's one-line read A single theorem now governs quadratic and cubic convergence for Rayleigh quotient iterations in eigenproblems, constrained optimization, and tensor eigenpairs.
desk verdict A genuinely unifying framework for Rayleigh quotient iteration with a real proof gap in the main theorem; worth refereeing, not desk-rejecting. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The central object is the vector Lagrangian $L(x,\lambda)$ together with its constraint $C(x)$, and the generalized Rayleigh quotient $R(x)$ that selects $\lambda$ at each iterate. The machinery that carries the argument is the projected generalized Hessian $\Pi_{L_\lambda} L_x$, obtained from a left inverse $L_\lambda^{-}$ of $L_\lambda$, restricted to the tangent space $JC(x)\eta=0$; uniform invertibility of this map is inequality (6.7). For the cubic statements the load-bearing object is the tensor $G(x)[\eta^{[2]}] = -\tfrac{1}{2} L_{xx}(x)[\eta^{[2]}] - L_{x\lambda}(x)[\eta, JR[\eta]] - \tfrac{1}{2} L_{\lambda\lambda}(x)[(JR[\eta])^{[2]}]$, which acts as the second covariant derivative of the Lagrangian along the iteration direction. The Schur form — solving $\zeta=-L_x^{-1}L_\lambda$ and $\xi=L_x^{-1}L$, then projecting with $(JC\,\zeta)^{-1}JC\,\xi$ — is the computational realization that links the abstract theorem to classical resolvent equations and to tensor and eigenvector algorithms.
What would settle it
Take a smooth Lagrangian and a curved constraint satisfying the stated hypotheses, for example $C(x)=x_1^2+x_2^2-1$ with $L(x,\lambda)=F(x)-x\lambda$ and a Rayleigh quotient $R$ with nonzero differential at the solution. Run iteration (6.8) in high-precision arithmetic from a nearby initial point and compute $Q_i=\log\|x_{i+1}-v\|/\log\|x_i-v\|$. If $Q_i$ does not approach $2$ for some such problem while all stated hypotheses hold, the missing normal-component control breaks the claimed quadratic rate; if $Q_i$ always approaches $2$, the control is automatic in these cases.
Extended reading notes
Core claim
On the paper's own terms, the discovery is that one construction governs all of these iterations. Let $L(x,\lambda)$ be a vector Lagrangian with constraint $C(x)=0$, and let $R$ be any $C^1$ Rayleigh quotient with $R(v)=\mu$ at a solution. Choose a left inverse $L_\lambda^{-}$ of $L_\lambda$, set $\Pi_{L_\lambda}=I - L_\lambda L_\lambda^{-}$, and solve the projected Hessian equation $\Pi_{L_\lambda} L_x(x,R(x))\,\eta = -\Pi_{L_\lambda} L(x,R(x))$ on the tangent space $JC(x)\eta=0$, updating by a retraction $x_{i+1}=r(x_i,\eta)$. This iteration converges quadratically to $(v,\mu)$ whenever the projected Hessian is uniformly invertible on the tangent space. If $L$ is $C^3$ and the retraction is second-order, the Rayleigh–Chebyshev step that appends the correction $T$ defined by $\Pi_{L_\lambda} L_x T = \Pi_{L_\lambda} G$ converges cubically, and when the projection of $G(v)$ vanishes the unmodified RQI is already cubic. The paper also gives the Schur-form realization of these equations when $L_x$ is invertible, which in the eigenvector case reproduces the classical resolvent formula and in constrained optimization reproduces the Riemannian Newton update equations.
Load-bearing premise
The load-bearing premise is uniform invertibility of the projected Hessian on the tangent space, inequality (6.7), together with the unstated premise that the displacement from the current point to the solution has a normal component controlled at second order, so that the tangent-space lower bound can honestly be applied to it.
Editorial extensions
If this is right
- For the eigenvector problem with the standard Rayleigh quotient, the generalized RQI is exactly classical RQI, and the condition $\Pi_{L_\lambda}(v)G(v)=0$ reproduces the known cubic convergence for normal matrices.
- For constrained optimization with $H=J_C^T$, the projected-Hessian equation is Riemannian Newton on the embedded manifold, giving a unified quadratic-convergence proof for that method and for feasible projected SQP.
- For real tensor eigenpairs, the Schur-form RQI is equivalent to the existing Newton-correction iteration on the sphere but avoids forming the projected Hessian; the paper reports about 16% faster runs, and about 34% faster with one redundant tensor evaluation removed.
- For nonnormal matrices, the two-sided RQI is recovered by writing left and right eigenvectors as a constrained Lagrangian, and its cubic convergence follows from the same $G$ criterion.
- A unitary version of the tensor RQI computes most complex eigenpairs of moderate tensors in minutes, and identifies the real pairs as a by-product without a homotopy run.
Reading between the lines
- The theorem's statement can be read as implicitly requiring a second-order normal estimate for the displacement $x_{i+1}-v$, since the proof applies the tangent-space lower bound (6.7) to it; if that estimate is not automatic on curved constraints, adding it as an explicit hypothesis is the cleanest repair.
- The freedom to choose any consistent $R$ suggests a design principle: choose the Rayleigh quotient so that the projected $G(v)$ vanishes, obtaining cubic convergence from plain RQI rather than paying for the Chebyshev correction; testing different left inverses on tensor eigenpairs would show whether such choices are practical.
- The paper's own numerical remark that a detailed global convergence analysis for the all-complex-eigenpair search is still needed, and that the last ten percent of pairs dominate runtime, points to the natural next problem of proving basin-of-attraction results or using RQI as a local accelerator inside a global polynomial solve.
- Because the Schur form becomes unstable when $L_x$ is nearly singular, the tangent form is the natural fallback near convergence; a hybrid that switches forms based on the conditioning of $L_x$ would extend the framework to harder nonnormal problems.
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper introduces a general framework for Rayleigh quotient iteration (RQI) applied to systems L(x, λ) = 0, C(x) = 0, where L is called a vector Lagrangian and C is a constraint. The authors define a generalized Rayleigh quotient R(x) and a generalized projected Hessian Π Lλ Lx, and state a convergence theorem (Theorem 6.2) asserting quadratic convergence of a generalized RQI, cubic convergence of a Rayleigh-Chebyshev variant, and cubic convergence of the plain RQI under a zero condition on a tensor G. They also provide a Schur-form realization of the iteration, apply the framework to tensor eigenpairs, constrained optimization, nonlinear eigenvalue problems, two-sided RQI, and Grassmann/Stiefel manifolds, and report numerical experiments including an algorithm for enumerating complex tensor eigenpairs.
Significance. If the convergence theorem is correct, the paper would provide a unified theory covering classical RQI, Riemannian Newton, and several recent tensor eigenvalue algorithms, while also producing a new Rayleigh-Chebyshev iteration. The numerical results, especially the complex tensor eigenpair enumeration, are potentially valuable and the framework is elegant. The paper ships with open-source code and the derivations are self-contained. However, the central proof has gaps that affect the claimed rates, and the complex enumeration claim relies on an unproven heuristic about random starting points. The framework is significant but the theoretical core needs repair before the claims can be accepted.
major comments (3)
- [§6, proof of Theorem 6.2, after Eq. (6.13)] The proof applies inequality (6.7) to the ambient difference x_{i+1} - v, but (6.7) is stated only for tangent vectors ψ in Null(JC(x)). Since both x_i and v lie on M, the vector x_{i+1} - v is generally not in Null(JC(x_i)): its normal component is O(||x_{i+1}-v||^2) but the estimate (1/C)||Π Lλ Lx(x)(x_{i+1}-v)|| bounds only the projected, tangent part. Without a separate bound for the normal component, the displayed conclusion ||x_{i+1}-v|| ≤ (1/C)||Π Lλ Lx(x)(x_{i+1}-v)|| does not follow. This gap affects the quadratic convergence claim in Theorem 6.2 and also feeds into the cubic-rate part.
- [§6, proof of Theorem 6.2, Chebyshev part] The proof uses the expansion r(x_i, τ) = x_i + τ + O(||τ||^3) for a second-order retraction. A second-order retraction only guarantees r(x_i, τ) = x_i + τ + (1/2) II_{x_i}(τ,τ) + O(||τ||^3), where II is the second fundamental form and the quadratic term is normal. The substitution x_i - v = (x_{i+1}-v) - τ - O(||τ||^3) ignores this O(||τ||^2) normal term, which enters the Taylor expansion of L at second order. The paper provides no argument that this term cancels with the normal component of x_i - v or is absorbed into the tensor G. Without such an estimate, the cubic convergence proof is incomplete.
- [§6, Eq. (6.15) and the A-term argument] Even after correcting the retraction expansion, the claim that the term A is dominated by Π Lλ Lx( x̂_i)(x_{i+1}-v) requires the operator Π Lλ Lx to be invertible uniformly on the relevant directions; the proof only states a crude bound ||A|| ≤ D||x_{i+1}-v|| ||x_i-v||. This bound also needs the normal component of x_{i+1}-v to be controlled, which is not established. The argument as written does not rigorously justify the transition from (6.14)-(6.15) to the final rate.
minor comments (4)
- [§7.2, complex tensor eigenpairs] The paper states that the algorithm enumerates all complex eigenpairs using random starting points and a deduplication table, but also acknowledges that 'a detailed global convergence analysis is still needed.' This heuristic is not a theorem, and the claim that all pairs are found rests on numerical experience; the presentation should make this limitation more prominent.
- [Throughout] There are numerous typographical and OCR artifacts, such as 'T aylor', 'iterati ons', 'coincident' for 'coincidence', and broken table formatting in §7.2. These should be cleaned up.
- [§7.1, Table] The comparison table for SO-NCM vs O-NCM reports only average times and improvement percentages without standard deviations or hardware details, making the claimed 16% improvement hard to assess; a brief description of the experimental setup would help.
- [Remark 6.3] The reference to 'the notebook TwoLeftInverses.ipynb' is informal; since the numerical experiments are part of the evidence, the manuscript should cite the specific code version or include the relevant results in the paper.
Circularity Check
No significant circularity: Theorem 6.2's convergence rates follow from Taylor expansion under explicitly stated hypotheses, not from fitted inputs or self-citations.
full rationale
The paper's central derivation (Theorem 6.2) is self-contained. The proof starts from the Taylor expansion 0 = L(v, mu) = L(x_i,R(x_i)) + Lx(...)(v-x_i) + Llambda(...)(mu - R(x_i)) + O(||(x_i,mu)-(v,R(x_i))||^2), applies the projection Pi_Llambda to remove the Lagrange multiplier term, and uses the defining Rayleigh equation Pi_Llambda Lx eta = -Pi_Llambda L to cancel the first-order residual; the estimate (6.7) on the projected Hessian then converts residual order into displacement order. The quadratic and cubic rates are obtained directly from the order of the Taylor remainder, with no parameter fitted to the target rates. The generalized Rayleigh quotient R is a hypothesis (R(v)=mu), not an output of the proof, and the 'G=0 gives cubic convergence' condition is computed from second derivatives of L and JR; the classical RQI, two-sided RQI, and Riemannian Newton cases are applications or benchmarks, not inputs. The only self-reference is the author's GitHub repository [22], which supplies numerical code and notebooks; it is not used to establish any theorem and therefore is not load-bearing. Whether the proof's use of (6.7) on x_{i+1}-v is fully justified on a curved constraint (normal-component control) is a correctness question, not a circularity one.
Assumptions & free parameters
assumptions (4)
- domain assumption Full-rank and invertibility hypotheses: L_lambda admits a C1 left inverse L_lambda^- and Pi_Llambda Lx is uniformly invertible on TM near v (inequality 6.7).
- standard math Retractions r of first and second order exist and satisfy xi+1 = r(xi, eta) with the stated tangent and normal behavior.
- domain assumption Smoothness assumptions: L and C are C2 or C3, and the Rayleigh quotient R is C1 or C2 in a neighborhood of the solution.
- ad hoc to paper For complex eigenpair enumeration, random starting points together with a deduplication table eventually visit all equivalence classes of eigenpairs.
Cite this review
Pith. "Pith review of Rayleigh Quotient Iteration, cubic convergence, and second covariant derivative." pith.science (2026). https://pith.science/paper/ANFDCJ23
@misc{pith2026190800639,
author = {Pith},
title = {Pith review of: Rayleigh Quotient Iteration, cubic convergence, and second covariant derivative},
year = {2026},
howpublished = {\url{https://pith.science/paper/ANFDCJ23}},
note = {Machine review of arXiv:1908.00639}
}
read the original abstract
We generalize the Rayleigh Quotient Iteration (RQI) to the problem of solving a nonlinear equation where the variables are divided into two subsets, one satisfying additional equality constraints and the other could be considered as (generalized nonlinear Lagrange) multipliers. This framework covers several problems, including the (linear\slash nonlinear) eigenvalue problems, the constrained optimization problem, and the tensor eigenpair problem. Often, the RQI increment could be computed in two equivalent forms. The classical Rayleigh quotient algorithm uses the {\it Schur form}, while the projected Hessian method in constrained optimization uses the {\it Newton form}. We link the cubic convergence of these iterations with a {\it constrained Chebyshev term}, showing it is related to the geometric concept of {\it second covariant derivative}. Both the generalized Rayleigh quotient and the {\it Hessian of the retraction} used in the RQI appear in the Chebyshev term. We derive several cubic convergence results in application and construct new RQIs for matrix and tensor problems.
Reference graph
Works this paper leans on
-
[1]
P. Absil, R. Mahony, R. Sepulchre, and P. V an Dooren , A Grassmann–Rayleigh quotient iteration for computing invariant subspaces, SIAM Rev., 44 (2002), pp. 57–73, https://doi.org/10.1137/S0036144500378648
-
[2]
P. Absil and J. Malick , Projection-like retractions on matrix manifolds , SIAM J. Optim, 22 (2012), pp. 135–158, https://doi.org/10.1137/100802529
- [3]
- [4]
-
[5]
P.-A. Absil and P. V an Dooren , Two-sided Grassmann-Rayleigh quotient iteration , Numer. Math., 114 (2010), pp. 549–571, https://doi.org/10.1007/s00211-009-0266-y
- [6]
-
[7]
S. Boyd and L. V andenberghe , Convex Optimization , Cambridge University Press, 2004, https://doi.org/10.1017/ CBO9780511804441
work page 2004
-
[8]
V. Candela and A. Marquina , Recurrence relations for rational cubic methods ii: The Che byshev method , Com- puting, 45 (1990), pp. 355–367, https://doi.org/10.1007/BF02238803
Show all 30 references
-
[9]
Cartwright and B
D. Cartwright and B. Sturmfels , The number of eigenvalues of a tensor , Linear Algebra and its Applications, 438 (2013), pp. 942 – 952, https://doi.org/https://doi.org/10.1016/j.laa.2011.05.040. Tensors and Multilinear Algebra
2013 doi
-
[10]
L. Chen, L. Han, and L. Zhou , Computing tensor eigenvalues via homotopy methods , SIAM Journal on Matrix Analysis and Applications, 37 (2016), pp. 290–319, https://doi.org/10.1137/15M1010725
2016 doi
-
[11]
Edelman, T
A. Edelman, T. A. Arias, and S. T. Smith , The geometry of algorithms with orthogonality constraints , SIAM J. Matrix Anal. Appl., 20 (1999), pp. 303–353, https://doi.org/10.1137/S0895479895290954
1999 doi
-
[12]
Gabay , Minimizing a differentiable function over a differential man ifold, Journal of Optimization Theory and Applications, 37 (1982), pp
D. Gabay , Minimizing a differentiable function over a differential man ifold, Journal of Optimization Theory and Applications, 37 (1982), pp. 177–219, https://doi.org/10.1007/BF00934767
1982 doi
-
[13]
Gabay , Reduced quasi-Newton methods with feasibility improvemen t for nonlinearly constrained optimization , Springer Berlin Heidelberg, Berlin, Heidelberg, 1982, pp
D. Gabay , Reduced quasi-Newton methods with feasibility improvemen t for nonlinearly constrained optimization , Springer Berlin Heidelberg, Berlin, Heidelberg, 1982, pp. 18–44, https://doi.org/10.1007/BFb0120946
1982 doi
-
[14]
Gabay and D
D. Gabay and D. Luenberger , Efficiently converging minimization methods based on the red uced gradient, SIAM J. Control Optim., 14 (1976), pp. 42–61, https://doi.org/10.1137/0314004
1976 doi
-
[15]
Güttel and F
S. Güttel and F. Tisseur , The nonlinear eigenvalue problem , Acta Numer., 26 (2017), p. 1–94, https://doi.org/10. 1017/S0962492917000034
2017
-
[16]
Jaffe, R
A. Jaffe, R. Weiss, and B. Nadler , Newton correction methods for computing real eigenpairs of symmetric ten- sors, SIAM Journal on Matrix Analysis and Applications, 39 (2018 ), pp. 1071–1094, https://doi.org/10.1137/ 17M1133312
2018
-
[17]
T. G. Kolda and J. R. Mayo , Shifted power method for computing tensor eigenpairs , SIAM Journal on Matrix Analysis and Applications, 32 (2011), pp. 1095–1124, https://doi.org/10.1137/100801482
2011 doi
-
[18]
Lancaster , Some applications of the Newton-Raphson method to non-line ar matrix problems , Proc
P. Lancaster , Some applications of the Newton-Raphson method to non-line ar matrix problems , Proc. Roy. Soc. London. Ser. A., 271 (1963), pp. 324–331, https://doi.org/10.1098/rspa.1963.0021
1963
-
[19]
L. H. Lim , Singular values and eigenvalues of tensors: a variational a pproach, in 1st IEEE International W orkshop on Computational Advances in Multi-Sensor Adaptive Proces sing, 2005., Dec 2005, pp. 129–132, https://doi.org/ 10.1109/CAMAP.2005.1574201
2005
-
[20]
J. H. Manton , Optimization algorithms exploiting unitary constraints , IEEE Transactions on Signal Processing, 50 (2002), pp. 635–650, https://doi.org/10.1109/78.984753
2002 doi
-
[21]
Mishra and R
B. Mishra and R. Sepulchre , Riemannian preconditioning, SIAM Journal on Optimization, 26 (2016), pp. 635–660, https://doi.org/10.1137/140970860
2016 doi
-
[22]
Nguyen , Project lagrange_rayleigh
D. Nguyen , Project lagrange_rayleigh. https://github.com/dnguyend/lagrange_rayleigh, 2019
2019
-
[23]
G. Ni, L. Qi, F. W ang, and Y. W ang , The degree of the e-characteristic polynomial of an even ord er tensor , Journal of Mathematical Analysis and Applications, 329 (20 07), pp. 1218 – 1229, https://doi.org/https://doi. org/10.1016/j.jmaa.2006.07.064
-
[24]
A. M. Ostrowski , On the convergence of the Rayleigh quotient iteration for th e computation of the characteristic roots and vectors. iii , Arch. Ration. Mech. Anal., 3 (1959), pp. 325–340, https://doi.org/10.1007/BF00284184
1959 doi
-
[25]
M. J. D. Powell , Algorithms for nonlinear constraints that use lagrangian f unctions, Math. Program., 14 (1978), pp. 224–248, https://doi.org/10.1007/BF01588967
1978 doi
-
[26]
Qi, Eigenvalues and invariants of tensors , Journal of Mathematical Analysis and Applications, 325 (2 007), pp
L. Qi, Eigenvalues and invariants of tensors , Journal of Mathematical Analysis and Applications, 325 (2 007), pp. 1363 – 1377, https://doi.org/https://doi.org/10.1016/j.jmaa.2006.02.071
2006 doi
-
[27]
K. Schreiber , Nonlinear eigenvalue problems: Newton-type methods and no nlinear Rayleigh functionals , PhD thesis, Technische Universitat Berlin, 2008, https://depositonce.tu-berlin.de/bitstream/11303/2168/2/Dokument_42. pdf
2008
-
[28]
R. A. Tapia , Diagonalized multiplier methods and quasi-newton methods for constrained optimization , J. Optim. Theory Appl., 22 (1977), pp. 135–194, https://doi.org/10.1007/BF00933161
1977 doi
-
[29]
Townsend, N
J. Townsend, N. Koep, and S. Weichw ald , Pymanopt: A python toolbox for optimization on manifolds us ing automatic differentiation , J. Mach. Learn. Res., 17 (2016), pp. 1–5, http://jmlr.org/papers/v17/16-177.html
2016
-
[30]
Z. Zhao, Z. Bai, and X. Jin , A Riemannian Newton algorithm for nonlinear eigenvalue pro blems, SIAM J. Matrix Anal. Appl., 36 (2015), pp. 752–774, https://doi.org/10.1137/140967994. 25 This manuscript is for review purposes only
2015 doi
Reviewed August 14, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.