Rayleigh Quotient Iteration, cubic convergence, and second covariant derivative

TL;DR

This paper extends Rayleigh Quotient Iteration to nonlinear equations with constraints, linking its cubic convergence to second covariant derivatives and applying it to matrix and tensor eigenproblems.

Contribution

It generalizes RQI for constrained nonlinear problems, establishes a geometric interpretation of cubic convergence, and develops new RQIs for matrices and tensors.

Findings

Established cubic convergence linked to second covariant derivatives.

Derived new RQIs for matrix and tensor eigenproblems.

Connected classical RQI forms with geometric concepts.

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.

Tables

Table 1. Table 1. Examples of residual errors for Rayleigh (RQI) and Rayleigh-Chebyshev (RC) iterations for tensor eigenpairs (Tensor, section   6.1.1 , n = 6 , m = 3 , d = 2 formulae-sequence 𝑛 6 formulae-sequence 𝑚 3 𝑑 2 n=6,m=3,d=2 ) and eigenvalue with a constant term (EigC) ( section   6.3 , n = 10 𝑛 10 n=10 ).

Iteration	Tensor RQI	Tensor RC	EigC RQI	EigC RC
1	2.594e-04	2.042e-05	6.986e-06	6.986e-06
2	3.080e-07	8.671e-14	2.368e-11	6.511e-15
3	2.509e-13	2.690e-39	3.500e-22	5.208e-42
4	1.049e-25	5.503e-115	9.982e-44

Equations

$$\mathcal{L}(X,\lambda )=\left[\begin{array}{c}\mathrm{L}(X,\lambda )\\ \mathrm{C}(X)\end{array}\right].$$

$$\mathrm{L}(X,\lambda )=$$

$$\mathrm{C}(X)=$$

$$\mathrm{L}(X,\lambda )=\nabla f(X)-{\mathrm{C}}^{\prime }(X{)}^{\mathsf{T}}\lambda =0.$$

$$\mathrm{L}(X,\lambda )=\mathrm{P}(\lambda )X=0.$$

$$\mathrm{L}(X,\lambda )=$$

$$X_{i+1}=\frac{(A-{\lambda }_{i}I{)}^{-1}{X}_i}{||(A-{\lambda }_{i}I{)}^{-1}{X}_i||}.$$

$$\mathfrak{r}(X,\eta )=X+\eta +O(|\eta {|}^2)$$

$$\mathrm{F}(X):=\mathrm{L}(X,X^{\mathsf{T}}AX)=AX-X{X}^{\mathsf{T}}AX=(I-X{X}^{\mathsf{T}})AX.$$

$$\Pi (X)\mathrm{F}(X)=\mathrm{F}(X).$$

$$\Pi (X_i){\mathrm{F}}^{\prime }(X_i){\eta }_i=-\mathrm{F}(X_i)$$

$${\mathrm{C}}^{\prime }(X_i){\eta }_i=0$$

$$\Pi (X)\mathrm{L}(X,\lambda )=0,$$

$$\mathrm{L}(X,\lambda )-\Pi (X)\mathrm{L}(X,\lambda )=0.$$

$${\mathrm{L}}_{\lambda }^{-}(X,\lambda )\mathrm{L}(X,\lambda )=0.$$

$$\Pi (X)\mathrm{L}(X,\mathcal{R}(X))=\mathrm{L}(X,\mathcal{R}(X)),$$

$$\Pi (X){\mathrm{L}}_{\lambda }(X,\mathcal{R}(X))=0.$$

$$\Pi (X){\mathrm{F}}^{\prime }(X)=\Pi (X)\{{\mathrm{L}}_X(X,\mathcal{R}(X))+{\mathrm{L}}_{\lambda }(X,\mathcal{R}(X)){\mathcal{R}}^{\prime }(X)\}=\Pi (X){\mathrm{L}}_X(X,\mathcal{R}(X)).$$

$$\begin{array}{c}PB\eta =-PF,\\ D\eta =0.\end{array}$$

$$\eta =-B^{-1}F+B^{-1}H(D{B}^{-1}H{)}^{-1}D{B}^{-1}F.$$

$$\eta =-X+\zeta (X^{\mathsf{T}}\zeta {)}^{-1}{X}^{\mathsf{T}}X=-X+(X^{\mathsf{T}}\zeta {)}^{-1}(A-\lambda I{)}^{-1}X$$

$$\eta =-Q_T(Q_{\Pi }^{\mathsf{T}}\Pi {\mathrm{L}}_X{Q}_T{)}^{-1}{Q}_{\Pi }^{\mathsf{T}}\mathrm{L},$$

$${\nabla }_{\eta }\mathrm{F}(X):=\Pi (X){\mathrm{F}}^{\prime }(X;\eta )={\mathrm{F}}^{\prime }(X;\eta )-{\Pi }^{\prime }(X,\eta ))\mathrm{F}(X).$$

$$\begin{array}{c}\mathrm{G}(X;{\eta }^{[2]})={\mathrm{F}}^{(2)}(X;{\eta }^{[2]})+{\mathrm{F}}^{\prime }(X){\mathfrak{r}}_{\eta \eta }(X,0;{\eta }^{[2]}).\end{array}$$

$$\Pi (X_i){\mathrm{F}}^{\prime }(X_i)X_i^{\tau }=\Pi (X_i)\mathrm{G}(X_i){\eta }_i^{[2]}$$

$$\begin{array}{c}{\mathrm{G}}_{\mathrm{L}}(X;{\eta }^{[2]})={\mathrm{L}}_{XX}(X,\lambda ){\eta }^{[2]}+2{\mathrm{L}}_{X\lambda }(X,\lambda )[\eta ,{\mathcal{R}}^{\prime }(X;\eta )]+{\mathrm{L}}_{\lambda \lambda }(X,\lambda ){\mathcal{R}}^{\prime }(X;\eta {)}^{[2]}+{\mathrm{L}}_X(X,\lambda ){\mathfrak{r}}_{\eta \eta }(X,0;{\eta }^{[2]}).\end{array}$$

$$\begin{array}{c}\mathrm{F}(X)=\mathrm{F}(X_0)+\sum _{j=1}^l\frac{1}{j!}{\mathrm{F}}^{(j)}(X_0;[X-X_0{]}^{[j]})+{\int }_0^1\frac{(1-t{)}^l}{l!}{\mathrm{F}}^{(l+1)}((1-t)X_0+tX;[X-X_0{]}^{[l+1]})dt\end{array}$$

$$\begin{array}{c}\mathrm{C}(X)=0,\\ {\mathrm{C}}^{\prime }(X)\eta =0.\end{array}$$

$$D_{\varphi }:(z,v)\mapsto (\varphi (z),{\varphi }^{\prime }(z)v)\in T\mathcal{M}$$

$$\mathtt{T}{\Omega }_{{\rho }_{\mathfrak{r}}}^{\mathfrak{r}}:=\{(X,\eta )|(X,\eta )\in T\mathcal{M},X\in {\Omega }^{\mathfrak{r}},|\eta |<{\rho }_{\mathfrak{r}}\}.$$

$$\mathtt{T}{\Omega }_{\varphi ,{\rho }_1,{\rho }_2}^{\mathfrak{r}}:=\{(\varphi (z),{\varphi }^{\prime }(z)v)||z|<{\rho }_1,|v|<{\rho }_2\}$$

$$A^{(2)}(X;{\eta }^{[2]})+A^{\prime }(X)\mathfrak{r}(X,.{)}^{(2)}(0;{\eta }^{[2]})$$

Full text

1. Introduction

Several important problems in mathematics could be reduced to solving a vector equation of the form $\mathcal{L}(X,\lambda)=0$, where the function $\mathcal{L}$ and the unknowns are divided into two groups as follows

[TABLE]

Here, $X$ and $\lambda$ are vector variables, defined on two vector spaces $\mathcal{E}$ and $\mathcal{E}_{L}$, respectively. The function $\mathrm{L}(X,\lambda)=0$ in the first equation involves all variables, $\mathrm{L}$ maps $E\times\mathcal{E}_{L}$ to $\mathcal{E}$. In the second group, the constraint $\mathrm{C}(X)=0$, involves only $X$, $\mathrm{C}$ maps $\mathcal{E}$ to the constraint vector space $\mathcal{E}_{L}$. The vector variable $\lambda\in\mathcal{E}_{L}$ plays the role of (generalized Lagrange) multipliers. Thus $\mathcal{L}$ maps $E\times\mathcal{E}_{L}$ to itself. This setup covers at least four classes of problems encountered in the literature:

a) The eigenvector problem:

[TABLE]

where $X$ is an $n\times 1$ vector, $A$ is an $n\times n$ matrix, $\lambda$ is a scalar. Here, $\mathcal{E}$ is the space of $n\times 1$ vectors and $\mathcal{E}_{L}$ is the base field, $\mathbb{R}$. A related problem has $\mathrm{L}(X,\lambda)=AX-X\lambda-\mathtt{b}=0$ where $\mathtt{b}\neq 0$ is a vector ([1] and references therein).

b) The constrained optimization problem. This is the problem of optimizing a cost function $f$ in a space $E$ under the constraint $\mathrm{C}(X)=0$ where $\mathrm{C}$ is a function with target space $\mathcal{E}_{L}$ and $\lambda\in\mathcal{E}_{L}$ are the Lagrange multipliers. Let $\mathrm{C}^{\prime}$ be the Jacobian of $\mathrm{C}$ then the Lagrangian multiplier equation is

[TABLE]

Note $\mathrm{C}^{\prime}(X)^{\mathsf{T}}$ is a map from $\mathcal{E}_{L}$ to $\mathcal{E}$. The system eq. 1.1 gives us the set of critical points.

c) The nonlinear eigenvalue problem:

[TABLE]

In the real case, we set $\mathcal{E}_{L}=\mathbb{R}$, with $\lambda\in\mathbb{R}$ is scalar. Here $\mathrm{P}$ is a matrix with polynomial entries in $\lambda$. While this is not in the form of eq. 1.1 we can impose the constraint $\mathrm{C}(X)=X^{\mathsf{T}}X-1$ (or $\mathrm{C}(X)=z^{\mathsf{T}}X-1$ for a fixed vector $z$). See [2] for a survey.

d) The tensor eigenpairs problem: Here, $\mathcal{E}=\mathbb{F}^{n}$ is a vector space over a base field $\mathbb{F}$, $\mathbb{F}=\mathbb{R}$ or $\mathbb{C}$. Set $\mathcal{E}_{L}=\mathbb{F}$ , thus $\lambda$ is scalar. Let $\mathcal{T}$ and $\mathcal{B}$ be two vector-valued functions from $\mathcal{E}$ to itself, with entries $\mathcal{T}_{i}$ and $\mathcal{B}_{i},1\leq i\leq n$ are homogeneous polynomials of degrees $m-1$ and $d-1$ respectively. The evaluations of $\mathcal{T}$ and $\mathcal{B}$ at $X\in\mathcal{E}$ are written $\mathcal{T}(X^{[m-1]})$ and $\mathcal{B}(X^{[d-1]})$. Set

[TABLE]

A popular constraint is $\mathrm{C}(X)=\frac{1}{2}(X^{\mathsf{T}}X-1)$ for the real case, and we will study $\mathrm{C}(X)=\frac{1}{2}(X^{*}X-1)$ for the complex case ($*$ is the Hermitian transpose). An important case is when $\mathcal{T}=1/m\hat{\mathcal{T}}^{\prime}$, where $\hat{\mathcal{T}}^{\prime}$ denotes the gradient of a scalar homogeneous polynomial $\hat{\mathcal{T}}$ of order $m$, and $\mathcal{B}(X)=X$. These eigenpairs could be used to determine if $\hat{\mathcal{T}}$ is nonnegative ([3]).

The Rayleigh quotient iteration (RQI) is among the most powerful methods to compute eigenvalues and vectors. For a vector $X$, the Rayleigh quotient is $\lambda=\frac{X^{\mathsf{T}}AX}{X^{\mathsf{T}}X}$, or $X^{\mathsf{T}}AX$ on the unit sphere. In the $i$-th step, with $\lambda_{i}$ computed from $X_{i}$ by this equation, the iteration computes

[TABLE]

It has cubic convergence when $A$ is normal and quadratic otherwise (for suitable initial points).

Similar iterations have been suggested for the remaining problems. Points satisfying $\mathrm{C}(X)=0$ are called feasible points. At the $i$-th step, where $X=X_{i}$ is feasible, a (vector-valued) function $\mathcal{R}(X)$ of $X$ is used to compute $\lambda$ for the iterative step, then an intermediate step $\hat{X}_{i+1}$ is produced by solving a linear equation depending on $\lambda$, and a feasibility perturbation is applied to produce the next feasible step $X_{i+1}$. The feasibility perturbation is often available for common constraints, thus the overall procedure consists of constructing the generalized Rayleigh quotient $\mathcal{R}$ and the iteration equation. In the literature, convergence results are obtained separately for each problem. We show here there is a common procedure and convergence analysis for these problems, including criteria for cubic convergence.

Recall an iteration $\{X_{i}\}$ converging to $X_{*}$ has order $k$ if $\lvert X_{i+1}-X_{*}\rvert\leq M\lvert X_{i}-X_{*}\rvert^{k}$ for some $M$, for large enough $i$. Quadratic and cubic convergence corresponds to $k=2$ and $k=3$. It is well-known Newton’s method has quadratic convergence order and the Chebyshev method [4] achieves cubic convergence. We link the cubic convergence of the Rayleigh quotient iteration with a constrained Chebyshev term. While very little geometric prerequisite is needed in the proof or in applications, we explain this Chebyshev term is essentially a second covariant derivative using two connections constructed from the problem formulation.

1.1. Generalized RQI: quadratic convergence

We denote the Jacobian of a nonlinear map $\mathrm{F}$ between vector spaces $V$ and $W$ by $\mathrm{F}^{\prime}$. For $X\in V$, $\mathrm{F}^{\prime}(X)$ is a linear operator between $V$ and $W$. Let $\mathtt{Lin}(V,W)$ be the space of linear operators from $V$ to $W$. We will assume sufficient smoothness. Consider a constraint $\mathrm{C}(X)=0$ defining a feasible set $\mathcal{M}$ such that $\mathrm{C}^{\prime}(X)$ is of full rank in the subset of interest, $\mathcal{M}$ is equipped with a feasibility perturbation $\mathfrak{r}$. Recall the tangent space of $\mathcal{M}$ at $X\in\mathcal{M}$ is the nullspace $T_{X}\mathcal{M}$, or $T_{X}$ for short, of $\mathrm{C}^{\prime}(X)$. We are mainly interested in feasibility perturbation from the tangent space. For $X\in\mathcal{M},\eta\in T_{X}$, instead of considering a map from $\mathcal{E}$ to $\mathcal{M}$, mapping $X+\eta$ to a point on $\mathcal{M}$, following [5], we consider a map $\mathfrak{r}$ sending $(X,\eta)$ to $\mathfrak{r}(X,\eta)\in\mathcal{M}$ for $\eta$ sufficiently small, with $\mathfrak{r}(X,0)=X$, ($0=0_{T_{X}}\in T_{X}$ is the zero vector). For the constraint $X^{\mathsf{T}}X=1$, we can use the rescaling $\mathfrak{r}(X,\eta)=\frac{1}{|X+\eta|}(X+\eta)=\frac{1}{(1+|\eta|^{2})^{1/2}}(X+\eta).$

Fix $X$, then $\mathfrak{r}(X,.):\eta\mapsto\mathfrak{r}(X,\eta)$ is a map from $T_{X}$ to $\mathcal{E}$. Assume it has the Taylor expansion

[TABLE]

or $\mathfrak{r}$ is a retraction, thus, the partial derivative $\mathfrak{r}_{\eta}$ of $\mathfrak{r}$ in $\eta$, considered as a linear map from $T_{X}$ to $\mathcal{E}$ satisfies $\mathfrak{r}_{\eta}(X,0)\eta=\eta$. This holds for the rescaling perturbation of the sphere and most feasibility perturbations in practice, and we can normalize a reasonable perturbation to this form.

The second ingredient of a generalized RQI is a projection function. In eqs. 1.2 and 1.5 of the classical RQI, with $\lambda=X^{\mathsf{T}}AX$

[TABLE]

The matrix $\Pi(X)=I-XX^{\mathsf{T}}$ is a projection, as $(I-XX^{\mathsf{T}})^{2}=I-XX^{\mathsf{T}}$. Thus $\Pi(X)\mathrm{F}(X)=\Pi(X)^{2}AX=\mathrm{F}(X)$. This seems incidental, but it is crucial in generalizing the RQI.

An operator $P\in\mathtt{Lin}(\mathcal{E},\mathcal{E})$ is called an affine projection if $P^{2}=P$. It is often constructed in two ways. If $A^{-}$ is a left inverse of a linear operator $A$, $A^{-}A=I$, then $AA^{-}$ is a projection. If $P$ is a projection then $I_{\mathcal{E}}-P$ is also a projection. If $\eta\in\mathcal{E}$ is in the range $\operatorname{Im}(P)$ of $P$ then $P\eta=\eta$. In the special case where $P$ is self-adjoint, $P$ is called an orthogonal projection. We call a smooth map $\Pi$ from $\mathcal{M}$ to $\mathtt{Lin}(\mathcal{E},\mathcal{E})$, the space of linear operators on $\mathcal{E}$ an (affine) projection function, or just projection when there is no confusion, if, for $X\in\mathcal{M}$, $\Pi(X)$ is an affine projection.

We will assume for a projection function $\Pi$, $\Pi(X)$ is of rank $\dim\mathcal{E}-\dim\mathcal{E}_{L}$ for $X\in\mathcal{M}$.

The following claim is the key to our generalized RQI to solve eq. 1.1.

Claim 1**.**

Consider a smooth feasible set $\mathcal{M}$ defined by $\mathrm{C}(X)=0$, $\mathrm{C}^{\prime}(X)$ is of full rank for $X\in\mathcal{M}$, and $\Pi:\mathcal{M}\to\mathtt{Lin}(\mathcal{E},\mathcal{E})$ is a projection function. Let $\mathrm{F}$ be a smooth function from $\mathcal{M}$ to $\mathcal{E}$ satisfying

[TABLE]

For an initial point $X_{0}\in\mathcal{M}$ sufficiently close to a solution $X_{*}$ of $\mathrm{F}(X)=0$, the iteration $X_{i+1}=\mathfrak{r}(X_{i},\eta_{i})$ with $\mathfrak{r}$ satisfies eq. 1.6 and $\eta_{i}\in T_{X_{i}}\mathcal{M}=T_{X_{i}}$ solving the system (1.8), (1.9) below

[TABLE]

is well-defined and converges quadratically to $X_{*}$ if the system is non-degenerated at $X_{*}$.

To use this claim to solve eq. 1.1, note for any projection $\Pi$, $\mathrm{L}(X,\lambda)=0$ is equivalent to

[TABLE]

If a function $\mathcal{R}(X)$ substituted to $\lambda$ satisfying eq. 1.11 for all $X$, consider $\mathrm{F}(X):=\mathrm{L}(X,\mathcal{R}(X))$, then $\Pi(X)\mathrm{F}(X)=\mathrm{F}(X)$ is satisfied by design, and we can use the iteration in the claim.

To simplify eq. 1.8, we construct $\Pi$ as below. Let $\mathrm{L_{\lambda}}(X,\lambda)$ be the partial derivative of $\mathrm{L}$ in $\lambda$, considered as a linear map from $\mathcal{E}_{L}$ to $\mathcal{E}$. Assume it is onto, which is the case when the inverse function theorem is satisfied for eq. 1.1, then it has a left inverse $\mathrm{L}_{\lambda}^{-}(X,\lambda)$ (any computationally efficient inverse or the Moore-Penrose inverse can be used). Solve for $\lambda=\mathcal{R}(X)$ from

[TABLE]

Then with $\Pi^{-}=\mathrm{L_{\lambda}}\mathrm{L}_{\lambda}^{-}$, $\Pi(X):=I_{\mathcal{E}}-\Pi^{-}(X,\mathcal{R}(X))$ is a projection, eq. 1.11 becomes eq. 1.13 and eq. 1.14 holds as $\mathrm{L_{\lambda}}=\mathrm{L_{\lambda}}\mathrm{L}_{\lambda}^{-}\mathrm{L_{\lambda}}$ by the left inverse assumption

[TABLE]

From eq. 1.14, the expression $\Pi(X)\mathrm{F}^{\prime}(X)$ in claim 1 with $\mathrm{F}(X)=\mathrm{L}(X,\mathcal{R}(X))$ reduces to

[TABLE]

Here $\mathrm{L}_{X}(X,\lambda)$ is the partial derivative in $X$. We will choose $\mathcal{R}$ solving eq. 1.12.

For the eigenvalue problem $\mathrm{L}(X,\lambda)=AX-X\lambda$ with constraint $X^{\mathsf{T}}X=1$, $\mathrm{L_{\lambda}}(X,\lambda)=-X$, take the left inverse $\mathrm{L}_{\lambda}^{-}=-X^{\mathsf{T}}$, then eq. 1.12 gives us $\lambda=X^{\mathsf{T}}AX$ and $\Pi(X)=I_{\mathcal{E}}-XX^{\mathsf{T}}$. Using the iteration in claim 1 for $\mathrm{F}(X)=AX-XX^{\mathsf{T}}AX$, we get the Rayleigh quotient iteration, since

Lemma 1.1** (Schur form).**

Let $P\in\mathtt{Lin}(\mathcal{E},\mathcal{E})$ be an affine projection, $B\in\mathtt{Lin}(\mathcal{E},\mathcal{E})$ be an invertible operator, $F\in\mathcal{E}$, and $D\in\mathtt{Lin}(\mathcal{E},\mathcal{E}_{L})$. Consider the linear equations in $\eta$

[TABLE]

If $H\in\mathtt{Lin}(\mathcal{E}_{L},\mathcal{E})$ is such that $PH=0$ and $DB^{-1}H$ is invertible, $\eta$ below solves eq. 1.15

[TABLE]

The proof is by direct substitution, see section A.1.

Continuing with the eigenvalue problem, with $\lambda=X^{\mathsf{T}}AX$, $\mathrm{F}(X)=AX-XX^{\mathsf{T}}AX$, $\mathrm{L}_{X}(X,\lambda)\eta=(A-\lambda I_{\mathcal{E}})\eta$ for $\eta\in\mathcal{E}$. In the lemma, set $B=A-\lambda I_{\mathcal{E}}$, $H:\lambda\mapsto X\lambda$, $D:\eta\mapsto X^{\mathsf{T}}\eta$, $F=AX-X\lambda$, then $B^{-1}F=X$ and $\zeta:=B^{-1}H=(A-\lambda I)^{-1}X$, hence

[TABLE]

with $(X^{\mathsf{T}}\zeta)^{-1}$ is a scalar. Then $X+\eta$ is proportional to $(A-\lambda I)^{-1}X$, and the retraction gives us the iteration $X_{i+1}=\|(A-\lambda I)^{-1}X_{i}\|^{-1}(A-\lambda I)^{-1}X_{i}$.

This “explains” why the Rayleigh quotient iteration has quadratic convergence and relates it to the Newton-type iteration of claim 1. This explanation for the classical Rayleigh quotient essentially appeared in [5, example 5]. Note $\Pi$ is hidden in the classical RQI.

In lemma 1.1, take $H=\mathrm{L_{\lambda}}(X,\mathcal{R}(X))$, we have the Schur form for claim 1. We will apply this result to the remaining three problems in section 6. The Schur form’s usage is restricted by the condition $\mathrm{L}_{X}(X,\mathcal{R}(X))$ is invertible. In optimization problems (section 6.2), eq. 1.8 is often solved directly. Let $Q_{\Pi}$ be an orthonormal basis of the nullspace of $\mathrm{L}_{\lambda}^{-}(X,\lambda)$ (same as the range of $\Pi(X)$), and $Q_{T}$ an orthonormal basis of $T_{X}$, then expressing $\mathrm{L}$ and $\Pi\mathrm{L}_{X}$ in these bases (evaluated at $(X,\lambda)$) we have

[TABLE]

the Newton form of the iteration. Note $Q_{\Pi}Q_{\Pi}^{\mathsf{T}}$ and $\Pi(X)$ are both projections to the range of $\Pi(X)$, hence, $Q_{\Pi}Q_{\Pi}^{\mathsf{T}}\Pi=\Pi,Q_{\Pi}Q_{\Pi}^{\mathsf{T}}\mathrm{L}=\mathrm{L}$ which helps verify $\Pi\mathrm{L}_{X}\eta=Q_{\Pi}Q_{\Pi}^{\mathsf{T}}\Pi\mathrm{L}_{X}\eta=-Q_{\Pi}Q_{\Pi}^{\mathsf{T}}\mathrm{L}=-\mathrm{L}$.

To summarize, we have a procedure to solve eq. 1.1 by a generalized RQI

1. Find a left inverse $\mathrm{L}_{\lambda}^{-}(X,\lambda)$ of $\mathrm{L}_{\lambda}(X,\lambda)$ such that $\mathrm{L}_{\lambda}^{-}(X,\lambda)\mathrm{L}(X,\lambda)=0$ is easy to solve for $\lambda$ as a function $\lambda=\mathcal{R}(X)$ in $X$. Choose a retraction $\mathfrak{r}$.

2. Determine if the Schur form or the Newton form solution of eq. 1.8 is preferable. For the Newton form, also set up $\Pi(X)=I_{E}-\mathrm{L}_{\lambda}(X,\lambda)\mathrm{L}_{\lambda}^{-}(X,\lambda)$ with $\lambda=\mathcal{R}(X)$.

3. Apply algorithm 1.

1.2. Geometric interpretation

The discussion in this section and section 5 are not required to follow the rest of the paper. However, we believe they are interesting interpretations of the results in geometric terms, clarifying several geometric concepts that are often presented more abstractly.

The projection function $\Pi$ on the feasible set $\mathcal{M}$ associates to each feasible point a vector space $\mathcal{E}_{X}:=\operatorname{Im}(\Pi(X))$. The collections $\mathcal{E}_{\Pi}=\cup_{X\in\mathcal{M}}\mathcal{E}_{X}$ and $T\mathcal{M}=\cup_{X\in\mathcal{M}}T_{X}$ are vector bundles in the geometric literature, assuming appropriate smoothness and constant rank of $\mathcal{E}_{X}$ and $T_{X}$. They could be considered smooth subsets of $\mathcal{E}^{2}$ of pairs $(X,\omega)$ with $X$ satisfying a nonlinear constraint $\mathrm{C}(X)=0$ not involving $\omega$, while the constraint on $\omega$ is linear. We have $\mathcal{E}_{\Pi}=T\mathcal{M}$ in the eigenvalue problem, but they are different for the RQI in section 6.1.1.

A function $\mathrm{F}$ from $\mathcal{M}$ to $\mathcal{E}$ satisfying the condition $\Pi(X)\mathrm{F}(X)=\mathrm{F}(X)$ is called a section of the bundle $\mathcal{E}_{\Pi}$, we require $\mathrm{F}(X)\in\mathcal{E}_{X}$. Our RQI framework could be understood as a construction of a projection function $\Pi$ defining a vector bundle $\mathcal{E}_{\Pi}$, (and the connection $\nabla$ below) and a section $\mathrm{F}(X)=\mathrm{L}(X,\mathcal{R}(X))$. A section $\mathtt{c}$ of the tangent bundle (a function to $\mathcal{E}$ with $\mathtt{c}(X)\in T_{X}$) is called a vector field, each $\mathtt{c}(X),X\in\mathcal{M}$ is a tangent vector at $X$.

For a section $\mathrm{F}$ and a tangent vector $\eta\in T_{X},X\in\mathcal{M}$, we define

[TABLE]

Here, $\Pi^{\prime}(X,\eta)$ and $\mathrm{F}^{\prime}(X,\eta)$ are directional derivatives of $\Pi$ and $\mathrm{F}$ in direction $\eta$, the second equality follows by differentiating $\Pi(X)\mathrm{F}(X)=\mathrm{F}(X)$ in direction $\eta$. If $\mathtt{c}$ is a vector field then $\nabla_{\mathtt{c}}\mathrm{F}:X\mapsto\Pi(X)\mathrm{F}^{\prime}(X;\mathtt{c}(X))$ is a section of $\mathcal{E}_{\Pi}$. The expression for $\nabla_{\mathtt{c}}\mathrm{F}$ shows it is a covariant derivative, or connection in the differential geometric sense [6, section 5.2]. Intuitively, a covariant derivative is a rule to take directional derivative of sections, resulting in sections. The Newton increment equation reads $\nabla_{\eta}\mathrm{F}(X)=-\mathrm{F}(X)$.

Thus, claim 1 is a Newton method on a vector bundle, see [5, 7, 8] and references therein. Newton method on the tangent bundle is well-studied in the Riemannian optimization literature. The novelty here is the identification of a generalized Rayleigh quotient with the construction of a vector bundle and a connection from the data of eq. 1.1 using a left-inverse, even for nonlinear multipliers, and the convergence analysis with a retraction (see [5, 6] for the tangent bundle case), versus geodesics in [7]. The Schur form iteration for the general case is also new.

1.3. Cubic convergence and Chebyshev iterations

For normal matrices, the classical RQI has cubic convergence and the two-sided iterations in [9, 10] also converge cubically for any matrix. The second goal of the paper is to clarify the condition for cubic convergence. We will define notations more properly in section 1.4, but briefly, the Hessian of a function $\mathrm{F}$ could be considered as a map valued in bilinear functions, we use the notation $\mathrm{F}^{(2)}(X;\eta^{[2]})=\mathrm{F}^{(2)}(X)\eta^{[2]}$ to denote the evaluation of this map at $X$, and of the bilinear function in both linear variables at $\eta$. The (partial) Hessian $\mathfrak{r}_{\eta\eta}(X,0;\eta^{[2]})$ of the retraction $\mathfrak{r}$ in the variable $\eta$ at $(X,0)$ and the directional derivative $\mathcal{R}^{\prime}(X;\eta)$ also appear in the analysis.

Claim 2**.**

For a function $\mathrm{F}$ satisfying $\Pi(X)\mathrm{F}(X)=\mathrm{F}(X)$ as in claim 1, and for a tangent vector $\eta\in T_{X}$ at $X$, set

[TABLE]

If $\Pi(X_{*})\mathrm{G}(X_{*};\eta^{[2]})=0$ for all $\eta\in T_{X_{*}}$ at a solution $X_{*}$, the RQI in claim 1 converges cubically to $X_{*}$ when $X_{0}$ is sufficiently close to $X_{*}$. If $\eta_{i}=X_{i}^{\mathcal{N}_{T}}$ is the RQI increment in claim 1, then the Rayleigh-Chebyshev iteration $X_{i+1}=\mathfrak{r}(X_{i},X_{i}^{\mathcal{N}_{C}})$ with $X_{i}^{\mathcal{N}_{C}}=\eta_{i}-\frac{1}{2}X_{i}^{\tau}\in T_{X}$ and $X_{i}^{\tau}\in T_{X_{i}}$ satisfies

[TABLE]

converges cubically for an initial point $X_{0}$ sufficiently close to $X_{*}$. If $\mathrm{F}=\mathrm{L}(X,\mathcal{R}(X))$, $\lambda=\mathcal{R}(X)$ arises from an RQI as in section 1.1, we can use $\mathrm{G}_{\mathrm{L}}$ in place of $\mathrm{G}$:

[TABLE]

The above is a constrained version of the Chebyshev iteration, see section 3. For the eigenvalue problems, expanding $\mathfrak{r}(X,\eta)=\frac{X+\eta}{(1+|\eta|^{2})^{1/2}}$ to Taylor series in $\eta$ to the second term gives $\mathfrak{r}_{\eta\eta}(X,0;\eta^{[2]})=-X|\eta|^{2}$. In section 6, we use this to verify cubic convergence for several eigenvalue RQIs. The term $\mathrm{F}^{\prime}(X)\mathfrak{r}_{\eta\eta}(X,0;\eta^{[2]})$ often vanishes at a solution $X_{*}$ for homogeneous tensor/eigenvalue problems. In section 5, we relate it to a connection constructed from $\mathfrak{r}$, and give an interpretation of the Chebyshev term as a second covariant derivative.

We apply the general analysis here to construct new RQIs for the tensor eigenpair problem in section 6.1. A unitary version of the Schur form RQI over $\mathbb{C}$ finds all complex eigenpairs and also identifies the real pairs as a by-product in the generic case, improving on previous work [11]. We find new complex pairs for the Motzkin polynomial [12, example 5.9], completing the eigenpair count for this classical example. We verify cubic convergence for several eigenvalue problems, propose a cubic convergence iteration for the generalized eigenvalue problem, and derive Chebyshev iterations in a few interesting cases. A natural question is if we can develop a framework for constrained homotopy continuation [13], using RQI-type iterations. We hope this work will lead to further research along this line.

1.4. Notations and outline

We mimic the convention of [14] for derivatives. We work with a base field $\mathbb{F}$ which is $\mathbb{R}$ or $\mathbb{C}$. We denote by $\mathbb{F}^{n\times m}$ the space of $n\times m$ matrices on $\mathbb{F}$, by $\mathtt{Lin}(V,W)$ the space of linear map between vector spaces $V$ and $W$. The zero vector in $V$ is denoted $0_{V}$, or just [math]. The inner product of two (column) vectors $X_{1},X_{2}$ is $X_{1}^{\mathsf{T}}X_{2}$ or $X_{1}^{*}X_{2}$, identifying $\mathbb{F}^{1\times 1}$ with $\mathbb{F}$.

For a vector-valued function $\mathrm{F}$ from an open subset $D\subset V$ to $W$, by $\mathrm{F}^{\prime}$ we denote the Jacobian or Fréchet derivative, and the directional derivative at $X\in D$ in direction $\eta$ is written $\mathrm{F}^{\prime}(X)\eta$ or $\mathrm{F}^{\prime}(X;\eta)$, the latter is preferable when we need grouping. Thus, if $D$ is an open subset in $V$, $\mathrm{F}^{\prime}$ is an operator-valued function, $\mathrm{F}^{\prime}(X)\in\mathtt{Lin}(V,W)$. Higher (partial) derivatives could be considered as a map from $D$ to the space of multilinear maps from $V$ to $W$. Thus, at $X\in D$, the $l$-th-order derivative of $\mathrm{F}$, denoted by $\mathrm{F}^{(l)}(X)$ is a $l$-linear map from $V$ to $W$. For $\eta_{1},\cdots,\eta_{l}\in V$, we denote its evaluation as $\mathrm{F}^{(l)}(X;\eta_{1},\cdots\eta_{l})$ or $\mathrm{F}^{(l)}(X;[\eta_{1},\cdots\eta_{l}])$. In the expressions $\mathrm{F}^{(l)}(X)\eta^{[l]}=\mathrm{F}^{(l)}(X;\eta^{[l]})$, $\eta^{[l]}$ denotes $\eta$ being repeated $l$ times. If $\Phi$ is an operator value function from $V$ to $\mathtt{Lin}(W_{1},W_{2})$, with $V,W_{1},W_{2}$ are vector spaces and $\Phi(X)\omega$ denotes its valuation at $X\in V$ operating on $\omega\in W_{1}$, then $\Phi^{(l)}(X;\eta_{1},\cdots\eta_{l})\omega$ denotes the $l$-th order derivative of $\Phi$ evaluated in direction $\eta_{1},\cdots\eta_{l}$ and operate at $\omega$. In general, we use the round brackets for base variable evaluations or for groupings of base and directional variables. The semicolumn separates the base variable(s) from the directional derivative variables. The $l$-terms multi-dimensional Taylor series expansion around $X_{0}$ is [14, NR 3.3-3]

[TABLE]

Partial derivatives are sometimes denoted by the usual subscript convention, eg $\mathrm{L}_{X},\mathrm{L_{\lambda}}$ when convenient. However, we will use a position-based notation for more complex formulas. If $\mathrm{F}$ is defined using two vector variables $X,Y\in D_{1}\times D_{2}$ for two sets $D_{1},D_{2}$, then we write $\mathrm{F}(X,.)$ for the map $Y\mapsto\mathrm{F}(X,Y)$ with fixed $X$ in the appropriate domain, and write $\mathrm{F}(.,Y)$ similarly. The partial derivatives in $Y$ evaluated at $Y$ are $\mathrm{F}(X,.)^{\prime}(Y)=\mathrm{F}_{Y}(X,Y)$, which is a linear operator, and $\mathrm{F}(X,.)^{\prime}(Y;\eta)$ denotes its valuation at $\eta$. Thus, for a function $\mathrm{F}$ and a retraction $\mathfrak{r}$, the directional derivative in direction $(0,\xi)$ of $\mathrm{F}(\mathfrak{r}(X,\eta))$ is $\mathrm{F}^{\prime}(\mathfrak{r}(X,\eta);\mathfrak{r}(X,.)^{\prime}(\eta;\xi))$ by the chain rule. The notation $\mathfrak{r}(X,.)^{\prime}(\eta)\xi=\mathfrak{r}(X,.)^{\prime}(\eta;\xi)$ is also used.

Denote by $A_{V_{0}\downarrow W_{0}}$ the restriction of a map $A\in\mathtt{Lin}(V,W)$ from two spaces $V$ and $W$ to subspaces $V_{0}\subset V,W_{0}\subset W$, with $A(V_{0})\subset W_{0}$. The inverse, if exists, is denoted by $A_{V_{0}\downarrow W_{0}}^{-1}$, shorthand for $(A_{V_{0}\downarrow W_{0}})^{-1}$.

To focus on the main ideas, we will assume sufficient smoothness, and use the word smooth for short. Most results hold for class $C^{3}$ or $C^{4}$.

The main idea of the convergence analysis is to require the retraction $\mathfrak{r}$ maps an open ball $B_{X_{f}}(\rho)$ of the tangent space of a feasible point $X_{f}$ (an initial or final point) to an open subset of $\mathcal{M}$ in section 2, then translate the problem to one of linear constraint, with a nontrivial retraction $\mathfrak{s}$. The linear constraint case in section 3 is the most technical, the proofs are in the appendices. The iterations in section 4 are as described in the introduction, with applications in section 6. The relationship with the second covariant derivative is in section 5.

2. Feasibility perturbation and retractions

We state the definition and prove the basic properties of a retraction in this section, illustrated with the example of the rescaling retraction of the sphere. We also define the rescaling retractions and compute their Taylor expansion in the tangent variable.

2.1. Definition and basic properties

Assume $\mathrm{C}:\mathcal{E}\to\mathcal{E}_{L}$ is a smooth map, with Jacobian $\mathrm{C}^{\prime}(X)$ which is onto $\mathcal{E}_{L}$ at any $X\in\mathcal{E}$ with $\mathrm{C}(X)=0$. Let $\mathcal{M}$ be the corresponding feasible set, the solution set of $\mathrm{C}(X)=0$. This is a multidimensional smooth surface, a manifold. Denote by $T\mathcal{M}$, its tangent bundle, the subset of $\mathcal{E}^{2}=\mathcal{E}\times\mathcal{E}$ of pairs $(X,\eta)$ such that

[TABLE]

For a fixed $X\in\mathcal{M}$, the vector space defined by $\mathrm{C}_{X}(X)\eta=0$ is the tangent space at $X$, denoted by $T_{X}\mathcal{M}$, or $T_{X}$ for short subsequently, and $\eta$ is called a tangent vector. In elementary physics, an element $\eta\in T_{X}\mathcal{M}$ is considered as a velocity vector of a particle moving smoothly on $\mathcal{M}$. The Jacobian of the constraints eq. 2.1 is a block diagonal matrix with two diagonal blocks $\mathrm{C}^{\prime}(X)$, hence is onto $\mathcal{E}_{L}^{2}$, thus $T\mathcal{M}$ is also a smooth manifold. In the next two paragraphs, we summarize the basic geometric facts required, well-known for surfaces and curves.

Locally, $\mathcal{M}$ is parametrized by open subsets of $E_{T}:=\mathbb{R}^{\dim\mathcal{E}-\dim\mathcal{E}_{L}}$. For $X\in\mathcal{M}$, a local parametrization around $X$ is defined as a smooth map $\phi$, injective, from an open subset $\Omega$ of $E_{T}$ into $\mathcal{E}$, with $X\in\phi(\Omega)\subset\mathcal{M}$, such that $\phi^{-1}$, defined on $\phi(\Omega)$, is continuous, and $\phi^{\prime}(z)$ is injective for $z\in\Omega$. The full rank assumption of $\mathrm{C}^{\prime}(X)$ implies there is a local parametrization around $X$.

Differentiating $\mathrm{C}(\phi(z))=0$ in direction $v\in\mathcal{E}_{T}$, we get $\mathrm{C}^{\prime}(\phi(z))\phi^{\prime}(z)v=0$, or $\phi^{\prime}(z)v$ is tangent to $\mathcal{M}$ at $\phi(z)$. The dimension count and the injective assumption of $\phi^{\prime}(z)$ imply $\phi^{\prime}(z)$ is bijective to $T_{\phi(z)}$. Let $\Omega_{2}$ be an open subset of $\mathcal{E}_{T}$, then for $(z,v)\in\Omega\times\Omega_{2}\subset\Omega\times\mathcal{E}_{T}$, define

[TABLE]

then $D_{\phi}$ is a local parametrization of $T\mathcal{M}\subset\mathcal{E}^{2}$.

A map $\mathfrak{r}$ from an open subset $T_{1}$ of $T\mathcal{M}$ to $\mathcal{M}$ is smooth at $(X,\eta)\in T_{1}$ if for a local parametrization $D_{\phi}$ on $\Omega\times\Omega_{2}$ around $(X,\eta)$ and $\phi_{1}$ on $\Omega_{1}$ around $\mathfrak{r}(X,\eta)$, the combined map $\phi^{-1}_{1}\circ\mathfrak{r}\circ D_{\phi}$ from $D_{\phi}^{-1}(T_{1})\cap\Omega\times\Omega_{2}$ to $\Omega_{1}$ is smooth. This is independent of parametrizations.

We are interested in a perturbation $\mathfrak{r}$ mapping a pair $(X,\eta)\in T\mathcal{M}$ to $\mathcal{M}$ for $\eta$ sufficiently small, with $\mathfrak{r}(X,0)=X$, since our Newton increments are tangent vectors.

Lemma 2.1**.**

Fix $X\in\mathcal{M}$. If $\mathtt{K}$ is a $C^{1}$ map from a neighborhood $\mathcal{U}\subset T_{X}$ of $\xi\in T_{X}$ to $\mathcal{M}\subset\mathcal{E}$ then for $\eta\in T_{X}$, $\mathtt{K}^{\prime}(\xi)$ as a map from $T_{X}$ to $\mathcal{E}$ satisfies $\mathtt{K}^{\prime}(\xi)\eta\in T_{\mathtt{K}(\xi)}$, or $\mathtt{K}^{\prime}(\xi)(\mathcal{U})\subset T_{\mathtt{K}(\xi)}$.

Proof.

Differentiate $\mathrm{C}(\mathtt{K}(\xi+t\eta))=0$ in $t$ at $t=0$, we get $\mathrm{C}^{\prime}(\mathtt{K}(\xi))\mathtt{K}^{\prime}(\xi;\eta)=0$ by the chain rule. ∎

Applying this to feasibility perturbations, for $\xi\in T_{X}$, then $\mathfrak{r}(X,.)^{\prime}(\xi)$ could be considered as a map from $T_{X}$ to $T_{\mathfrak{r}(X,\xi)}$. In particular, $\mathfrak{r}(X,.)^{\prime}(0)$ maps $T_{X}$ to itself. Note $\mathfrak{r}(X,.)(\xi)^{\prime}\in\mathtt{Lin}(T_{X},\mathcal{E})$, but we will denote $\mathfrak{r}(X,.)^{\prime}(\xi)^{-1}$ the inverse, if exists, of $\mathfrak{r}(X,.)^{\prime}(\xi)_{T_{X}\downarrow T_{\mathfrak{r}(X,\xi)}}$ (recall this means restricting the range of $\mathfrak{r}(X,.)^{\prime}(\xi)$ to the image).

Remark 2.1*.*

Technically, it is simplest to consider a retraction as a map from $T\mathcal{M}$ onto $\mathcal{M}$ satisfying certain requirements on the Taylor series. This works for the rescaling map on the sphere, as $(X+\eta)/\lvert X+\eta\rvert$ is defined for any $\eta\in T_{X}$. However, we do want to consider, for example, rescaling maps for other constraints where rescaling only works for tangent vectors close enough to $0_{X}$. The following makes explicit the requirement that the radius where a perturbation exists is not too small, required in iterations, otherwise, it is adapted from [5, 15].

Definition 2.2**.**

For a point $X_{f}\in\mathcal{M}$, a retraction around $X_{f}$ (or a retraction to $\Omega^{\mathfrak{r}}$) is a smooth map $\mathfrak{r}$ from an open subset $\mathtt{T}_{\mathfrak{r}}\subset T\mathcal{M}$ to $\mathcal{M}$, with range in $\mathcal{M}$ containing a neighborhood $\Omega^{\mathfrak{r}}$ of $X_{f}$, and there is $\rho_{\mathfrak{r}}>0$ such that $\mathtt{T}_{\mathfrak{r}}$ contains

[TABLE]

If $\mathfrak{r}(X,0)$ is defined then we require $\mathfrak{r}(X,0)=X$ and $\mathfrak{r}(X,.)^{\prime}(0)=I_{T_{X}}$.

Remark 2.3*.*

Let $\phi$ be a local parametrization of $\mathcal{M}$ with $\phi(0)=X_{f}$. Requiring the existence of $\Omega^{\mathfrak{r}}$ and $\rho_{\mathfrak{r}}$ such that $\mathfrak{r}$ is defined in $\mathtt{T}\Omega^{\mathfrak{r}}_{\rho_{\mathfrak{r}}}$ is equivalent to requiring the range of $\mathfrak{r}$ containing $\{(\phi(z)|\;|z|\leq\rho_{1}\}$ and the domain of $\mathfrak{r}$ contains a set $\mathtt{T}\Omega^{\mathfrak{r}}_{\phi,\rho_{1},\rho_{2}}$ for some $\rho_{1},\rho_{2}>0$ where

[TABLE]

This follows as $\phi^{\prime}(z)$ and $\phi^{\prime}(z)^{-1}$ exist and are continuous for a parametrization, thus bounded uniformly on a bounded set, so we can translate a bound on $|\eta|=|\phi^{\prime}(z)v|$ to a bound on $v$, and vice versa.

We can normalize a feasibility perturbation to a retraction. If $\mathfrak{r}$ is feasibility perturbation satisfying the conditions of the retraction except for $\mathfrak{r}(X,.)^{\prime}(0)$ is only assumed to be invertible, then $\hat{\mathfrak{r}}(X,\eta):=\mathfrak{r}(X,(\mathfrak{r}(X,.)^{\prime}(0))^{-1}\eta)$ is a retraction since $\hat{\mathfrak{r}}(X,.)^{\prime}(0;\eta)=\mathfrak{r}(X,.)^{\prime}(0;(\mathfrak{r}(X,.)^{\prime}(0))^{-1}\eta)=\eta.$

If $\mathcal{M}$ is the unit sphere, for a feasible point $X_{f}$, the rescaling of $X_{f}+\eta$ for $\eta\in T_{X_{f}}$ maps $T_{X_{f}}$ bijectively to the hemisphere centered at $X_{f}$, and this map parametrizes the hemisphere.

We now generalize this result, showing $\mathfrak{r}(X_{f},.)$ parametrizes an open set around $X_{f}$.

Proposition 2.1**.**

Assume $\mathfrak{r}$ is a retraction around a feasible point $X_{f}\in\mathcal{M}$.

1) There is a radius $\rho>0$ such that the map $\mathtt{K}_{f}:=\mathfrak{r}(X_{f},.):\eta\mapsto\mathfrak{r}(X_{f},\eta)$ from the ball $B_{X_{f}}(\rho):=\{\eta\in T_{X_{f}}||\eta|<\rho\}$ to its image $C_{X_{f}}(\rho)\subset\mathcal{M}$ is a local parametrization. The ball $B_{X_{f}}(\rho)\subset T_{X_{f}}$ is called a retraction ball, $\mathtt{K}_{f}(B_{X_{f}}(\rho))=C_{X_{f}}(\rho)\subset\mathcal{M}$ is called the retraction cap.

2) Thus, the set $\mathtt{T}\Omega^{\mathfrak{r}}_{\mathtt{K}_{f},\rho_{1},\rho_{2}}$ in eq. 2.4 exists, for each $X_{f}\in\Omega$, there are radii $\rho_{1}>0$, $\rho_{2}>0$ such that $\mathtt{K}_{f}$ maps $B_{X_{f}}(\rho_{1})$ to $C_{X_{f}}(\rho_{1})$ with invertible Jacobian, and for all $X\in C_{X_{f}}(\rho_{1})$, and $\delta\in T_{X_{f}}$ with $|\delta|<\rho_{2}$, if $X=\mathfrak{r}(X_{f},\xi)$ then $\mathfrak{r}(X,\mathtt{K}_{f}^{\prime}(\xi;\delta))$ and $\mathtt{K}_{f}^{-1}\mathfrak{r}(X,\mathtt{K}_{f}^{\prime}(\xi;\delta))$ exist.

Proof.

Let $\phi$ be a local parametrization of $\mathcal{M}$ around $X_{f}$. By remark 2.3, there are radii $r_{0},\Delta_{0}$ such that $\phi(B_{0})$ is in the range of $\mathfrak{r}$, where $B_{0}$ is the ball radius $r_{0}$ in $\mathcal{E}_{T}=\mathbb{R}^{\dim\mathcal{E}-\dim\mathcal{E}_{L}}$, and $\mathtt{T}\Omega^{\mathfrak{r}}_{\phi,r_{0},\Delta_{0}}$ exists.

Since $\mathtt{K}_{f}:=\mathfrak{r}(X_{f},.)$ is continuous, $\mathtt{K}_{f}^{-1}(\phi(B_{0}))$ is open in $T_{X_{f}}$, thus there is a neighborhood $\mathcal{U}\subset T_{X_{f}}$ of $0_{T_{X_{f}}}$ such that $\mathtt{K}_{f}(\mathcal{U})\subset\phi(B_{0})$. Consider the map $h:\mathcal{U}\times\mathcal{E}_{T}\to\mathcal{E}_{T}$, $h(\eta,z)=\phi^{-1}(\mathtt{K}_{f}(\eta))-z$ for $(\eta,z)\in\mathcal{U}\times\mathcal{E}_{T}$. Since $\mathtt{K}_{f}^{\prime}(0)\xi=\xi$, by the chain rule, the partial derivative $h_{\eta}(0,0)$ is $\xi\mapsto\phi^{\prime}(0)^{-1}\xi\in T_{X_{f}}$ for $\xi\in T_{X_{f}}$, where $\phi^{\prime}(0)$ is invertible as a map to $T_{X_{f}}$. Thus, the implicit function theorem [14, 5.2.4] applies, we have open balls $B_{1}\subset B_{0}$ centered at $0_{\mathcal{E}_{T}}$, $B_{\eta}\subset T_{X_{f}}$ at $0_{T_{X_{f}}}$, and a unique function $\eta_{f}(z)$ with $h(\eta_{f}(z),z)=0$ for $z\in B_{0},\eta_{f}(z)\in B_{\eta}$.

We have $\phi^{-1}(\mathtt{K}_{f}(\eta_{f}(z))=z$ or $\mathtt{K}_{f}(\eta_{f}(z))=\phi(z)$. If $Y\in\phi(B_{1})$, then $\eta=\eta_{f}\circ\phi^{-1}(Y)$ satisfies $\mathtt{K}_{f}(\eta)=\phi(\phi^{-1}(Y))=Y$. Since $\mathtt{K}_{f}$ is continuous, $\mathtt{K}_{f}^{-1}(\phi(B_{1}))$ is open, and its intersection with $B_{1}$ is open, containing a ball of radius $\rho$. Restricting to $B_{X_{f}}(\rho)$, the uniqueness of the implicit function shows $\mathtt{K}_{f}$ is one-to-one, with continuous inverse $\eta_{f}\circ\phi^{-1}$. The remaining properties of parametrizations are verified. This proves 1).

Item 2) follows from remark 2.3. ∎

We offer the first clue why the (partial) Hessian of $\mathfrak{r}$ appears in claim 2.

Proposition 2.2**.**

If $A$ is a smooth function on $\mathcal{E}$ with values in a vector space $V$, assume $X\in\mathcal{M}$ and $\eta\in T_{X}\mathcal{M}$ then $A^{\prime}(X)\eta$ is only dependent on the values of $A$ on $\mathcal{M}$. If $\mathfrak{r}$ is a retraction then

[TABLE]

is also only dependent on the values of $A$ on $\mathcal{M}$.

Proof.

The first statement is well-known. Let $\phi$ be a local parametrization around $X$ defined on $\Omega\subset\mathcal{E}_{T}$, and $\eta=\phi^{\prime}(0)v$ for some $v\in\mathcal{E}_{T}$. Consider the function $g(t)=A(\phi(tv))$ for $t$ sufficiently small. Then $\dot{g}(0)=A^{\prime}(\phi(0))\phi^{\prime}(0)v=A^{\prime}(X)\eta$, but $g$ depends on values of $A$ in $\phi(\Omega)\subset\mathcal{M}$ only.

If $\mathfrak{r}$ is a retraction, consider $f(t)=A(\mathfrak{r}(X,t\eta))$. The second derivative at [math] of $f$ is

[TABLE]

by the chain rule. On the other hand, $A(\mathfrak{r}(X,t\eta))$ is dependent on the values of $A$ on $\mathcal{M}$ only. ∎

We can get similar statements for higher derivatives. The function $\mathfrak{s}(\xi,\delta)$ below expresses a retraction to $X=\mathfrak{r}(X_{f},\xi)\in\mathcal{M},\xi\in T_{X_{f}}$ in terms of a retraction to a chosen point $X_{f}\in\mathcal{M}$. It appears when we translate an iteration on $\mathcal{M}$ to an iteration on $T_{X_{f}}$ by a change of variable.

Proposition 2.3**.**

For $X_{f}\in\subset\mathcal{M}$, set $\mathtt{K}_{f}:=\mathfrak{r}(X_{f},.)$, then with $\rho_{1},\rho_{2}>0$ as in 2) proposition 2.1,

[TABLE]

exists for $|\xi|<\rho_{1}$ and $|\delta|<\rho_{2}$ and is a retraction to $T_{X_{f}}\in\mathcal{E}$. Moreover,

[TABLE]

Proof.

Recall for $\phi_{1},\phi_{2}\in T_{X}$, $\mathfrak{r}(X,.)^{\prime}(\phi_{1};\phi_{2})$ is $\lim_{t\to 0}\frac{1}{t}(\mathfrak{r}(X,\phi_{1}+t\phi_{2})-\mathfrak{r}(X,\phi_{1}))$. Rewrite the equation for $\mathfrak{s}$ then differentiate it with respect to $\delta$ in direction $\epsilon\in T_{X_{f}}$ twice

[TABLE]

Set $\delta$ to zero, we get the first three terms of the power series expansion of $\mathfrak{s}(\xi,\epsilon)$,

[TABLE]

where we used $\mathfrak{r}(\mathtt{K}_{f},\xi),.)^{\prime}(0;\mathtt{K}_{f}^{\prime}(\xi;\epsilon))=\mathtt{K}_{f}^{\prime}(\xi;\epsilon)$ as $\mathfrak{r}$ is a retraction, note $(\mathtt{K}_{f}^{\prime}(\xi))^{-1}$ exists in the retraction ball. The statement follows from the Taylor series expansion. ∎

2.2. The rescaling retraction on the unit sphere

The computation of $\mathtt{K}_{f}^{\prime}$ for the unit sphere is known in [6, section 8.1].

Proposition 2.4**.**

Consider the sphere $\mathcal{M}=S^{n-1}\in\mathbb{R}^{n}$ ($n\geq 2$), defined by the equation $X^{\mathsf{T}}X=1$. For $X_{f}\in S^{n-1}$, the retraction $\mathtt{K}_{f}(\xi)=\mathfrak{r}(X_{f},\xi)=\frac{1}{(1+|\xi|^{2})^{1/2}}(X_{f}+\xi)$ defined on $T_{\mathfrak{r}}=TS^{n-1}$ maps $\xi\in T_{X_{f}}$ to the hemisphere $\{X\in S^{n-1}|\;X_{f}^{\mathsf{T}}X>0\}$. For $\rho>0$, let $C_{X_{f}}(\rho)$ be the retraction cap, then

[TABLE]

Any $\rho_{1},\rho_{2}>0$ satisfying condition 2.15 below satisfies the first condition of 2), proposition 2.1, that is $\mathfrak{s}(X_{f};\xi,\delta)$ exists if $|\xi|<\rho_{1},|\delta|<\rho_{2}$.

[TABLE]

Proof.

The Taylor series expansion of $(1+|\xi|^{2})^{-1/2}$ gives eq. 2.10. Since $X_{f}^{\mathsf{T}}\mathfrak{r}(X_{f},\xi)=\frac{1}{(1+|\xi|^{2})^{1/2}}>0$, and the function $(1+\rho^{2})^{-1/2}$ is decreasing in $\rho$, the image of $\mathtt{K}_{f}$ is in the hemisphere $X_{f}^{T}X>0$ and (2.11) follows. We can verify directly $\mathtt{K}_{f}^{-1}X$ on the right-hand side of eq. 2.12 satisfying $X_{f}^{\mathsf{T}}(\mathtt{K}_{f}^{-1}X)=0$, and $X_{f}+\mathtt{K}_{f}^{-1}X=\frac{1}{X_{f}^{\mathsf{T}}X}X$ is proportional to $X$, thus $\frac{1}{|X_{f}+\mathtt{K}_{f}^{-1}X|}(X_{f}+\mathtt{K}_{f}^{-1}X)=X$.

A routine derivative gives $\mathtt{K}_{f}^{\prime}$. Since $\mathfrak{r}(\mathtt{K}_{f}(\xi),\mathtt{K}_{f}^{\prime}(\xi,\delta))$ is proportional to $\mathtt{K}_{f}(\xi)+\mathtt{K}_{f}^{\prime}(\xi,\delta)$, we assume the scaling factor is of the form $c(1+|\xi|^{2})^{1/2}$, then using $X_{f}^{\mathsf{T}}\xi=X_{f}^{\mathsf{T}}\delta=0$

[TABLE]

which gives us the formula for $\mathfrak{s}$. It exists if $f_{1}=1+|\xi|^{2}-\xi^{\mathsf{T}}\delta>0$. If $|\delta|=\rho$ is constant, then $f_{1}$ is smallest if $\delta=\frac{\rho}{|\xi|}\xi$, with $f_{1}=1+|\xi|^{2}-\rho|\xi|$. The region defined by $1+\rho_{1}^{2}-\rho_{1}\rho_{2}>0,\rho_{1}>0,\rho_{2}>0$ could be divided into three subregions, the union of two is described in eq. 2.15, characterized by the condition that if $0<z_{1}<\rho_{1},0<z_{2}<\rho_{2}$ then $1+z_{1}^{2}-z_{1}z_{2}>0$. The remaining subregion with $\rho_{1}\geq 1,\rho_{2}\geq 2$ does not have this property to ensure $\mathfrak{s}(\xi,\delta)$ exists, if $|\xi|=z_{1},|\delta|=\rho=z_{2}$. ∎

2.3. Rescaling retractions

We generalize the rescaling retraction on the sphere to feasible sets with a single constraint.

Proposition 2.5** (Rescaling retraction).**

Assume $\mathcal{E}_{L}=\mathbb{R}$, and $\mathrm{C}^{\prime}(X;X)\neq 0\in\mathcal{E}_{L}=\mathbb{R}$ near $X_{f}\in\mathcal{M}$. There is a neighborhood $\Omega^{\mathfrak{r}}\subset\mathcal{M}$ of $X_{f}$ and $\rho_{\mathfrak{r}}>0$ such that the equation in $\gamma\in\mathbb{R}$

[TABLE]

for $(X,\eta)\in\mathtt{T}\Omega^{\mathfrak{r}}_{\rho_{\mathfrak{r}}}\subset T\mathcal{M}$ (in eq. 2.3) has an implicit function solution with $\gamma=1$ at $\eta=0$. Then $\mathfrak{r}(X,\eta)=\gamma(X+\eta)$ is a retraction defined on $\mathtt{T}\Omega^{\mathfrak{r}}$. We have the Taylor expansion

[TABLE]

Proof.

Let $\phi$ be a local parametrization from $\Omega_{T}\subset\mathcal{E}_{T}=\mathbb{R}^{\dim(\mathcal{E})-\dim(\mathcal{E}_{L})}$ of $\mathcal{M}$ near $X_{f}$. Consider $h(z,v,\gamma)=\mathrm{C}(\gamma(\phi(z)+\phi^{\prime}(z)v))\in\mathcal{E}_{L}$ defined on $(z,v,\gamma)\in\Omega_{T}\times E_{T}\times\mathbb{R}$. Then $h(0,0,1)=\mathrm{C}(X_{f})=0$, $h_{\gamma}(0,0,1)=\mathrm{C}^{\prime}(X_{f};X_{f})\neq 0$. Thus, the implicit function theorem applies, giving us a ball $D=D_{1}\times D_{2}$ in $\Omega_{T}\times\mathcal{E}_{T}$ of points $\{(z,v)\}$ with $|z|<\rho_{1},|v|<\rho_{2}$ and $\delta_{\gamma}>0$ such that $\gamma$ can be solved uniquely as a function $\gamma_{\phi}(z,v)$ with $\gamma_{\phi}(0,0)=1$ for $(z,v)\in D$ and $|\gamma-1|<\delta_{\gamma}$. The conditions of remark 2.3 for $\mathfrak{r}(X,\eta)=\gamma(X+\eta)$ follows by properties of the implicit function. By local uniqueness of the implicit solution, $\gamma_{\phi}(z,0)=1$, or $\mathfrak{r}(X,0)=X$ for $X\in\phi(D_{1})$.

Composing with $\phi^{-1}$, $\gamma=\gamma_{\phi}\circ\phi^{-1}$ is a function of $X$ and $\eta$, its partial derivatives are evaluated by implicit function rule for $g(\gamma,X,\eta)=\mathrm{C}(\gamma(X+\eta))$. Hence, we can compute $\mathfrak{r}(X,.)^{\prime}$ using

[TABLE]

where $\xi\in T_{X}$, hence $\mathrm{C}^{\prime}(X,\xi)=0$. Consider a fixed $X$, set $\gamma_{\eta}=\gamma(X,.)^{\prime},\gamma_{\eta\eta}=\gamma(X,.)^{(2)}$. Take directional derivative of $\mathrm{C}(\gamma(X+\eta))=0$ in $\eta$ in tangent directions $\xi_{1},\xi_{2}$ consecutively

[TABLE]

At $\eta=0,\gamma=1$, $\gamma_{\eta}(0;\xi)=0$ for $\xi\in T_{X}\mathcal{M}$ since $\mathrm{C}^{\prime}(X,\xi)=0$ as before, the above gives us

[TABLE]

Thus, equation (2.17) follows from $\mathfrak{r}(X,.)^{(2)}(0;\eta^{[2]})=\gamma_{\eta\eta}(0;\eta^{[2]})X$ and the Taylor formula. ∎

The requirement $\mathrm{C}^{\prime}(X;X)\neq 0$ means the vector $X$ is not tangent to $\mathcal{M}$. Note that when $\mathrm{C}(X)=f(X)-1$ and $f$ is a homogeneous function of degree $n\neq 0$, Euler’s identity implies $\mathrm{C}^{\prime}(X;X)=nf(X)=n\neq 0$, the equation $\mathrm{C}(\gamma(X+\eta))=0$ has a simple solution $\gamma=f(X+\eta)^{-1/n}$ if $\eta$ is sufficiently small. The sphere corresponds to $f(X)=X^{\mathsf{T}}X$ and $n=2$. In general, $\gamma$ is a solution of a scalar equation.

3. Linear constraints and unconstrained case with a projection

This section is the most technical in the paper. The main proofs are deferred to the appendix.

Using the result of section 2, we translate the constrained problem to one of linear constraints, by considering a reference point $X_{f}$, and solve $\mathcal{F}(\xi)=\mathrm{F}(\mathfrak{r}(X_{f},\xi))=0$ instead of $\mathrm{F}(X)=0$, for $\xi\in T_{X_{f}}$ with the projection $\Phi(\xi)=\Pi(\mathfrak{r}(X_{f},\xi))$, $\Phi(\xi)\mathcal{F}(\xi)=\mathcal{F}(\xi)$, and the retraction $\mathfrak{s}$ in proposition 2.3. The feasible set is now called $T_{0}$ ($=T_{X_{f}}$). The variable becomes $\xi$ instead of $X$.

A local linear constrained problem becomes an unconstrained problem if we parametrized the constraint set. Hence, we can consider the unconstrained problem over a vector space $T_{0}$.

Thus, we will start with two vector spaces $T_{0}$ and $\mathcal{E}$, with a function $\mathcal{F}$ from an open subset $\Omega$ from $T_{0}$ to $\mathcal{E}$, with Jacobian $\mathcal{F}^{\prime}$. We assume $\Omega$ to be convex, and $\dim T_{0}\leq\dim\mathcal{E}$.

When $\dim T_{0}=\dim\mathcal{E}$, for $\xi_{*}\in\Omega$, if $\mathcal{F}(\xi_{*})=0$ with $\mathcal{F}^{\prime}(\xi_{*})$ invertible, an inverse function $\mathcal{F}_{\mathfrak{I}}$ exists, $\mathcal{F}_{\mathfrak{I}}(\mathcal{F}(\xi))=\xi$ near $\xi_{*}$. For $Y\in\mathcal{U}\subset\mathcal{E}$ where the inverse function is defined, the $k$ terms Taylor expansion of $\mathcal{F}_{\mathfrak{I}}$ at $\mathcal{F}(\xi_{i})$ (assumed to be in $\mathcal{U}$) is

[TABLE]

The Taylor coefficients for $\mathcal{F}_{\mathfrak{I}}$ could be computed by the implicit-valued function theorem, or by power series substitution. We have $\mathcal{F}_{\mathfrak{I}}^{\prime}(Y)\eta=(\mathcal{F}^{\prime}(\xi))^{-1}\eta$ with $\xi=\mathcal{F}_{\mathfrak{I}}(Y)$ and

[TABLE]

Following [4], we can construct an iteration with order-$k$ convergence as follows. For $Y=0$, $\mathcal{F}_{\mathfrak{I}}(0)=\xi_{*}$, and if we set $\xi_{i+1}=\xi_{i}+\sum_{j=1}^{k-1}\frac{1}{j!}\mathcal{F}_{\mathfrak{I}}^{(j)}(\mathcal{F}(\xi_{i});[-\mathcal{F}(\xi_{i})]^{[j]})$, the residual $\xi_{*}-\xi_{i+1}$, from the above expansion is $O(|\mathcal{F}(\xi_{i})|^{k})=O(|\xi_{i}-\xi_{*}|^{k})$ if we consider a bounded region and $\mathcal{F}$ has bounded derivatives. Thus, $k$-th order convergence follows from the Taylor series of $\mathcal{F}_{\mathfrak{I}}$ by design. The case $k=2$ corresponds to the Newton method with increment $\delta_{2}=-(\mathcal{F}^{\prime}(\xi))^{-1}\mathcal{F}(\xi)$, and the case $k=3$ corresponds to the Chebyshev method, with step $\xi^{\mathrm{C}}=\xi+\delta_{2}-\frac{1}{2}(\mathcal{F}^{\prime}(\xi))^{-1}\mathcal{F}^{(2)}(\xi;\delta_{2}^{[2]})$.

When $\dim T_{0}<\dim\mathcal{E}$ and $\mathcal{F}^{\prime}(\xi_{*})$ is not invertible, assume there is an affine projection $\Phi$ such that $\Phi(\xi)\mathcal{F}(\xi)=\mathcal{F}(\xi)$ for $\xi\in\Omega$. We also assume $\dim(\operatorname{Im}(\Phi(\xi)))=\dim(T_{0})$. In proposition 3.1, we define a left inverse $\mathcal{F}^{\natural}$ of $\mathcal{F}^{\prime}$ and compute derivatives of $\mathcal{F}_{\diamond}$, playing the role of $\mathcal{F}_{\mathfrak{I}}$ (compare eqs. 3.8 and 3.9 with $\mathcal{F}_{\mathfrak{I}}^{\prime}$ and $\mathcal{F}_{\mathfrak{I}}^{(2)}$). The proof is in section A.3.

Recall the notation $A_{V_{0}\downarrow W_{0}}$ denoting the restriction of a map $A$ from two spaces $V$ and $W$ to subspaces $V_{0}\subset V,W_{0}\subset W$, with inverse (if exists) denoted by $A_{V_{0}\downarrow W_{0}}^{-1}$, shorthand for $(A_{V_{0}\downarrow W_{0}})^{-1}$.

Proposition 3.1**.**

Assume the map $\mathcal{F}$ from $\Omega\subset T_{0}$ to $\mathcal{E}$ and the projection function $\Phi:\Omega\mapsto\mathtt{Lin}(\mathcal{E},\mathcal{E})$ are smooth. Assume for $\xi\in\Omega$

[TABLE]

Set $\mathcal{E}_{\xi}:=\operatorname{Im}(\Phi(\xi))\subset\mathcal{E}$. If $\Phi(\xi)\mathcal{F}^{\prime}(\xi)$ is a bijection from $T_{0}$ to $\mathcal{E}_{\xi}$, for $\omega\in\mathcal{E}$, let $\mathcal{F}^{\natural}(\xi)\omega\in T_{0}$ be the unique solution $\eta$ of $\Phi(\xi)\mathcal{F}^{\prime}(\xi)\eta=\Phi(\xi)\omega$, thus $\mathcal{F}^{\natural}(\xi)_{\mathcal{E}_{\xi}\downarrow T_{0}}=(\Phi(\xi)\mathcal{F}^{\prime}(\xi))_{T_{0}\downarrow\mathcal{E}_{\xi}}^{-1}$ and

[TABLE]

1) For each $\xi$, the linear map $\mathcal{F}^{\natural}(\xi)\in\mathtt{Lin}(\mathcal{E},T_{0})$, if exists, is a left inverse of $\mathcal{F}^{\prime}(\xi)$ and

[TABLE]

Assume $\mathcal{F}^{\natural}(\xi)$ is defined in a domain $D_{1}\subset\Omega$, then it is differentiable and for $\xi,\eta\in T_{0}$ and $\omega\in\mathcal{E}$

[TABLE]

2) Assume $\mathcal{F}^{\natural}(\xi)$ is defined in a domain $D_{1}\subset\Omega$. Then for $(\xi,Y)\in D_{1}\times\mathcal{E}$

[TABLE]

For $\xi_{f}\in\mathcal{E}$, there is an implicit-function solution $\xi=\mathcal{F}_{\diamond}(Y)$ for the second equation near $(\xi_{f},\mathcal{F}(\xi_{f}))\in D_{1}\times\mathcal{E}$, with $Y\in\mathcal{E}$, $\xi\in T_{0}$, where $\mathcal{F}_{\diamond}$ is of class $C^{k}$ if $\mathcal{F}$ and $\Phi$ are. If $(\xi,Y)$ satisfies eq. 3.6 then $(\xi,Z)$ also satisfies eq. 3.6 for $Z=\Phi(\xi)Y$. For $\omega\in\mathcal{E}$, we have

[TABLE]

Assume $\mathcal{F}$ and $\Phi$ are smooth, then for $\omega\in\mathcal{E}$

[TABLE]

3) If $(I_{T_{0}}-\mathcal{F}^{\natural}(\xi)\Phi^{\prime}(\xi;.;Y))^{-1}$ exists and is bounded in a region $D_{1}$, and $\mathcal{F}$ and $\Phi$ are of class $C^{k}$, $k\geq 1$ with bounded derivatives, then $\mathcal{F}_{\diamond}$ is also of class $C^{k}$ with bounded derivatives.

We now formulate the convergence theorem with a nontrivial retraction. We will see the retraction in proposition 2.3 will appear in RQIs. For the case of the sphere eq. 2.14, it has the form $\mathfrak{s}_{s}:(\xi,\delta)\mapsto\xi+\frac{1+|\xi|^{2}}{1+\xi^{\mathsf{T}}(\xi-\delta)}\delta$, for $\xi,\delta\in T_{0}$, and is only defined for $\delta$ satisfying eq. 2.15. Thus, we need to address retractions of this type.

Theorem 3.1**.**

Assume the function $\mathcal{F}$, the projection $\Phi$ are smooth satisfying $\Phi(\xi)\mathcal{F}(\xi)=\mathcal{F}(\xi)$ for $\xi\in\Omega\subset T_{0}$ as in eq. 3.1. 1) Assume $(\Phi(\xi)\mathcal{F}^{\prime}(\xi))_{T_{0}\downarrow\mathcal{E}_{\xi}}^{-1}$ is defined in a ball $B(\xi_{*},\rho)$ in $T_{0}$, with $\mathcal{F}(\xi_{*})=0$. For $\rho,\rho_{2}>0$, assume the map $\mathfrak{s}$ from $D_{1}=B(\xi_{*},\rho)\times B(0,\rho_{2})\subset T_{0}\times T_{0}$ to $T_{0}$ of class $C^{2}$ has the Taylor expansion in $\delta$

[TABLE]

such that there is a constant $c_{\mathfrak{s},2}$ with $|\mathfrak{s}_{2}(\xi,\delta)|\leq c_{2}|\delta|^{2}$ in $D_{1}$. Then there exists a radius $\rho_{3}$ such that if the starting point $\xi_{0}$ satisfies $|\xi_{0}-\xi_{*}|<\rho_{3}$ then the iteration defined by

[TABLE]

is well-defined and converges quadratically to $\xi_{*}$. 2) If $\mathfrak{s}$ is of class $C^{3}$ and the 3-term Taylor series of $\mathfrak{s}$ in $\delta$ is given by

[TABLE]

with $|\mathfrak{s}_{3}(\xi,\delta)|<c_{\mathfrak{s},3}|\delta|^{3}$ and $\mathcal{F}$ is of class $C^{3}$. If $\Phi(\xi_{*})\mathcal{G}(\xi_{*},\delta^{[2]})=0$ for all $\delta\in T_{0}$ where

[TABLE]

then the Newton iteration in 1) converges cubically to $\xi_{*}$. 3) With the smoothness assumption in 2), but we do not assume $\Phi(\xi_{*})\mathcal{G}(\xi_{*},\delta^{[2]})=0$, then there exists a radius $\rho_{4}$ such that for $|\xi_{0}-\xi_{*}|<\rho_{4}$, the following iteration converges cubically to $\xi_{*}$

[TABLE]

The proof is in section A.4.

4. Constrained equations and Rayleigh/Rayleigh-Chebyshev iterations

We now consider a feasible set $\mathcal{M}$ defined by a full range constraint $\mathrm{C}(X)=0$ with a retraction $\mathfrak{r}$. Let $\Pi$ be an affine projection function from $\mathcal{M}$ to $\mathtt{Lin}(\mathcal{E},\mathcal{E})$, $\Pi(X)^{2}=\Pi(X)$ for $X\in\mathcal{M}$. Denote $\mathcal{E}_{X}=\operatorname{Im}(\Pi(X)),T_{X}=\operatorname{Null}(\mathrm{C}^{\prime}(X))$. In the subset of $\mathcal{M}$ considered, we assume

[TABLE]

We try to solve $\mathrm{F}(X)=0$ for a function $\mathrm{F}$ from $\mathcal{M}$ to $\mathcal{E}$ satisfying the condition

[TABLE]

If $(\Pi(X)\mathrm{F}^{\prime}(X))_{T_{X}\downarrow E_{X}}$ is invertible, (i.e. the system in claim 1 is nondegenerate), define

[TABLE]

Proposition 4.1**.**

Fix a feasible point $X_{f}\in\mathcal{M}$, let $\mathfrak{r}$ be a retraction around $X_{f}$ with a retraction cap $\mathrm{C}_{X_{f}}(\rho)\subset\mathcal{M}$, parametrized by a retraction ball $B_{X_{f}}(\rho)\subset T_{X_{f}}$ via $\mathtt{K}_{f}:=\mathfrak{r}(X_{f},.)$. Let $\mathrm{F}$ be a map from $\mathcal{M}$ to $\mathcal{E}$, $\Pi$ be an affine projection. For $X=\mathfrak{r}(X_{f},\xi)\in C(X_{f},\rho)$, consider the Newton increment and step

[TABLE]

assume they are well-defined and $X^{\mathcal{N}}\in C_{X_{f}}(\rho)$. For $\xi\in B_{X_{f}}(\rho)$, set $\Phi(\xi):=\Pi(\mathfrak{r}(X_{f},\xi))$, $\mathcal{F}(\xi):=\mathrm{F}(\mathfrak{r}(X_{f},\xi))$. Then the Newton-type iteration $\xi^{\mathcal{N}_{f}}$ of the problem $\mathcal{F}(\xi)=0$

[TABLE]

with $\mathfrak{s}$ defined in eq. 2.5 corresponds to the Newton iteration $X^{\mathcal{N}}$. That means

[TABLE]

Proof.

From the chain rule, $\mathcal{F}^{\prime}(\xi)$ and its left-inverse $\mathcal{F}^{\natural}(\xi)$ in eq. 3.2 are

[TABLE]

where $(\mathfrak{r}(X_{f},.)^{\prime}(\xi))$ moves to the front in the inverse. For eqs. 4.8, 4.6 and 4.7, with $X=\mathfrak{r}(X_{f},\xi)$

[TABLE]

We now translate the theorems in section 3 to results on constrained iterations.

Theorem 4.1** (Local convergence).**

Assume $\mathrm{F}$ and $\Pi$ satisfy eqs. 4.1 and 4.2. Assume $\mathrm{F}(X_{*})=0$ and $(\Pi(X_{*})\mathrm{F}^{\prime}(X_{*}))_{T_{X_{*}}\downarrow E_{X_{*}}}$ is invertible. Let $\mathfrak{r}$ be a retraction around $X_{*}$. Then there is a radius $\rho>0$ defining a retraction ball $B_{X_{*}}(\rho)\subset T_{X_{*}}$ and cap $C_{X_{*}}(\rho)\subset\mathcal{M}$ such that for $X_{0}\in C_{X_{*}}(\rho)$, the iteration $X_{i+1}=X_{i}^{\mathcal{N}}$ defined in proposition 4.1 is well-defined and converges to $X_{*}$.

If $\Pi(X_{*})\mathrm{G}(X_{*};\eta^{[2]})=0$ for $\eta\in T_{X_{*}}$, where $\mathrm{G}$ is a tensor-valued function on $\mathcal{M}$ defined by

[TABLE]

then the Newton iteration $X_{i+1}=X_{i}^{\mathcal{N}}$ converges cubically. Define the Chebyshev increment

[TABLE]

Then there is $\rho_{C}>0$ such that for $X_{0}\in C_{X_{*}}(\rho_{C})$, the Rayleigh-Chebyshev iteration $X_{i+1}=\mathfrak{r}(X_{i},X_{i}^{\mathcal{N}_{C}})$ converge cubically to $X_{*}$.

Proof.

Using proposition 4.1, the Newton case follows from theorem 3.1. Set $\mathcal{F}(\xi)=\mathrm{F}(\mathfrak{r}(X_{f},\xi))$ for a feasible point $X_{f}$, then differentiate $\mathcal{F}^{\prime}(\xi)\delta=\mathrm{F}^{\prime}(\mathfrak{r}(X_{f},\xi))\mathfrak{r}(X_{f},.)^{\prime}(\xi;\delta)$ in direction $\delta$

[TABLE]

if $X=\mathfrak{r}(X_{f},\xi),\eta=\mathfrak{r}(X_{f},.)^{\prime}(\xi;\delta)$ and $\mathcal{G}$ defined in eq. 3.13. Cubic convergence of the Newton algorithm if $\mathrm{G}$ vanishes follows. The statement for Rayleigh-Chebyshev also follows. ∎

4.1. Rayleigh quotient iteration for constrained equations

As mentioned, for the equation eq. 1.1, if $\mathrm{L_{\lambda}}(X,\lambda)$ is onto, as a linear map from $\mathcal{E}_{L}$ to $\mathcal{E}$, then there is a left-inverse $\mathrm{L}^{-}_{\lambda}(X,\lambda)$ of $\mathrm{L_{\lambda}}(X,\lambda)$. If it is constructed smoothly, then with the Rayleigh quotient $\lambda=\mathcal{R}(X)$ solving $\mathrm{L}^{-}_{\lambda}(X,\lambda)\mathrm{L}(X,\lambda)=0$, the projection $\Pi(X):=I_{\mathcal{E}}-\mathrm{L}_{\lambda}(X,\mathcal{R}(X))\mathrm{L}_{\lambda}^{-}(X,\mathcal{R}(X))$ satisfies

[TABLE]

Let $\mathrm{F}(X)=\mathrm{L}(X,\mathcal{R}(X))$, then $\Pi(X)\mathrm{F}(X)=\mathrm{F}(X)$, thus, the results of the previous section apply. Using $\Pi(X)\mathrm{L_{\lambda}}(X,\mathcal{R}(X))=0$ to simplify $\Pi(X)\mathrm{F}^{\prime}(X)$ and $\Pi(X)\mathrm{G}(X)$, the following theorem follows from the chain rule. The second derivatives of $\mathrm{L}$ are bilinear maps, $\mathrm{L}_{XX}$ operates on two copies of $\mathcal{E}$ via the variable $\eta$ below, $\mathrm{L}_{X\lambda}$ operates on $\mathcal{E}\times\mathcal{E}_{L}$ via $\eta\in\mathcal{E}$ and $\mathcal{R}^{\prime}(X;\eta)\in\mathcal{E}_{L}$, and $\mathrm{L}_{\lambda\lambda}$ on $\mathcal{R}^{\prime}(X;\eta)^{[2]}\in\mathcal{E}_{L}^{2}$.

Theorem 4.2**.**

If the Rayleigh quotient $\lambda=\mathcal{R}(X)$ and the projection $\Pi$ satisfies eq. 4.13, with $\mathrm{L}$ and $\mathcal{R}$, $\mathrm{C}$, $\mathfrak{r}$, $\Pi$ are smooth, then with $\lambda=\mathcal{R}(X)$, the Newton increment $X^{\mathcal{N}_{T}}\in T_{X}$ of the RQI iteration $X_{i+1}=\mathfrak{r}(X_{i},X^{\mathcal{N}_{T}})$ is the solution to the equation

[TABLE]

and the convergence criteria for $\mathrm{F}=\mathrm{L}(X,\mathcal{R}(X))=0$ in theorem 4.1 apply. Let

[TABLE]

Assuming $\mathrm{C}$, $\mathfrak{r}$, $\Pi$, $\mathrm{L},\mathcal{R}$ are smooth. For $\eta\in T_{X}$, set $\nu=\mathfrak{r}(X,.)^{(2)}(0;\eta^{[2]})$ and

[TABLE]

then $\Pi(X)\mathrm{G}_{\mathrm{L}}(X;\eta^{[2]})=0$ if and only if $\Pi(X)\mathrm{G}(X;\eta^{[2]})=0$ in eq. 4.11. Thus, if $\Pi(X_{*})\mathrm{G}_{\mathrm{L}}(X_{*};\eta^{[2]})=0$ for all $\eta\in T_{X_{*}}$ and $\mathrm{F}^{\natural}(X_{*})$ exists for a solution $X_{*}$ of $\mathrm{F}(X)=0$, the RQI above converges cubically. If $\mathrm{F}^{\natural}(X_{*})$ exists then the Rayleigh-Chebyshev iteration $X_{i+1}=\mathfrak{r}(X_{i},X^{\mathcal{N}_{C}})$ with

[TABLE]

converges cubically to $X_{*}$ if $X_{0}\in C_{X_{*}}(\rho)$ for some $\rho>0$. If $\mathrm{L}_{X}(X,\lambda)\in\mathtt{Lin}(\mathcal{E},\mathcal{E})$ is invertible and $\mathrm{C}^{\prime}(X)\mathrm{L}_{X}(X,\lambda)^{-1}\mathrm{L_{\lambda}}(X,\lambda)$ is invertible, $X^{\mathcal{N}_{T}}$ and $X^{\tau}$ could be computed in Schur form

[TABLE]

If $\mathrm{L}$ is affine in $\lambda$, $\mathrm{L}(X,\lambda)=A(X)-H(X)\lambda$ for two functions $A(X)$ and $H(X)$, then $\mathrm{L}_{X}(X,\lambda)=A^{\prime}(X)-H^{\prime}(X)\lambda$. Dropping the variable name $X$ for brevity, if $H^{-}$ is a left inverse of $H$, we can take $\mathrm{L}^{-}_{\lambda}$ to be $-H^{-}$, giving us a Rayleigh quotient $\lambda=H^{-}A$ and $\mathrm{F}=A-HH^{-}A=\Pi A$ with $\Pi=I_{E}-HH^{-}$. Thus, in lemma 1.1, we can use $\mathrm{A}$ instead of $\mathrm{F}$ in the Schur form solution

[TABLE]

5. Cubic convergence and second covariant derivative

In section 1.2, with the projection function $\Pi$, we defined the vector bundle $E_{\Pi}=\cup_{X\in\mathcal{M}}\operatorname{Im}(\Pi(X))$ with a connection $(\nabla_{\mathtt{c}}\mathrm{F})(X)=\Pi(X)\mathrm{F}^{\prime}(X,\mathtt{c}(X))\in T_{X}$ for a section $\mathrm{F}$, $\Pi(X)\mathrm{F}(X)=\mathrm{F}(X)$ and a vector field $\mathtt{c}$, $\mathtt{c}(X)\in T_{X}$. In the unconstrained case, the Chebyshev iteration

[TABLE]

converge cubically, as seen in section 3. It is natural to ask if there is a relationship between the Rayleigh-Chebyshev term and $\nabla$ in the constrained case. The answer is yes, with the concept of a second covariant derivative.

If $\mathrm{C}$ is the constrained function for $\mathcal{M}$, for two vector fields $\mathtt{c}_{1},\mathtt{c}_{2}$, functions from $\mathcal{M}$ to $\mathcal{E}$ satisfying $\mathrm{C}^{\prime}(X)\mathtt{c}_{i}(X)=0$, $i=1,2,X\in\mathcal{M}$, then $\mathtt{c}_{2}^{\prime}(X;\mathtt{c}_{1}(X))$ is not a vector field, as taking derivative of $\mathrm{C}^{\prime}(X)\mathtt{c}_{2}(X)=0$ in direction $\mathtt{c}_{1}(X)$ does not give $\mathrm{C}^{\prime}(X)\mathtt{c}_{2}^{\prime}(X;\mathtt{c}_{1}(X))=0$ but

[TABLE]

An adjustment of the form $\mathtt{c}_{2}^{\prime}(X,\mathtt{c}_{1}(X))+\Gamma(X,\mathtt{c}_{1}(X),\mathtt{c}_{2}(X))$ is needed to get a vector field, called a covariant derivative of $\mathtt{c}_{2}$ at $X$ along $\mathtt{c}_{1}$. As seen in section 1.2, a projection function to $T\mathcal{M}$ defines a covariant derivative (connection). A retraction $\mathfrak{r}$ also defines a covariant derivative $\nabla^{\mathfrak{r}}$:

Proposition 5.1**.**

For a retraction $\mathfrak{r}$ and two vector fields $\mathtt{c}_{1},\mathtt{c}_{2}$, $\nabla_{\mathtt{c}_{1}}^{\mathfrak{r}}\mathtt{c}_{2}$ defined by

[TABLE]

is a vector field, thus $\nabla^{\mathfrak{r}}$ is a covariant derivative.

Proof.

Differentiate $\mathrm{C}^{\prime}(\mathfrak{r}(X,t\mathtt{c}_{1}(X)))\mathfrak{r}(X,.)^{\prime}(t\mathtt{c}_{1}(X);\mathtt{c}_{2}(X))=0$ (from lemma 2.1) in $t$ then set $t=0$ we get

[TABLE]

Subtracting from eq. 5.1, we get

[TABLE]

This shows $\nabla_{\mathtt{c}_{1}}^{\mathfrak{r}}\mathtt{c}_{2}(X)$ is tangent to $\mathcal{M}$ at $X$. ∎

The association of a connection with a retraction is discussed in [6, section 8.1.1]. Here, we give an explicit formula. A naive attempt for a geometric expression of the Chebyshev term is to consider an adjustment like $\frac{1}{2}\nabla_{\eta}(\nabla_{\eta}\mathrm{F})$ , with $\nabla$ is the connection from $\Pi$ in section 1.2 and $\eta$ is the Newton increment. To evaluate the outer $\nabla_{\eta}$, we need the inner to be a section, for this we need to replace $\eta$ with a vector field. For two vector fields $\mathtt{c}_{1},\mathtt{c}_{2}$

[TABLE]

The last term involves $\mathtt{c}_{2}^{\prime}(X;\mathtt{c}_{1}(X))$, hence dependent on the extensions of $\eta$ to vector fields. But $\nabla_{\nabla_{\mathtt{c}_{1}}^{\mathfrak{r}}\mathtt{c}_{2}}\mathrm{F}(X)=\Pi(X)\mathrm{F}^{\prime}(X;\nabla_{\mathtt{c}_{1}}^{\mathfrak{r}}(\mathtt{c}_{2}X))$ is well-defined, hence $\nabla^{2}_{\mathtt{c}_{1},\mathtt{c}_{2}}\mathrm{F}(X)$ given below is not dependent on derivatives of $\mathtt{c}_{1}$ or $\mathtt{c}_{2}$, it is the second covariant derivative using $\nabla$ and $\nabla^{\mathfrak{r}}$

[TABLE]

Therefore, we can define $\frac{1}{2}\nabla^{2}_{\eta,\eta}\mathrm{F}(X)$ as above, for vector fields satisfying $\mathtt{c}_{1}(X)=\mathtt{c}_{2}(X)=\eta$. [16] uses this expression as the third-order adjustment for the geodesic retraction. Compare with eq. 4.11, we have an extra term $\Pi(X)\Pi^{\prime}(X,\eta)\mathrm{F}^{\prime}(X)\eta$. Differentiate $\Pi(Y)\Pi^{\prime}(Y;\mathtt{c}(Y))\mathrm{F}(Y)=0$ (eq. A.1) in direction $\mathtt{c}(Y)$, $Y\in\mathcal{M}$

[TABLE]

At a solution $X_{*}$, $\mathrm{F}(X_{*})=0$, the first three terms are zero, hence $\Pi(X_{*})\Pi^{\prime}(X_{*},\eta)\mathrm{F}^{\prime}(X_{*})\eta=0$ for any tangent vector $\eta$ at $X_{*}$ if $\mathtt{c}(X_{*})=\eta$. Thus, this term does not affect the order of convergence.

6. Applications

The codes for these sections are found in folder colab in [17], implemented in Python and Julia (for two examples showing cubic convergence of Rayleigh-Chebyshev iterations).

6.1. Complex tensor eigenpairs

For our purpose, a tensor is a vector-valued homogeneous function. An example of the real tensor eigenpair problem is finding critical points of a scalar homogeneous function $\hat{\mathcal{T}}$ under the constraint $\hat{\mathcal{B}}=1$. Here, $\hat{\mathcal{B}}$ is a scalar homogeneous polynomial, which leads to the equation $\operatorname{grad}\hat{\mathcal{T}}=\lambda\operatorname{grad}\hat{\mathcal{B}}$. We will focus on equations of the form $\mathcal{T}=\lambda\mathcal{B}$ explained below, the constraints imposed may be unrelated to $\mathcal{B}$.

We consider a complex-valued homogeneous vector function $\mathcal{T}$, which may not be a gradient of a scalar function. Each entry of $\mathcal{T}$ is a scalar homogeneous polynomial of order $m-1$. The coefficients $(t_{i_{1}\cdots i_{m}})_{1\leq i_{j}\leq n,1\leq j\leq m}$ of entry $\mathcal{T}_{i_{1}},i_{1}=1\cdots n$ can be arranged to what is called a multidimensional array (also called a tensor), we call $\mathcal{T}$ a tensor of dimension $n$ order $m$. For $k$ vectors $X_{n-k+1}\cdots X_{n}$ of size $n$ in the base field, contraction means taking sum over the last $k$ indices of $(t_{i_{1}\cdots i_{m}})$ multiply by entries of $X_{j}$’s, $j=n-k+1\cdots n$, the uncontracted entries are written $I$. Write $X^{[k]}$ for $X$ repeated $k$ times. Abuse of notation, $\mathcal{T}(X^{[m-1]})=\mathcal{T}(I,X^{[m-1]})$ is the vector homogeneous function evaluated at $X$, its Jacobian is $(m-1)\mathcal{T}(I,I,X^{[m-2]})$, while $\hat{\mathcal{T}}(X^{[m]})$ is the scalar function in the symmetric case, where $\operatorname{grad}\hat{\mathcal{T}}=m\mathcal{T}$ (We mostly follows [11]).

Let $\mathcal{B}$ be a homogeneous polynomial vector function of order $d-1$, the tensor eigenpair problem in eq. 1.4 solves for pairs $(X,\lambda)$ satisfying $\mathcal{T}(Z^{[m-1]})=\lambda\mathcal{B}(Z^{[d-1]})$.

We will focus on the case $d=2$, much of what is discussed here holds for $d\geq 2$. The number of complex eigenpairs is given as by $\sum_{i=0}^{n-1}(m-1)^{i}(d-1)^{n-i-1}$ ([13, 12]), with certain counting convention described there. An important case is the Z-pairs [3], [18], where $\mathcal{B}(Z)=Z$

[TABLE]

with the normalization $Z^{\mathsf{T}}Z=1$ in the real case [3]. For the complex case, we normalize by $Z^{*}Z=1$ ($*$ is the Hermitian transpose), and call the pairs unitary Z-pairs (UZ-pairs).

There are a few methods to find all real eigenpairs [13, 19], with the first citation using the homotopy method is state-of-the-art, which also finds complex pairs. In [11], the authors proposed a RQI-type algorithm, called O-NCM, which empirically finds all real Z-pairs in less time than the homotopy method, however, there is no certification that all real pairs are found, as there is no formula to count real pairs. Using the RQI method described in this paper, we extend the O-NCM method to the complex case, using the count $\sum(m-1)^{i}$ for complex (Z-pairs) to decide if we have found all eigenpairs, if all pairs are distinct, up to scaling.

If $m>2$, we can assume $\lambda$ is real, or $\mathcal{E}_{L}=\mathbb{R}$, by scaling $\lambda$ by $\exp(-i\theta)$ and $Z$ by $\exp(-i\theta/(m-2))$, with $\theta$ is the angle of $\lambda$ in polar coordinates. We treat $Z$ as a real vector of dimension $2n$. From $\mathrm{C}(Z)=Z^{*}Z-1$, we have $\mathrm{C}^{\prime}(Z)\eta=2\text{Re}(Z^{*}\eta)$, thus a tangent vector satisfies $\text{Re}(Z^{*}\eta)=0$. We have $\mathrm{L_{\lambda}}(Z,\lambda)(\delta)=-Z\delta$, $\delta\in\mathbb{R}$. A left inverse is $\mathrm{L_{\lambda}}^{-}:\omega\mapsto-\text{Re}(Z^{*}\omega)\in\mathcal{E}_{L},\omega\in\mathcal{E}$, thus

[TABLE]

To use theorem 4.2, we solve for the Rayleigh quotient $\lambda=\mathcal{R}(Z)$ from $\mathrm{L_{\lambda}}^{-}(X,\lambda)\mathrm{L}(Z,\lambda)=0$, or $\text{Re}(Z^{*}\mathcal{T}(Z^{[m-1]})-\lambda\text{Re}Z^{*}Z)=0$

[TABLE]

The projection is $\Pi(Z)\omega=\omega-Z\text{Re}(Z^{*}\omega)$ for $\omega\in\mathcal{E}$. The Newton increment $Z^{\mathcal{N}_{T}}$ satisfies

[TABLE]

The Newton form eq. 1.17 is the complex version of the algorithm O-NCM op. cit. The tangent space is the same as the range $E_{Z}$ of $\Pi(Z)$, both are nullspace of $\omega\mapsto\text{Re}Z^{*}\omega$. In eq. 4.21, with $A=\mathcal{T}(Z^{[m-1]})$, $H=Z$, set $\zeta=\mathrm{L}_{X}^{-1}H,\nu=\mathrm{L}_{X}^{-1}A$, the Schur form solution is

[TABLE]

The corresponding Schur form in the real case is different from algorithm NCM in [11]. The authors found O-NCM outperforms NCM. Performance comparison between our RQI Newton and Schur forms is inconclusive, however, we focus on the Schur form for ease of implementation.

In algorithm 2, the initial point $Z_{0}$ on the unitary sphere $Z_{0}^{*}Z_{0}=1$ is chosen randomly. To make sure we only count distinct pairs under the equivalent relation, we keep a table tracking all eigenpairs we have found within our search, identifying pairs where $Z$ is scaled by a $(m-2)$-th root of $1$. This algorithm finds all real pairs if the pairs are isolated. The approach has the advantage of a quadratic convergence algorithm for an individual pair. It does not work well with nonisolated zero, but some remedies could be applied. The speed to find all pairs depends on the distribution of the pairs; generally it works well for randomly generated tensors.

In UZPairsEigenTensor.ipynb in [17] (containing calculations for all examples here), we compare with examples in [19]. For several examples in that paper, RQI outperforms by a large time factor. For nonisolated eigenpairs or infinite pairs, both approaches need modifications. We will not go into details, but note the case where $\mathcal{T}$ is the gradient of the Motzkin’s polynomial $\hat{\mathcal{T}}(Z^{[6]})=z_{0}^{4}z_{1}+z_{0}^{2}z_{1}^{4}+z_{2}^{6}-3z_{0}^{2}z_{1}^{2}z_{2}^{2}$ [12, example 5.9], we found new complex pairs. The authors found 25 eigenpairs out of the 31 expected pairs. The 6 missing pairs are complex pairs, (satisfying $Z^{\mathsf{T}}Z=0$, the authors used the E-pair normalization $Z^{\mathsf{T}}Z=1$) with eigenvalues $\frac{1}{12}$, eigenvectors of the form $(z_{0},-\bar{z}_{0},0)$, $4z_{0}^{4}=-1$ (2 equivalent pairs) and $\frac{3}{16}$, with eigenvectors $(\pm\frac{\sqrt{-1}}{2},\pm\frac{\sqrt{-1}}{2},\frac{\sqrt{2}}{2})$ (4 equivalent pairs). Up to equivalence, we found $6$ real pairs with eigenvalue [math], while the paper counted 14. A perturbation, adding a small diagonal tensor shows there are indeed 14 ($8$ complex and $6$ real) pairs of small eigenvalues that collapse to the $6$ real pairs, the eigenvectors $(1,0,0)$ and $(0,1,0)$ each has multiplicity $5$, and each eigenvector $(\pm\frac{1}{\sqrt{3}},\pm\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}})$ has multiplicity $1$ for the count of $14$. See [13, table 7] or [17] for the remaining pairs.

6.1.1. $\mathcal{B}$ eigenpairs

For the general problem $\mathrm{L}(X,\lambda)=\mathcal{T}(X^{[m-1]})-\mathcal{B}(X^{[d-1]})\lambda$, it is easy to derive several RQIs based on different constraints, e.g. by using a linear constraint $\mathrm{C}^{\mathsf{T}}Z=1$, where $\mathrm{C}$ is a randomly generated real vector [13], thus $\mathcal{E}_{L}=\mathbb{C}$, which we will not go into detail.

We discuss briefly the case where $\mathcal{T},\mathcal{B}$ are symmetric, $\mathcal{T}=(1/m)\hat{\mathcal{T}}^{\prime},\mathcal{B}=(1/d)\hat{\mathcal{B}}^{\prime}$ where $\hat{\mathcal{T}},\hat{\mathcal{B}}$ are scalar homogeneous polynomials of order $m$ and $d$, over $\mathbb{R}$. Assuming $\hat{\mathcal{B}}(X^{[d]})=1$ has a nonempty solution set. We can use the constraint $\hat{\mathcal{B}}(X^{[d]})=1$ with the rescaling/projection retraction discussed in section 2.3. The tangent space at $X$ is defined by the equation $\mathcal{B}(X^{[d-1]})^{\mathsf{T}}\eta=0$. We can take $\mathrm{L_{\lambda}}^{-}(X,\lambda)=X^{\mathsf{T}}$, resulting in the Rayleigh quotient $\mathcal{R}(X)=\hat{\mathcal{T}}(X^{[m]})$. The projection is defined by $\Pi(X)\omega=\omega-\mathcal{B}(X^{[d-1]})X^{\mathsf{T}}\omega$, hence $\operatorname{Im}(\Pi(X))=\mathcal{E}_{X}$ is defined by $X^{\mathsf{T}}\omega=0$. The RQI step can be computed as above.

To analyze cubic convergence, $\mathcal{R}^{\prime}(X_{*})\eta=m\mathcal{T}(X_{*}^{[m-1]})^{\mathsf{T}}\eta=m\lambda\mathcal{B}(X_{*}^{[d-1]})^{\mathsf{T}}\eta=0$, at an eigenpair $(X_{*},\mathcal{R}(X_{*})$ and $\eta\in T_{X_{*}}$. By eq. 2.17, $\mathfrak{r}_{\eta\eta}(X,0)\eta^{[2]}$ is proportional to $X$,

[TABLE]

From eq. 4.16, only $\Pi\mathrm{L}_{XX}$ remains in $\mathrm{G}_{\mathrm{L}}$. Thus, for $m=d=2$, we have a cubically convergence iteration for the generalized eigenvalue problem $AX-BX\lambda=0$ for symmetric matrices $A$ and $B$ ($\mathcal{T}(X)=AX,\mathcal{B}=BX$).

[TABLE]

if $B$ is invertible, where $X_{rscl}$ is the rescaling of $X\in\mathbb{R}^{n}$ to satisfy $X^{\mathsf{T}}BX=1$. Table 1 compares RQI and Rayleigh-Chebyshev iteration for $m=3,d=2$ (see also [17]).

6.2. Constrained optimization

Recall in this problem, $\mathrm{L}(X,\lambda)=\nabla f(X)-\mathrm{C}^{\prime}(X)^{\mathsf{T}}\lambda$ for a real-valued function $f$, $\mathrm{L_{\lambda}}(X,\lambda)=-\mathrm{C}^{\prime}(X)^{\mathsf{T}}$. Let $\mathrm{C}_{X}=\mathrm{C}^{\prime}(X)$ then the Moore-Penrose left inverse is $\mathrm{L_{\lambda}}^{-}(X)=-(\mathrm{C}_{X}\mathrm{C}_{X}^{\mathsf{T}})^{-1}\mathrm{C}_{X}$, and the equation $\mathrm{L_{\lambda}}^{-}(X)\mathrm{L}(X,\lambda)=0$ implies $\lambda=(\mathrm{C}_{X}\mathrm{C}_{X}^{\mathsf{T}})^{-1}\mathrm{C}_{X}\nabla f(X)$ which is known classically. This implies $\Pi(X)=I_{\mathcal{E}}-\mathrm{C}_{X}^{\mathsf{T}}(\mathrm{C}_{X}\mathrm{C}_{X}^{\mathsf{T}})^{-1}\mathrm{C}_{X}$. Newton’s method of this type has been studied before [20, 21], using the Newton form to solve for $X^{\mathcal{N}_{T}}$. Let $Q$ be an orthogonal basis of the range of $\Pi(X)$, (which is the nullspace of $\mathrm{C}_{X}$), constructed from the QR decomposition of $\mathrm{C}_{X}^{\mathsf{T}}$. Then $Q^{\mathsf{T}}Q=I$ and $\Pi(X)=QQ^{\mathsf{T}}$. In eq. 1.17, $Q_{\Pi}=Q_{T}=Q$, and the Newton increment reduces to

[TABLE]

with $Q^{\mathsf{T}}\nabla f(X)$ called the projected gradient, $Q^{\mathsf{T}}\mathrm{L}_{X}(X,\lambda)Q$ the projected Hessian (the manifold literature calls $\Pi(X)\nabla f(X)$ the projected gradient and $\Pi(X)\mathrm{L}_{X}(X,\lambda)$ the projected Hessian). The Schur complement method is also used, see [22, Chapter 10].

6.3. Eigenvector and eigenvector with a constant term

We have shown our RQI with the constraint $X^{\mathsf{T}}X-1=0$ is the classical RQI for the eigenvalue problem, with $\mathrm{L}_{X}(X,\lambda)\eta=(A-\lambda I)\eta$, $\mathrm{L_{\lambda}}(X,\lambda)=-X$. With $\mathcal{R}(X)=X^{\mathsf{T}}AX$, $\mathcal{R}^{\prime}(X)\eta=\eta^{\mathsf{T}}AX+X^{\mathsf{T}}A\eta=\eta^{\mathsf{T}}(A+A^{\mathsf{T}})X$. Thus,

[TABLE]

When $A$ is normal, at an eigenvector $(X_{*},\lambda)$ we have $AX_{*}=A^{\mathsf{T}}X_{*}=\lambda X_{*}$, hence $\mathcal{R}^{\prime}(X_{*})\eta=2\lambda\eta^{\mathsf{T}}X_{*}=0$, thus $\mathrm{G}(X_{*})=0$, therefore we have cubic convergence.

When $A$ is not normal, at an eigenvector $X_{*}$, only $2\eta^{\mathsf{T}}A^{\mathsf{T}}X_{*}$ may be nonzero. To the extent we have tested, the Rayleigh-Chebyshev iteration seems not competitive for generic matrices, as although it eventually has cubic convergence, its area of convergence seems smaller than RQI.

For nonnormal matrices, the two-sided Rayleigh quotient algorithm in [9] has cubic convergence. It is a special case of the two-sided nonlinear eigenvalue algorithm in section 6.4.

For $\mathtt{b}\in\mathbb{R}^{n},\mathtt{b}\neq 0$, an RQI for $\mathrm{L}(X,\lambda)=AX-\lambda X-\mathtt{b}=0$ could be derived similarly. $\mathrm{L}_{X},\mathrm{L_{\lambda}},\mathrm{L_{\lambda}^{-}}$ are as above, the Rayleigh quotient is $\lambda=X^{\mathsf{T}}AX-X^{\mathsf{T}}\mathtt{b}$. In the resulting Schur form iteration, $\mathrm{L}_{X}^{-1}\mathrm{L}=X-(A-\lambda I)^{-1}\mathtt{b}$. Define $Z_{unit}=Z/\|Z\|$ for $Z\in\mathbb{R}^{n}$ then

[TABLE]

The Chebyshev iteration is more stable than the eigenvector case as $(A-\lambda I)$ is often invertible near a solution. While requiring less iterations, it still does not outperform the RQI in total time. This does not rule out applications to special matrices. The expression for $\mathrm{G}_{\mathrm{L}}$ is unchanged, with $\eta$ in eq. 6.6 and $\mathcal{R}^{\prime}(X,\eta)=\eta^{\mathsf{T}}(A+A^{\mathsf{T}})X-\eta^{\mathsf{T}}\mathtt{b}$. The RQI does not converge cubically even for symmetric matrices as $\eta^{\mathsf{T}}AX=\eta^{T}\mathtt{b}$ near a solution. The Chebyshev increment is

[TABLE]

and the Chebyshev iteration is $X_{i+1}=(X_{i}+\eta^{\mathrm{C}})_{unit}$. The last two columns of table 1 shows the RQI and Rayleigh-Chebyshev iterations for this case.

6.4. Nonlinear eigenvalue problem

A nonlinear RQI appeared in algorithm 4.9 in [2], proposed in [10]. We consider $\mathcal{E}=\mathbb{R}^{n}$, $\mathcal{E}_{L}=\mathbb{R}$. Recall $\mathrm{L}(X,\lambda)=\mathrm{P}(\lambda)X$ where $\mathrm{P}$ is a square matrix with polynomial entries.

Assume the unit sphere constraint $\mathrm{C}(X)=\frac{1}{2}(X^{\mathsf{T}}X-1)$. For $\eta\in\mathcal{E}$, we have $\mathrm{C}^{\prime}(X;\eta)=X^{\mathsf{T}}\eta$ and $\mathrm{L}_{X}(X,\lambda)\eta=\mathrm{P}(\lambda)\eta$ for $\eta\in\mathcal{E}$, $\mathrm{L_{\lambda}}(X,\lambda)\delta=\mathrm{P}^{\prime}(\lambda)X\delta$ for $\delta\in\mathbb{R}$. Assuming $X^{\mathsf{T}}\mathrm{P}^{\prime}(\lambda)X\neq 0$, $\mathrm{L_{\lambda}^{-}}(X,\lambda)\omega=\frac{X^{\mathsf{T}}}{X^{\mathsf{T}}\mathrm{P}^{\prime}(\lambda)X}$ is a left inverse of $\mathrm{L_{\lambda}}(X,\lambda)$. The equation $\mathrm{L_{\lambda}^{-}}(X,\lambda)\mathrm{L}(X,\lambda)=0$ for the Rayleigh quotient implies $X^{\mathsf{T}}\mathrm{P}(\lambda)X=0$, and $\lambda=\mathcal{R}(X)$ is solved from here. The projection is

[TABLE]

Applying algorithm 1 in Schur form, $\mathrm{L}_{X}(X,\lambda)^{-1}\mathrm{L}(X,\lambda)=X,\mathrm{L}_{X}(X,\lambda)^{-1}\mathrm{L_{\lambda}}(X,\lambda)=P(\lambda)^{-1}P^{\prime}(\lambda)X$,

[TABLE]

Thus, $\eta+X$ is proportional to $P(\lambda)^{-1}P^{\prime}(\lambda)X$, and together with the retraction, the iteration could be given as $X_{i+1}=\frac{P(\lambda_{i})^{-1}P^{\prime}(\lambda_{i})X_{i}}{|P(\lambda_{i})^{-1}P^{\prime}(\lambda_{i})X_{i}|}$. For cubic convergence analysis, from eq. 4.16

[TABLE]

At a solution, $P(\lambda_{*})X_{*}=0$, thus the last term vanishes. From the implicit function theorem

[TABLE]

If $\mathrm{P}$ is normal, $\mathcal{R}^{\prime}(X_{*},\eta)=0$ at a solution, hence $\mathrm{G}_{\mathrm{L}}(X_{*})=0$, implying cubic convergence.

The two-sided iteration. Consider $\mathcal{E}=\mathbb{R}^{2n}$ and $\mathcal{E}_{L}=\mathbb{R}^{2}$, $X=\begin{bmatrix}u^{\mathsf{T}}&v^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}\in\mathcal{E}$, $u,v\in\mathbb{R}^{n}$ and $\vec{\lambda}=[\lambda_{1},\lambda_{2}]^{\mathsf{T}}\in\mathcal{E}_{L}$. Set $\mathrm{C}(X)=\frac{1}{2}[u^{\mathsf{T}}u-1,v^{\mathsf{T}}v-1]^{\mathsf{T}}\in\mathbb{R}^{2}$ and

[TABLE]

In matrix form, $\mathrm{L}_{X}(X,\vec{\lambda})=\hat{\mathrm{P}}(\vec{\lambda})$, $\mathrm{L_{\lambda}}$ and a left inverse $\mathrm{L_{\lambda}^{-}}$ are

[TABLE]

assuming $v^{\mathsf{T}}\mathrm{P}^{\prime}(\lambda_{1})u,u^{\mathsf{T}}\mathrm{P}^{\prime}(\lambda_{2})v\neq 0$. Then, $\vec{\lambda}=\mathcal{R}(X)$ satisfies $\mathrm{L_{\lambda}^{-}}(X,\vec{\lambda})\mathrm{L}(X,\vec{\lambda})=0$

[TABLE]

We have $\mathrm{C}^{\prime}(X)=\operatorname{bdiag}(u^{\mathsf{T}},v^{\mathsf{T}})\in\mathbb{R}^{2\times 2n}$ where $\operatorname{bdiag}$ denotes a rectangular block diagonal matrix and $\mathrm{L}_{X}^{-1}\mathrm{L}=X$. In algorithm 1

[TABLE]

where $\eta$ is the Newton increment. Thus, components of $X+\eta$ are proportional to $\mathrm{P}(\lambda_{1})^{-1}\mathrm{P}(\lambda_{1})u$ and $\mathrm{P}(\lambda_{2})^{\mathsf{T}})^{-1}\mathrm{P}(\lambda_{2})^{\mathsf{T}}v$. Hence, the Newton step is

[TABLE]

For $\lambda\in\mathbb{R}$, $v^{\mathsf{T}}\mathrm{P}(\lambda)u=u^{\mathsf{T}}\mathrm{P}(\lambda)^{\mathsf{T}}v$, so in eq. 6.10 we can choose $\lambda_{1}=\lambda_{2}$, the resulting iteration converges to left and right eigenvectors of the same eigenvalue. We have cubic convergence in this case (if $\lambda_{1}\neq\lambda_{2}$, the iteration may still converge quadratically to different eigenvalues). In fact, if $\lambda(u,v)$ solves $v^{\mathsf{T}}\mathrm{P}(\lambda)u=0$, differentiate in the tangent direction $\eta=(\eta_{u},\eta_{v})$

[TABLE]

At an eigenvalue $\lambda_{*}$ with right and left eigenvectors $u_{*},v_{*}$, the first two terms vanish, hence $\lambda^{\prime}(X_{*},\eta)=0$ if $v^{\mathsf{T}}_{*}\mathrm{P}^{\prime}(\lambda_{*})u_{*}\neq 0$. For $X_{*}=(u_{*},v_{*}),\vec{\lambda}_{*}=(\lambda_{*},\lambda_{*})$, we have $\mathrm{G}_{\mathrm{L}}(X_{*},\lambda_{*})=0$ since

[TABLE]

7. Concluding remarks

We gave an effective procedure to derive Rayleigh quotient iterations for constrained equations with nonlinear multipliers and provided a clear analysis of second and third-order convergence. The theory developed here explains classical cubic convergence results. It also gives an effective algorithm to find all complex tensor eigenpairs. In future work, we will provide a Kantorovich’s type theorem [8]. The non-isolated zeros case could also be studied.

Appendix A Various proofs

A.1. Proof of lemma 1.1

By direct substitution, using $PH=0$

[TABLE]

A.2. A lemma for projection functions

Lemma A.1**.**

Let $T_{0}$ and $\mathcal{E}$ be two vector spaces, $\Omega\subset T_{0}$ be an open subset, $\Phi:\Omega\to\mathtt{Lin}(\mathcal{E},\mathcal{E})$ be a smooth projection function. If $Z\in\operatorname{Im}(\Phi(\xi))\subset\mathcal{E}$, hence $\Phi(\xi)Z=Z$, then for $\eta\in T_{0}$

[TABLE]

Proof.

Differentiate $\Phi(\xi)^{2}Z=Z$ in direction $\eta$, we get $\Phi(\xi)\Phi^{\prime}(\xi;\eta)Z+\Phi^{\prime}(\xi;\eta)\Phi(\xi)Z=0$. Apply $\Phi(\xi)$ to both sides we get $2\Phi(\xi)\Phi^{\prime}(\xi;\eta)Z=0.$ ∎

A.3. Proof of proposition 3.1

The first two equalities in 1) are immediate. Identify $\mathcal{E}$ and $T_{0}$ with $\mathbb{R}^{n}$ and $\mathbb{R}^{k}$ so we can define transpose, if $\mathcal{F}^{\natural}$ is defined, then $\Phi(\xi)\mathcal{F}^{\prime}(\xi)$ is of rank $k$, and $B(\xi):=(\Phi(\xi)\mathcal{F}^{\prime}(\xi))^{\mathsf{T}}\Phi(\xi)\mathcal{F}^{\prime}(\xi)\in\mathtt{Lin}(T_{0},T_{0})$ is invertible. The equation for $\mathcal{F}^{\natural}$ could be written $B(\xi)\mathcal{F}^{\natural}(\xi)\omega=(\Phi(\xi)\mathcal{F}^{\prime}(\xi))^{\mathsf{T}}\Phi(\xi)\omega$ for $\omega\in\mathcal{E}$. Since $B$ is differentiable, $\xi\mapsto B(\xi)^{-1}$ is differentiable, thus $\mathcal{F}^{\natural}(\xi)$ is differentiable.

For eq. 3.5, differentiated $\Phi(\xi)\mathcal{F}^{\prime}(\xi)\mathcal{F}^{\natural}(\xi)\omega=\Phi(\xi)\omega$ in $\xi$ in direction $\eta$ then apply $\mathcal{F}^{\natural}(\xi)$

[TABLE]

For 2), since $\mathcal{F}^{\natural}(\xi)\Phi(\xi)=\mathcal{F}^{\natural}(\xi)$, $\Phi(\xi)Y-\mathcal{F}(\xi)=0$ implies

[TABLE]

The converse follows by simplifying $\Phi(\xi)\mathcal{F}^{\prime}(\xi)\mathcal{F}^{\natural}(\xi)[Y-\mathcal{F}(\xi)]=0$, using $\Phi(\xi)\mathcal{F}(\xi)=\mathcal{F}(\xi)$.

Consider $G:(\xi,Y)\mapsto\mathcal{F}^{\natural}(\xi)[Y-\mathcal{F}(\xi)]$, mapping $D_{1}\times T_{0}$ to $T_{0}$. If $(\xi,Y)$ satisfies $\mathcal{F}(\xi)=\Phi(\xi)Y$

[TABLE]

for $\eta\in T_{0}$, where using eq. 3.5 to expand $(\mathcal{F}^{\natural})^{\prime}$, we note $\mathcal{F}^{\natural}(\xi)[Y-\Phi(\xi)Y]=0$, hence the first and last terms in eq. 3.5 vanish, then use eq. A.1. To apply the implicit function theorem for $G$, note that $\mathcal{F}^{\natural}(\xi)\Phi^{\prime}(\xi;.;Z)=0$ for $Z=\mathcal{F}(\xi)$, hence, $I_{T_{0}}-\mathcal{F}^{\natural}(\xi)\Phi^{\prime}(\xi;.;Z)=I_{T_{0}}$ is invertible, and near $Z$ it is still invertible, thus we have eq. 3.7 from

[TABLE]

Differentiate the above in $Y$ in direction $\omega_{1}\in\mathcal{E}$, noting $\xi=\mathcal{F}_{\diamond}(Y)$ (and not constant) we have

[TABLE]

For $Y=\mathcal{F}(\xi)=:Z,\omega_{1}=\omega$, using eqs. 3.8 and A.1 to simplify the above

[TABLE]

Let $\nu=\Phi^{\prime}(\xi;\mathcal{F}^{\natural}(\xi)\omega)\mathcal{F}(\xi)-\mathcal{F}^{\prime}(\xi)\mathcal{F}^{\natural}(\xi)\omega$. Differentiate $\Phi(\xi)\mathcal{F}(\xi)=\mathcal{F}(\xi)$ in direction $\mathcal{F}^{\natural}(\xi)\omega$

[TABLE]

Thus, $\nu=-\Phi(\xi)\mathcal{F}(\xi)^{\prime}\mathcal{F}^{\natural}(\xi)\omega$. By eq. A.1, $\mathcal{F}^{\natural}(\xi)\Phi^{\prime}(\xi;\mathcal{F}^{\natural}(\xi)\omega)\nu=0$. Let $\psi=(\omega+\Phi^{\prime}(\xi;\mathcal{F}^{\natural}(\xi)\omega)Z)$. Use eq. 3.5 and $\mathcal{F}^{\natural}(\xi)\psi=\mathcal{F}^{\natural}(\xi)\omega$ to expand $\mathcal{F}_{\diamond}^{(2)}(Z;\omega^{[2]})$ further to

[TABLE]

This proves eq. 3.9. Higher derivatives of $\mathcal{F}_{\diamond}$ could be computed by repeated differentiation of $\mathcal{F}_{\diamond}^{\prime}$. Inductively, its derivatives could be expressed algebraically in terms of derivatives of $\mathcal{F},\Phi$ and powers of $(I_{T_{0}}-\mathcal{F}^{\natural}(\xi)\Phi^{\prime}(\xi;.)Y)^{-1}$, giving us the $C^{k}$ condition with bounded derivatives.

A.4. Proof of theorem 3.1

From the mean value theorem, if $\mathcal{F}^{\prime}(\xi)$ is bounded by $b_{\mathcal{F}^{\prime}}$ in a domain then

[TABLE]

For 1), from the Taylor formula of $\mathcal{F}_{\diamond}$ at $\xi_{i}$ in 2) of proposition 3.1, set $R_{2}$ to be the Taylor remainder

[TABLE]

in a ball in $D_{1}$ where $\mathcal{F}_{\diamond}$ is defined with bounded derivatives. Boundedness and bilinearity of $\mathcal{F}_{\diamond}^{(2)}$ implies there is $c_{\mathcal{F}_{\diamond}}>0$ such that $|\mathcal{F}_{\diamond}^{(2)}(\xi;\delta^{[2]})|\leq c_{\mathcal{F}_{\diamond}}|\delta|^{2}$ for $\delta\in T_{0}$. Since $\xi_{i+1}=\mathfrak{s}(\xi_{i},\delta_{2,i})$,

[TABLE]

provided $\delta_{2,i}$ is defined and $|\delta_{2,i}|<\rho_{2}$. Thus, in $D_{1}$, $|\delta_{2}|=|\mathcal{F}^{\natural}(\xi)\mathcal{F}(\xi)|\leq\|\mathcal{F}^{\natural}(\xi)\||\mathcal{F}(\xi)|\leq b_{\mathcal{F}^{\natural}}b_{\mathcal{F}^{\prime}}|\xi-\xi_{*}|$ if $\|\mathcal{F}^{\natural}(\xi)\|$ and $|\mathcal{F}^{\prime}(\xi)|$ are bounded by $b_{\mathcal{F}^{\natural}}$ and $b_{\mathcal{F}^{\prime}}$, respectively. Together, we have

[TABLE]

For $\xi_{i}$ to be well-defined for all $i$, we also need $|\delta_{2,i}|<\rho_{2}$. Thus, if $|\xi_{0}-\xi_{*}|<\rho_{3}$ with

[TABLE]

then $\delta_{2,i}$ and $\xi_{i}$ are well-defined, $|\xi_{i}-\xi_{*}|$ decreases with $i$, $\{\xi_{i}\}$ converges quadratically to $\xi_{*}$.

For 2), write $\mathfrak{s}_{2}(\xi_{i},\delta_{2,i})=\frac{1}{2}\mathfrak{s}(\xi_{i},.)^{(2)}(0,\delta_{2,i}^{[2]})+\mathfrak{s}_{3}(\xi_{i},\delta_{2,i})$ and add $\frac{1}{2}\mathcal{F}^{\natural}(\xi_{*})\mathcal{G}(\xi_{*},\delta_{2,i}^{[2]})=0$ in eq. A.2

[TABLE]

if we bound $\mathcal{F}_{\diamond}^{(3)}(.,[\delta_{1},\delta_{2},\delta_{3}])\leq c_{\mathcal{F}_{\diamond,3}}|\delta_{1}||\delta_{2}||\delta_{3}|$, hence the integral by $c_{\mathcal{F}_{\diamond,3}}\int_{0}^{1}(1-t)^{2}|\mathcal{F}(\xi_{i})||\mathcal{F}(\xi_{i})|^{2}dt$ using the mean value theorem for $\mathcal{F}_{\diamond}^{(2)}$ and $\mathcal{F}(\xi_{*})=0$, and similarly for the terms involving $\mathfrak{s}_{3}$ and $|\mathfrak{s}(\xi_{*},.)^{[2]}-\mathfrak{s}(\xi_{i},.)^{[2]}|$, eventually get an expression of the form $c|\xi_{*}-\xi_{i}|^{3}$. It remains to estimate the expression $Q$ below using eq. 3.9, adding $-2\mathcal{F}^{\natural}(\xi_{i})\Phi^{\prime}(\xi_{i},\mathcal{F}^{\natural}(\xi_{i})\mathcal{F}(\xi_{i}))\mathcal{F}(\xi_{i})=0$ from eq. A.1

[TABLE]

The first term is zero since $\mathcal{F}(\xi_{*})=0$, and the other terms are bounded by the mean value theorem, for example, for the last term, consider the mean value of $f(t)=\mathcal{F}^{\natural}(\xi_{*})\Phi^{\prime}(\xi_{t},\mathcal{F}^{\natural}(\xi_{t})\mathcal{F}(\xi_{i}))\mathcal{F}(\xi_{i})$ with $\xi(t)=\xi_{i}+t(\xi_{*}-\xi_{i})$, which gives us a term $(\xi_{*}-\xi_{i})$, together with the two terms $\mathcal{F}(\xi_{i})$ and boundedness of derivatives giving an $O(|\xi_{*}-\xi_{i}|^{3})$ estimate.

For 3), with $R_{3}$ denotes the third-order Taylor remainder of $\mathcal{F}_{\diamond}$ centered at $\mathcal{F}(\xi_{i})$, from eqs. 3.8 and 3.9, (the term with $2\mathcal{F}^{\natural}\Phi^{\prime}\mathcal{F}$ is zero by eq. A.1)

[TABLE]

Let $\mathfrak{s}_{3}$ be the third Taylor remainder for $\mathfrak{s}$, expand $\delta_{3,i}$, add and subtract $\frac{1}{2}\mathcal{F}^{\natural}(\xi_{i})\Phi^{(2)}(\xi_{i},\delta_{2,i}^{[2]})\mathcal{F}(\xi_{i})$

[TABLE]

We show each term of the last expression above is $O(|\xi_{i}-\xi_{*}|^{3})$. For $R_{3}$, in a sufficiently small ball, there is a constant $c_{\mathcal{F}_{\diamond},3}$ such that $|\mathcal{F}_{\diamond}^{(3)}(\xi,\delta^{[3]})|\leq c_{\mathcal{F}_{\diamond},3}|\delta|^{3}$ for $\delta\in T_{0}$ and

[TABLE]

For the next term, bound $\|\mathcal{F}^{\natural}(\xi_{i})\|$ by a constant in a ball, $\Phi^{(2)}(\xi,\delta_{2,i}^{[2]})$ by $b_{\Phi,2}|\delta_{2,i}|^{2}$ for a constant $b_{\Phi,2}$, and $\mathcal{F}(\xi_{i})$ by $b_{\mathcal{F}^{\prime}}|\xi_{i}-\xi_{*}|$. Since $\mathfrak{s}(\xi_{i},.)^{(2)}$ is bilinear, expand $\delta_{3,i}=\delta_{2,i}+(\delta_{3,i}-\delta_{2,i})$, noting $|\delta_{2,i}|<c_{3}|\xi_{*}-\xi_{i}|$, $|\delta_{3,i}-\delta_{2,i}|=\frac{1}{2}|\mathcal{F}^{\natural}(\xi)\mathcal{G}(\xi_{i};\delta_{2,i}^{[2]})\leq c_{4*}|\delta_{2,i}|^{2}\leq c_{4}|\xi_{*}-\xi_{i}|^{2}$ for constants $c_{3},c_{4}$, the third term is bounded by

[TABLE]

The $\mathfrak{s}_{3}$ term is third order. The choice of the initial point for the series to be well-defined is similar to the quadratic case.