A Globally Convergent Inexact Newton-Like Cayley Transform Method for Inverse Eigenvalue Problems

Abstract

We propose an inexact Newton method for solving inverse eigenvalue problems (IEP). This method is globalized by employing the classical backtracking techniques. A global convergence analysis of this method is provided and the R-order convergence property is proved under some mild assumptions. Numerical examples demonstrate that the proposed method is very effective in solving the IEP with distinct eigenvalues.

1. Introduction

In the present paper, we consider inverse eigenvalue problems (IEP) which are defined as follows. Let𝐜=(𝑐₁,𝑐₂,…,𝑐_𝑛)^⊤∈ℝ^𝑛and let {𝐴_𝑖}^𝑛_𝑖=0be a sequence of real symmetric𝑛×𝑛matrices. Define 𝐴(𝐜)=𝐴₀+^𝑛_𝑖=1𝑐_𝑖𝐴_𝑖 (1) and denote its eigenvalues by{𝜆_𝑖(𝐜)}^𝑛_𝑖=1with the increasing order𝜆₁(𝐜)≤𝜆₂(𝐜)≤⋯≤𝜆_𝑛(𝐜). Given𝑛real numbers{𝜆^∗_𝑖}^𝑛_𝑖=1which are arranged in increasing order𝜆^∗₁≤𝜆^∗₂≤⋯≤𝜆^∗_𝑛, the IEP is to find a vector𝐜^∗∈ℝ^𝑛such that 𝜆_𝑖𝐜^∗=𝜆^∗_𝑖,𝑖=1,2,…,𝑛. (2) Such vector𝐜^∗is called a solution of the IEP. This type of inverse problem arises in a variety of applications, for instance, the inverse Toeplitz eigenvalue problem [1, 2], inverse Sturm-Liouville’s problem, inverse vibrating string problem, and the pole assignment problem; see [3–5] and the references therein for more details on these applications.

Define𝑓∶ℝ^𝑛→ℝ^𝑛by 𝑓(𝐜)=𝜆₁(𝐜)−𝜆^∗₁,𝜆₂(𝐜)−𝜆^∗₂,…,𝜆_𝑛(𝐜)−𝜆^∗_𝑛^⊤. (3) Then, solving IEP (2) is equivalent to solving the nonlinear equation𝑓(𝐜)=0onℝ^𝑛. It is clear that𝐜^∗is a solution of the IEP if and only if𝐜^∗is a solution of the equation𝑓(𝐜)=0. Based on this equivalence, Newton’s method can be applied to the IEP, and it converges quadratically [6]. As it is known, each iteration of Newton’s method involves solving a complete eigenproblem for the matrix𝐴(𝐜). To overcome this drawback, different Newton-like methods have been proposed and studied [7, 8]. To alleviate the over-solving problem, Bai et al. presented in [9] an inexact Cayley transform method for solving the nonlinear system𝑓(𝐜)=0. To avoid the computation of the approximate Jacobian equations, Shen and Li proposed in [10, 11] Ulm-like methods for solving the IEPs. However, all these numerical methods for solving the IEP converge only locally.

In this paper, we study the numerical methods with global convergence property for solving the IEP. Since the IEP is a nonlinear equation, we review some classical work on solving the general nonlinear equation 𝑓(𝐜)=0. Among the inexact Newton-type methods where a line search procedure is used, an inexact Newton backtracking method was proposed in [12]. It performed backtracking along the inexact Newton step, and computational results on a large set of test problems have shown its robustness and efficiency [13, 14].

The purpose of the present paper is motivated by the inexact Newton backtracking method proposed in [12], to propose an inexact Newton-type method which combines with the Cayley transform method for solving the IEP. In the backtracking procedure, we use the Rayleigh quotient instead of the classical merit function and therefore reduce the computational cost. Under the classical assumption, which is also used in [7, 9, 10], that the given eigenvalues are distinct and the Jacobian matrix 𝐽(𝐜^∗) is invertible, we show that this method is globally convergent. Some numerical examples are reported to illustrate the effectiveness of the proposed method with distinct eigenvalues.

The paper is organized as follows. In Section 2, a global inexact Newton-type algorithm is proposed. The global convergence analysis is given in Section 3. And finally in Section 4, some numerical examples are given to confirm the numerical effectiveness and the good performance of our algorithm.

2. A Globally Inexact Newton-Like Cayley Transform Method

In this section, we present our algorithm. Letℝ^𝑛×𝑛denote the set of all real𝑛×𝑛matrices. Let ‖⋅‖ and ‖⋅‖_𝐹 denote the 2-norm and the Frobenius norm inℝ^𝑛, respectively. The induced 2-norm inℝ^𝑛×𝑛is also denoted by ‖⋅‖; that is, ‖𝐴‖∶=sup_{𝐱∈ℝ𝑛,𝐱≠0}‖𝐴𝐱‖‖𝐱‖,𝐴∈ℝ^𝑛×𝑛. (4) Then, we have‖𝐴‖≤‖𝐴‖_𝐹for any 𝐴∈ℝ^𝑛×𝑛. For 𝐱∈ℝ^𝑛 and a positive number𝑟, we use𝐁(𝐱,𝑟)to stand for the open ball with radius 𝑟 and center 𝐱. Let𝜆₁(𝐜)≤𝜆₂(𝐜)≤⋯≤𝜆_𝑛(𝐜)be the eigenvalues of matrix𝐴(𝐜) and let{𝐪_𝑖(𝐜)}^𝑛_𝑖=1be the normalized eigenvectors corresponding to{𝜆_𝑖(𝐜)}^𝑛_𝑖=1. Define𝐽(𝐜)=([𝐽(𝐜)]_𝑖𝑗)by [𝐽(𝐜)]_𝑖𝑗∶=𝐪_𝑖(𝐜)^⊤𝐴_𝑗𝐪_𝑖(𝐜),1≤𝑖,𝑗≤𝑛. (5) Let{𝜆^∗_𝑖}^𝑛_𝑖=1be given with𝜆^∗₁≤𝜆^∗₂≤⋯≤𝜆^∗_𝑛and write𝝀^∗=(𝜆^∗₁,𝜆^∗₂,…,𝜆^∗_𝑛)^⊤.

The Cayley transform method for computing approximately the eigenproblem of the matrix𝐴(𝐜)was proposed in [6] and was applied in [9, 10]. We now recall this method and then apply it to our algorithm. Suppose that𝐜^∗is a solution of the IEP. Then, there exists an orthogonal matrix𝑄_∗such that𝑄^⊤_∗𝐴𝐜^∗𝑄_∗=diag𝜆^∗₁,𝜆^∗₂,…,𝜆^∗_𝑛≜Λ^∗. (6) Assume that𝐜, 𝝆, and𝑃are the current approximations of 𝐜^∗, 𝝀^∗, and 𝑄_∗, respectively. Define 𝑒^𝑌∶=𝑃^⊤𝑄_∗, where 𝑌 is a skew-symmetric matrix. Then, (6) can be rewritten as 𝑃^⊤𝐴𝐜^∗𝑃=𝑒^𝑌Λ^∗𝑒^−𝑌=1𝐈+𝑌+2𝑌²+⋯Λ^∗1𝐈−𝑌+2𝑌²+⋯=Λ^∗+𝑌Λ^∗−Λ^∗𝑌+𝑂‖𝑌‖². (7) Based on (7), we define the new approximation𝐜^newof𝐜^∗by neglecting the second-order terms in𝑌: 𝑃^⊤𝐴𝐜^new𝑃=Λ^∗+𝑌Λ^∗−Λ^∗𝑌. (8) By equating the diagonal elements in (8), we have 𝜆^∗_𝑖=𝐩_𝑖(𝐜)^⊤𝐴𝐜^new𝐩_𝑖(𝐜),𝑖=1,2,…,𝑛, (9) where{𝐩_𝑖(𝐜)}^𝑛_𝑖=1are the column vectors of𝑃. Thus, once we get𝐜^newby solving the Jacobian equation, we can obtain𝑌by equating the off-diagonal elements in (8); that is, [𝑌]_𝑖𝑗𝐩∶=^𝑖(𝐜)^⊤𝐴𝐜^new𝐩_𝑗(𝐜)𝜆^∗_𝑗−𝜆^∗_𝑖,1≤𝑖≠𝑗≤𝑛. (10)

In order to update the new approximation𝑃^newof𝑄_∗, we construct an orthogonal matrix𝑈using Cayley's transform 1𝑈∶=𝐈+2𝑌1𝐈−2𝑌⁻¹, (11) and set𝑃^new=𝑃𝑈; that is, we can obtain𝑃^newby solving 1𝐈+2𝑌𝑃^new=1𝐈−2𝑌𝑃. (12) Finally, the new approximations of eigenvalues can be obtained by 𝜌^new_𝑖∶=𝐩_𝑖𝐜^new⊤𝐴𝐜^new𝐩_𝑖𝐜^new,𝑖=1,2,…,𝑛, (13) where{𝐩_𝑖(𝐜^new)}^𝑛_𝑖=1are the column vectors of𝑃^new.

Note that (12) can be computed as follows. Compute𝐻∶=(𝐼−(1/2)𝑌)𝑃^⊤and let𝐡_𝑖be the𝑖th column of𝐻at first. Then, solve𝐰_𝑖, 𝑖=1,2,…,𝑛, iteratively from the𝑛linear systems: 1𝐼+2𝑌𝐰_𝑖=𝐡_𝑖,𝑖=1,2,…,𝑛. (14) Finally, set𝑃^new∶=[𝐰₁,…,𝐰_𝑛]^⊤. Since𝑃is an orthogonal matrix and𝑌is skew-symmetric matrix, we see that𝑃^newmust be orthogonal. To maintain the orthogonality of𝑃^new, (14) cannot be solved inexactly. One could expect that it requires only few iterations to solve each system of (14). This is due to the fact that as{𝐜^𝑘}converges to𝐜^∗,‖𝑌_𝑘‖converges to zero; see [6,(3.64)]. Consequently, the coefficient matrix on the left-hand side of (14) approaches the identity matrix in the limit.

For solving the general nonlinear equation𝑓(𝐱)=0, linesearch techniques [15] are often used to enlarge the convergence basin of a locally convergent method. They are based on a globally convergent method for a problem of the formmin_𝐱∈ℝ^𝑛𝑀(𝐱), where𝑀is an appropriately chosen merit function whose global minimum is a zero of𝑓. In these cases, for a given direction𝐬∈ℝ^𝑛, we have the iteration form 𝐱_𝑘+1=𝐱_𝑘+𝛼𝐬, where𝛼∈(0,1]is such that𝑀(𝐱_𝑘+𝛼𝐬)<𝑀(𝐱_𝑘). The existence of such an𝛼is ensured if there exists an𝛼₀>0such that𝑀(𝐱_𝑘+𝛼𝐬)<𝑀(𝐱_𝑘)for all𝛼<𝛼₀.

In typical linesearch strategies, the step length𝛼is chosen by using so-called backtracking approach. Among the backtracking method, inexact Newton backtracking methods (INB) [12] is a globally convergent process where the𝑘th iteration of an inexact Newton method is embedded in a backtracking strategy. The merit function𝑀of INB usually used𝑀∶=‖𝑓‖; see, for example, [12–14, 16]. Thanks to (3), this will involve computing𝜆_𝑖(𝐜^𝑘)of𝐴(𝐜^𝑘)which are costly to compute. Our intention here is to replace them by the Rayleigh quotient (see (13)). In Section 3 we will show that this replacement retain superlinear and global convergence.

The details of our algorithm for solving the IEP are specified as Algorithm 1. In Step 7, the following sufficient decrease in the merit function‖𝝆(𝐜)−𝝀^∗‖based on the Rayleigh quotient is provided: 𝝆𝐜^𝑘+Δ𝐜^𝑘−𝝀^∗≤1−𝜉1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗. (15) The while loop in Step 7 is also called backtracking loop below. Note that if we set 𝜂¹_𝑘𝝆𝐜∶=^𝑘−𝝀^∗𝝀^∗,𝜂²_𝑘𝝆𝐜∶=^𝑘−𝝀^∗𝝆𝐜^𝑘−1−𝝀^∗, (16) when 𝐜^𝑘 is away from an accumulation point 𝐜^∗, then𝜂¹_𝑘>1or𝜂²_𝑘>1may be possible. Accordingly, the forcing term𝜂_𝑘given in (20) guarantees that𝜂_𝑘∈[0,1)for each𝑘.

3. Convergence Analysis

In this section, we analyze the global behavior of Algorithm 1. We will show that if the given eigenvalues are distinct and if there exists an accumulation point𝐜^∗of{𝐜^𝑘}such that the Jacobian matrix𝐽(𝐜^∗)is invertible, then the iterations are guaranteed to remain near 𝐜^∗ and 𝝆(𝐜^∗)−𝝀^∗=0, 𝐜^𝑘→𝐜^∗as 𝑘→∞. Furthermore, for𝑘sufficiently large, we have the equality𝐜^𝑘+1=𝐜^𝑘+Δ𝐜^𝑘. Thus, we obtain that the ultimate rate of convergence is𝛽which depends on the choices of the𝜂_𝑘given in (20).

Algorithm 1 (Inexact Newton-Like Backtracking Cayley Transform Method for IEP). For any 𝐜⁰∈ℝ^𝑛, 𝜂_max∈[0,1), 𝜉∈(0,1), 0<𝜃_min<𝜃_max<1. Compute the orthonormal eigenvectors {𝐪_𝑖(𝐜⁰)}^𝑛_𝑖=1 and the eigenvalues {𝜆_𝑖(𝐜⁰)}^𝑛_𝑖=1 of 𝐴(𝐜⁰). Let 𝑃₀∶=𝐩₁𝐜⁰,𝐩₂𝐜⁰,…,𝐩_𝑛𝐜⁰=𝐪₁𝐜⁰,𝐪₂𝐜⁰,…,𝐪_𝑛𝐜⁰,𝝆𝐜⁰∶=𝜌₁𝐜⁰,𝜌₂𝐜⁰,…,𝜌_𝑛𝐜^0⊤=𝜆₁𝐜⁰,𝜆₂𝐜⁰,…,𝜆_𝑛𝐜^0⊤. (17) For𝑘=0,1,2,…until convergence do the following.
Step 1. Form[𝐽_𝑘]_𝑖𝑗=𝐩_𝑖(𝐜^𝑘)^⊤𝐴_𝑗𝐩_𝑖(𝐜^𝑘)for1≤𝑖,𝑗≤𝑛.
Step 2. SolveΔ𝐜^𝑘inexactly from the approximate Jacobian equation: 𝐽_𝑘Δ𝐜^𝑘+𝝆𝐜^𝑘−𝝀^∗=0 (18) such that 𝐽_𝑘Δ𝐜^𝑘+𝝆𝐜^𝑘−𝝀^∗≤𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗, (19) where𝜂₀∈(0,1)and 𝜂_𝑘∶=min𝝆(𝐜^𝑘)−𝝀^∗𝛽𝝀^∗𝛽,𝝆(𝐜^𝑘)−𝝀^∗𝛽𝝆(𝐜^𝑘−1)−𝝀^∗𝛽,𝜂_max,𝛽∈(1,2],𝑘=1,2,…. (20)
Step 3. SetΔ𝐜^𝑘=Δ𝐜^𝑘and𝜂_𝑘=𝜂_𝑘.
Step 4. Form the skew-symmetric matrix𝑌_𝑘by (10) with𝐜^new=𝐜^𝑘+Δ𝐜^𝑘.
Step 5. Compute matrix 𝑃_𝑘∶=[𝐩₁(𝐜^𝑘+Δ𝐜^𝑘),…,𝐩_𝑛(𝐜^𝑘+Δ𝐜^𝑘)] by solving (12).
Step 6. Compute 𝝆(𝐜^𝑘+Δ𝐜^𝑘)=(𝜌₁(𝐜^𝑘+Δ𝐜^𝑘),…,𝜌_𝑛(𝐜^𝑘+Δ𝐜^𝑘))^⊤ by (13).
Step 7. While ‖𝝆(𝐜^𝑘+Δ𝐜^𝑘)−𝝀^∗‖>(1−𝜉(1−𝜂_𝑘))‖𝝆(𝐜^𝑘)−𝝀^∗‖ do the following: choose𝜃∈𝜃_min,𝜃_max,thenupdateΔ𝐜^𝑘⟵𝜃Δ𝐜^𝑘,𝜂_𝑘⟵1−𝜃1−𝜂_𝑘. (21)
Step 8. Set𝐜^𝑘+1=𝐜^𝑘+Δ𝐜^𝑘. As Steps 4–6, compute, respectively, the new approximations 𝑃_𝑘+1∶=[𝐩₁(c^𝑘+1),…,𝐩_𝑛(𝐜^𝑘+1)]and𝝆(𝐜^𝑘+1)∶=(𝜌₁(𝐜^𝑘+1),…,𝜌_𝑛(𝐜^𝑘+1))^⊤.

It is worth noting that if{𝐜^𝑘}has no accumulation point, {𝐜^𝑘}has one or more accumulation points and the Jacobian matrix is singular at each of them, or the vectorΔ𝐜^𝑘computed by solving the Jacobian equation (18) is such thatΔ𝐜^𝑘=0, then our algorithm fails.

It is clear that if𝝆(𝐜^𝑘)−𝝀^∗→0as𝑘→∞and𝐜^∗is an accumulation point of{𝐜^𝑘}, then 𝝆(𝐜^∗)−𝝀^∗=0. Let{𝑃_𝑘}be generated by Algorithm 1 (see Step 7) and define𝐸_𝑘∶=𝑃_𝑘−𝑄^∗for each𝑘=0,1,…. The following lemma is taken from [17, Lemma 2].

Lemma 2 (see [17]). For any 𝐜, 𝐜∈ℝ^𝑛, one has ‖𝐴(𝐜)−𝐴(𝐜)‖≤𝐿‖𝐜−𝐜‖, (22) where𝐿∶=(∑^𝑛_𝑖=1‖𝐴_𝑖‖²)^1/2.

Based on Lemma 2, the following lemma is a straightforward application of [9, Lemma 4].

Lemma 3. Assume that𝐜^∗is an accumulation point of{𝐜^𝑘}. Let the given eigenvalues{𝜆^∗_𝑖}^𝑛_𝑖=1be distinct. Then 𝝆𝐜^𝑘−𝝀^∗≤√𝑛𝐿𝐜^𝑘−𝐜^∗+𝜇₂𝐸_𝑘, (23) for any𝑘∈ℕ, where𝜇₂∶=2√𝑛⋅‖𝐴(𝐜^∗)‖.

As shown in [18, Theorem 2.3], in the case when the given eigenvalues{𝜆^∗_𝑖}^𝑛_𝑖=1are distinct, the eigenvalues of𝐴(𝐜)are distinct too for any point𝐜in some neighborhood of𝐜^∗. It follows that the function𝑓(⋅)defined in (3) is analytic in the same neighborhood. However, if𝐜is not near the solution𝐜^∗, the analyticity of the function𝑓(⋅)cannot be guaranteed.

For any symmetric matrix𝑋∈ℝ^𝑛×𝑛, setΛ(𝑋)∶=diag(𝜆₁(𝑋),…,𝜆_𝑛(𝑋)), where𝜆_𝑖(𝑋), 𝑖=1,2,…,𝑛 are the eigenvalues of𝑋. As proved by D. Sun and J. Sun [19, Theorem 4.7], Λ(⋅) is a strongly semismooth function. Based on this result, for any𝐜∈ℝ^𝑛and𝜀>0, there exists𝛿>0sufficiently small such that ‖𝝆(𝐜)−𝝆(𝐜)−𝐽(𝐜)(‖𝐜−𝐜)‖≤𝜀𝐜−𝐜‖²,forany𝐜∈𝐁(𝐜,𝛿). (24)

The following lemma says that the backtracking loop in Step 7 of Algorithm 1 terminates after a finite number of steps.

Lemma 4. There exists𝜂_𝑘^min∈[0,1)such that, for any𝜂_𝑘∈[𝜂_𝑘^min,1), there isΔ𝐜^𝑘satisfying 𝝆𝐜^𝑘+Δ𝐜^𝑘−𝝀^∗≤1−𝜉1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗. (25)

Proof. By using the strong semismoothness of all eigenvalues of a real symmetric matrix [19], for any given𝜀>0, there exists𝛿>0sufficiently small such that 𝝆𝐜^𝑘+Δ𝐜−𝝆𝐜^𝑘−𝐽_𝑘Δ𝐜≤𝜀‖Δ𝐜‖ (26) whenever‖Δ𝐜‖≤𝛿. Choose ()𝜀=1−𝜉1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗Δ𝐜^𝑘, (27) and set 𝜂_𝑘^min∶=max𝜂_𝑘,1−1−𝜂_𝑘𝛿Δ𝐜^𝑘. (28) For any𝜂_𝑘∈[𝜂_𝑘^min,1), letΔ𝐜^𝑘∶=((1−𝜂_𝑘)/(1−𝜂_𝑘))Δ𝐜^𝑘. Then, by the definition of𝜂_𝑘^mingiven in (28), we have Δ𝐜^𝑘=1−𝜂_𝑘1−𝜂_𝑘Δ𝐜^𝑘≤1−𝜂_𝑘^min1−𝜂_𝑘Δ𝐜^𝑘≤1−𝜂_𝑘𝛿Δ𝐜^𝑘⋅11−𝜂_𝑘Δ𝐜^𝑘=𝛿. (29)
On the other hand, by (19), one gets that 𝝆𝐜^𝑘−𝝀^∗+𝐽_𝑘Δ𝐜^𝑘≤𝜂_𝑘−𝜂_𝑘1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗+1−𝜂_𝑘1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗+𝐽_𝑘Δ𝐜^𝑘≤𝜂_𝑘−𝜂_𝑘1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗+1−𝜂_𝑘1−𝜂_𝑘⋅𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗=𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗. (30) This together with (26) and (27) yields that 𝝆𝐜^𝑘+Δ𝐜^𝑘−𝝀^∗≤𝝆𝐜^𝑘+Δ𝐜^𝑘−𝝆𝐜^𝑘−𝐽_𝑘Δ𝐜^𝑘+𝝆𝐜^𝑘−𝝀^∗+𝐽_𝑘Δ𝐜^𝑘≤𝜀Δ𝐜^𝑘+𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗=𝜀1−𝜂_𝑘1−𝜂_𝑘Δ𝐜^𝑘+𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗=(1−𝜉)1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗+𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗=1−𝜉1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗, (31) and the proof is completed.

Next, we give sufficient conditions for Algorithm 1 not to break down in the backtracking loop in Step7.

Lemma 5. If𝝆(𝐜^𝑘)−𝝀^∗≠0and there existsΓ>0for which Δ𝐜^𝑘≤Γ1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗, (32) then the backtracking loop terminates.

Proof. For constantΓin (32) and the given𝜉∈(0,1), choosing𝜀=(1−𝜉)/Γ, there exists𝛿>0sufficiently small such that (26) holds whenever‖Δ𝐜‖≤𝛿. We choose𝜂_𝑘∈[𝜂_𝑘,1]satisfying 1−𝜂_𝑘<𝛿Γ𝝆𝐜^𝑘−𝝀^∗. (33) It follows that Δ𝐜^𝑘≤Γ1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗<𝛿, (34) which gives𝐜^𝑘+Δ𝐜^𝑘∈𝐁(𝐜^𝑘,𝛿). Thus, we have 𝝆𝐜+Δ𝐜^𝑘−𝝆𝐜^𝑘−𝐽_𝑘Δ𝐜^𝑘≤1−𝜉ΓΔ𝐜^𝑘. (35) This together with (19) gives that 𝝆𝐜^𝑘+Δ𝐜^𝑘−𝝀^∗≤𝝆𝐜^𝑘−𝝀^∗+𝐽_𝑘Δ𝐜^𝑘+𝝆𝐜^𝑘+Δ𝐜^𝑘−𝝆𝐜^𝑘−𝐽_𝑘Δ𝐜^𝑘≤𝜂_𝑘𝝆𝐜^𝑘−𝜆^∗+(1−𝜉)1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗=1−𝜉1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗. (36) It follows that the backtracking loop terminates. This completes the proof.

The next lemma gives condition under which (32) is satisfied.

Lemma 6. Assume that the Jacobian matrix𝐽_𝑘is invertible and set𝑀_𝑘∶=‖𝐽_𝑘⁻¹‖. Then there existsΓsuch that (32) holds.

Proof. By using condition (19), one has that Δ𝐜^𝑘≤𝐽_𝑘⁻¹𝐽_𝑘Δ𝐜^𝑘≤𝑀_𝑘𝝆𝐜^𝑘−𝝀^∗+𝝆𝐜^𝑘−𝝀^∗+𝐽_𝑘Δ𝐜^𝑘≤𝑀_𝑘1+𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗=𝑀_𝑘1+𝜂_𝑘1−𝜂_𝑘1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗≤𝑀_𝑘1+𝜂_max1−𝜂_max1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗. (37)
We finish the proof by takingΓ∶=𝑀_𝑘(1+𝜂_max)/(1−𝜂_max).

Lemmas 5 and 6 yield the result below.

Corollary 7. Assume that𝝆(𝐜^𝑘)−𝝀^∗≠0and𝐽_𝑘are invertible. Set𝑀_𝑘∶=‖𝐽_𝑘⁻¹‖and Γ∶=𝑀_𝑘(1+𝜂_max)/ (1−𝜂_max). Then, the backtracking loop in Step 7 of Algorithm 1 terminates with 1−𝜂_𝑘≥min1−𝜂_𝑘,𝛿𝜃_minΓ𝝆𝐜^𝑘−𝝀^∗, (38) for any𝛿>0small enough such that (26) holds whenever‖Δ𝐜‖≤𝛿.

Proof. Suppose that𝜂_𝑘is the final value determined by the while-loop. If𝜂_𝑘=𝜂_𝑘, then (38) is trivial. Assume that𝜂_𝑘≠𝜂_𝑘; that is, the body of the while-loop has been executed at least once. Denoting the penultimate value by𝜂⁻_𝑘, then it follows from (33) that1−𝜂⁻_𝑘≥𝛿/Γ‖𝝆(𝐜^𝑘)−𝝀^∗‖. Thus, we have 1−𝜂_𝑘=𝜃1−𝜂⁻_𝑘≥𝜃_min1−𝜂⁻_𝑘≥𝛿𝜃_minΓ𝝆𝐜^𝑘−𝝀^∗. (39) This completes the proof.

Lemma 8. Assume that𝐜^∗is an accumulation point of{𝐜^𝑘}such that there exists a constantΓindependent of𝑘for which Δ𝐜^𝑘≤Γ1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗, (40) whenever𝐜^𝑘is sufficiently near𝐜^∗and𝑘is sufficient large. Then𝐜^𝑘→𝐜^∗as𝑘→∞.

Proof. Suppose that 𝐜^𝑘↛𝐜^∗. Then, there exists 𝛿>0 sufficiently small such that there are infinitely many 𝑘 for which 𝐜^𝑘∉𝐁(𝐜^∗,𝛿) and (40) hold whenever𝐜^𝑘∈𝐁(𝐜^∗,𝛿)for𝑘sufficiently large.
Since𝐜^∗is an accumulation point of{𝐜^𝑘}, there exists subsequence {𝐜^𝑘𝑖}⊂{𝐜^𝑘} such that 𝐜^𝑘𝑖∈𝐁(𝐜^∗,𝛿/𝑖) for 𝑖 sufficiently large. Choose ℓ_𝑖>0 satisfying𝑘_𝑖+ℓ_𝑖<𝑘_𝑖+1and𝐜^{𝑘𝑖+ℓ𝑖}∉𝐁(𝐜^∗,𝛿). It follows that 𝝆𝐜^{𝑘𝑖+𝑗−1}−𝝀^∗≤1Γ1−𝜂_{𝑘𝑖+𝑗−1}𝝆𝐜^{𝑘𝑖+𝑗−1}−𝝀^∗−𝝆𝐜^{𝑘𝑖+𝑗}−𝝀^∗,𝑗=1,…,ℓ_𝑖. (41) Then, by (40), we have, for𝑖sufficiently large, 𝛿2≤𝐜^{𝑘𝑖+ℓ𝑖}−𝐜^𝑘𝑖≤^{𝑘𝑖+ℓ𝑖−1}_{𝑘=𝑘𝑖}Δ𝐜^𝑘≤^{𝑘𝑖+ℓ𝑖−1}_{𝑘=𝑘𝑖}Γ1−𝜂_𝑘𝝆𝐜^𝑘−𝝀^∗≤Γ𝑡𝝆𝐜^𝑘𝑖−𝝀^∗−𝝆𝐜^𝑘𝑖+1−𝝀^∗. (42) Note that𝐜^𝑘𝑖→𝐜^∗as𝑖→∞. It follows that‖𝝆(𝐜^𝑘𝑖)−𝝀^∗‖−‖𝝆(𝐜^𝑘𝑖+1)−𝝀^∗‖→0, which is a contradiction. Therefore,𝐜^𝑘→𝐜^∗as𝑘→∞.

Lemma 9. Assume that 𝐜^∗ is an accumulation point of {𝐜^𝑘} such that 𝐽(𝐜^∗) is invertible. Set 𝑀=‖𝐽(𝐜^∗)⁻¹‖ and let 𝛿₁ be such that 0<𝛿₁≤1/2𝜇₁𝑀, where 𝜇₁∶=2𝑛⋅max_{1≤𝑗≤𝑛}‖𝐴_𝑗‖. Suppose that ‖𝐸_𝑘‖≤𝛿₁ for 𝑘 sufficiently large. Then, for all 𝑘 sufficiently large, 𝐽_𝑘⁻¹≤𝐽𝐜^∗−11−𝜇₁𝛿₁𝐽(𝐜^∗)⁻¹≤2𝑀. (43)

Proof. By the definitions of[𝐽_𝑘]_𝑖𝑗and[𝐽(𝐜^∗)]_𝑖𝑗, for all𝑘sufficiently large, 𝐽_𝑘𝑖𝑗−𝐽𝐜^∗_𝑖𝑗≤𝐩_𝑖𝐜^𝑘⊤𝐴_𝑗𝐩_𝑖𝐜^𝑘−𝐩_𝑖𝐜^𝑘⊤𝐴_𝑗𝐪_𝑖𝐜^∗+𝐩_𝑖𝐜^𝑘⊤𝐴_𝑗𝐪_𝑖𝐜^∗−𝐪_𝑖𝐜^∗⊤𝐴_𝑗𝐪_𝑖𝐜^∗≤𝐩_𝑖𝐜^𝑘⊤𝐴_𝑗𝐩_𝑖𝐜^𝑘−𝐪_𝑖𝐜^∗+𝐩_𝑖𝐜^𝑘⊤−𝐪_𝑖𝐜^∗⊤𝐴_𝑗𝐪_𝑖𝐜^∗=2𝐴_𝑗𝐩_𝑖𝐜^𝑘−𝐪_𝑖𝐜^∗,1≤𝑖,𝑗≤𝑛. (44) Then, we have, for all𝑘sufficiently large, 𝐽_𝑘−𝐽𝐜^∗≤𝐽_𝑘−𝐽(𝐜^∗)_𝐹≤2𝑛⋅max_{1≤𝑗≤𝑛}𝐴_𝑗⋅max_{1≤𝑖≤𝑛}𝐩_𝑖𝐜^𝑘−𝐪_𝑖𝐜^∗. (45) Noting that 𝐩_𝑖(𝐜^𝑘)−𝐪_𝑖(𝐜^∗)is the𝑖th column of𝐸_𝑘, then‖𝐩_𝑖(𝐜^𝑘)−𝐪_𝑖(𝐜^∗)‖≤‖𝐸_𝑘‖for𝑖=1,…,𝑛. So, for all𝑘sufficiently large, 𝐽𝐜^∗−1𝐽_𝑘−𝐽𝐜^∗≤𝜇₁𝛿₁1𝑀≤2<1. (46) It follows from Banach lemma that𝐽_𝑘is invertible and‖𝐽_𝑘⁻¹‖≤2𝑀for all𝑘sufficiently large. This completes the proof.

Lemmas 6, 8, and 9 yield the result below.

Corollary 10. Assume that 𝐜^∗ is an accumulation point of {𝐜^𝑘} such that 𝐽(𝐜^∗) is invertible. Set 𝑀∶=‖𝐽(𝐜^∗)⁻¹‖ and let 𝛿₁ be determined by Lemma 9. Assume that ‖𝐸_𝑘‖≤𝛿₁ for 𝑘 sufficiently large. Then 𝐜^𝑘→𝐜^∗ as 𝑘→∞.

Proof. By Lemma 9, we know𝐽_𝑘is invertible and‖𝐽_𝑘⁻¹‖≤2𝑀for all𝑘sufficiently large. From the proof of Lemma 6, (40) holds for the constantΓ=2𝑀(1+𝜂_max)/(1−𝜂_max)independent of𝑘. Therefore, 𝐜^𝑘→𝐜^∗ as 𝑘→∞ follows from Lemma 8.

Corollary 11. Assume that𝐜^∗is an accumulation point of{𝐜^𝑘}such that𝐽(𝐜^∗)is invertible. Then, 𝝆(𝐜^𝑘)−𝝀^∗→0as𝑘→∞. Moreover, for all𝑘sufficiently large, one has𝜂_𝑘=𝜂_𝑘.

Proof. IfΔ𝐜^𝑘is computed in the backtracking loop, the backtracking terminates with𝜂_𝑘such that (38) holds. Since𝐜^𝑘→𝐜^∗by Corollary 10, we have𝐜^𝑘∈𝐁(𝐜^∗,𝛿)for all𝑘sufficiently large. Thus, the series ∑^∞_𝑘=0(1−𝜂_𝑘) is divergent. Then, we have _{1≤𝑖<𝑘}1−𝜉1−𝜂_𝑖≤_{1≤𝑖<𝑘}𝑒^{−𝜉(1−𝜂𝑖)}=𝑒^{−𝜉⋅∑𝑘−1𝑖=1(1−𝜂𝑖)}⟶0,𝑘⟶∞, (47) and so 𝝆𝐜^𝑘−𝝀^∗≤1−𝜉1−𝜂_𝑘−1𝝆𝐜^𝑘−1−𝝀^∗≤𝝆𝐜⁰−𝝀^∗∏_{0≤𝑖<𝑘}1−𝜉1−𝜂_𝑖⟶0,𝑘⟶∞. (48) This together with (38) yields that𝜂_𝑘=𝜂_𝑘for all𝑘sufficiently large.

Lemma 12. Assume that𝐜^∗is an accumulation point of{𝐜^𝑘}such that𝐽(𝐜^∗)is invertible. Then, there exists𝛿₂>0such that 𝝆𝐜^𝑘−𝝀^∗−𝐽_𝑘𝐜^𝑘−𝐜^∗≤𝜇₃𝐸_𝑘²_𝐹 (49) whenever𝐜^𝑘∈𝐁(𝐜^∗,𝛿₂)and𝑘sufficiently large, where𝜇₃∶=2max_{1≤𝑖≤𝑛}|𝜆^∗_𝑖|.

Proof. Set𝐻_𝑘∶=𝑄(𝐜^∗)^⊤𝑃_𝑘−𝐈. Then, we have, for all𝑘sufficiently large, 𝑃^⊤_𝑘𝐴𝐜^∗𝑃_𝑘=𝐈+𝐻_𝑘^⊤Λ^∗𝐈+𝐻_𝑘=Λ^∗+Λ^∗𝐻_𝑘+𝐻^⊤_𝑘Λ^∗+𝐻^⊤_𝑘Λ^∗𝐻_𝑘. (50) It follows that, for all𝑘sufficiently large, 𝐩_𝑖𝐜^𝑘⊤𝐴𝐜^∗𝐩_𝑖𝐜^𝑘=𝜆^∗_𝑖+2𝜆^∗_𝑖𝐻_𝑘𝑖𝑖+^𝑛_𝑗=1𝜆^∗_𝑗𝐻_𝑘²_𝑖𝑗,𝑖=1,2,…,𝑛. (51) Set 𝑣_𝑖∶=2𝜆^∗_𝑖𝐻_𝑘𝑖𝑖+^𝑛_𝑗=1𝜆^∗_𝑗𝐻_𝑘²_𝑖𝑗,𝑖=1,2,…,𝑛, (52) and write𝐯=(𝑣₁,𝑣₂,…,𝑣_𝑛)^⊤. It follows that𝐽_𝑘𝐜^∗=𝝀^∗−𝐚^𝑘+𝐯, where 𝐚^𝑘_𝑖=𝐩_𝑖𝐜^𝑘⊤𝐴₀𝐩_𝑖𝐜^𝑘,𝑖=1,2,…,𝑛. (53) In view of that𝐽_𝑘𝐜^𝑘=𝝆(𝐜^𝑘)−𝐚^𝑘for any𝑘∈ℕ, one has that, for all𝑘sufficiently large, 𝐯=𝝆𝐜^𝑘−𝝀^∗−𝐽_𝑘𝐜^𝑘−𝐜^∗. (54) Since 𝐈+𝐻_𝑘𝐻^⊤_𝑘+𝐻_𝑘𝐻^⊤_𝑘=𝐈+𝐻_𝑘𝐈+𝐻_𝑘^⊤=𝑄𝐜^∗⊤𝑃_𝑘𝑃^⊤_𝑘𝑄𝐜^∗=𝐈, (55) we have 𝐻_𝑘𝑖𝑖1=−2^𝑛_𝑗=1𝐻_𝑘²_𝑖𝑗,𝑖=1,2,…,𝑛. (56)
It follows from (52) that, for all𝑘sufficiently large, 𝑣²_𝑖=2𝜆^∗_𝑖𝐻_𝑘𝑖𝑖+^𝑛_𝑗=1𝜆^∗_𝑗𝐻_𝑘²_𝑖𝑗²≤24𝜆^∗_𝑖²𝐻_𝑘²_𝑖𝑖+^𝑛_𝑗=1𝜆^∗_𝑗𝐻_𝑘²_𝑖𝑗². (57) Combining (56) with (57), one has that, for all𝑘sufficiently large, ^𝑛_𝑖=1𝑣²_𝑖≤2^𝑛_𝑖=14𝜆^∗_𝑖²𝐻_𝑘²_𝑖𝑖+^𝑛_𝑖=1^𝑛_𝑗=1𝜆^∗_𝑗𝐻_𝑘²_𝑖𝑗²=2^𝑛_𝑖=1𝜆^∗_𝑖^2𝑛_𝑗=1𝐻_𝑘²_𝑖𝑗²+^𝑛_𝑖=1^𝑛_𝑗=1𝜆^∗_𝑗𝐻_𝑘²_𝑖𝑗²≤2max_{1≤𝑖≤𝑛}𝜆^∗_𝑖²⋅^𝑛_𝑖=1^𝑛_𝑗=1𝐻_𝑘²_𝑖𝑗²+max_{1≤𝑖≤𝑛}𝜆^∗_𝑖²⋅^𝑛_𝑖=1^𝑛_𝑗=1𝐻_𝑘²_𝑖𝑗²≤4max_{1≤𝑖≤𝑛}𝜆^∗_𝑖²⋅^𝑛𝑛_𝑖=1𝑗=1𝐻_𝑘²_𝑖𝑗²=4max_{1≤𝑖≤𝑛}𝜆^∗_𝑖²⋅𝐻_𝑘⁴_𝐹, (58) which gives ‖𝐯‖≤2max_{1≤𝑖≤𝑛}𝜆^∗_𝑖⋅𝐻_𝑘²_𝐹. (59) Note that, 𝐻_𝑘𝐹=𝑄𝐜^∗⊤𝑃_𝑘−𝐈_𝐹=𝑃_𝑘−𝑄(𝐜^∗)_𝐹=𝐸_𝑘𝐹. (60) Therefore, we obtain from (54) that (49) holds for all𝑘sufficiently large. This completes the proof.

Lemma 13. Assume that𝐜^∗is an accumulation point of{𝐜^𝑘}such that𝐽(𝐜^∗)is invertible. Then, there exist𝛿₃>0,𝛿₄∈(0,𝛿₁]sufficiently small such that‖𝐸_𝑘‖≤𝛿₄and 𝐜^𝑘+1−𝐜^∗≤2𝑀𝜇₃𝐸_𝑘²_𝐹+2𝜇₃^𝛽𝝆(𝐜^𝑘)−𝝀^∗𝛽 (61) whenever𝐜^𝑘∈𝐁(𝐜^∗,𝛿₃)and𝑘sufficiently large, where𝛿₁is determined by Lemma 9.

Proof. Thanks to Corollary 11, we have𝜂_𝑘=𝜂_𝑘andΔ𝐜^𝑘=Δ𝐜^𝑘for all𝑘sufficiently large. Thus, the Jacobian equation (18) is equivalent to, for all𝑘sufficiently large, 𝐽_𝑘𝐜^𝑘+1−𝝀^∗+𝐚^𝑘=0. (62) Assume that the residual of this approximate Jacobian equation is defined by𝐫^𝑘; that is, for all𝑘sufficiently large, 𝐫^𝑘=𝐽_𝑘𝐜^𝑘+1−𝝀^∗+𝐚^𝑘. (63) This together with𝐽_𝑘𝐜^∗−𝝀^∗+𝐚^𝑘=𝐯gives, for all𝑘sufficiently large, 𝐽_𝑘𝐜^∗−𝐜^𝑘+1=𝐯−𝐫^𝑘. (64) By (19), Lemmas 9 and 12, we obtain, for all𝑘sufficiently large, 𝐜^𝑘+1−𝐜^∗≤𝐽_𝑘⁻¹‖𝐯‖+𝐫^𝑘≤2𝑀𝜇₃𝐸_𝑘²_𝐹+𝝆𝐜^𝑘−𝝀^∗𝛽+1𝝀^∗𝛽. (65) It follows from Lemma 3 that 𝝆𝐜^𝑘−𝝀^∗≤√𝑛𝐿𝐜^𝑘−𝐜^∗+𝜇₂𝐸_𝑘,∀𝑘∈ℕ. (66) Thus, we can choose𝛿₃>0and0<𝛿₄≤𝛿₁sufficiently small such that √𝑛𝐿𝛿₃+𝜇₂𝛿₄≤1, (67) whenever‖𝐜^𝑘−𝐜^∗‖≤𝛿₃and‖𝐸_𝑘‖_𝐹≤𝛿₄. Therefore, combining this with the definition of𝜇₃given in Lemma 12, (61) follows.

Lemma 14 (see [6]). There exist two positive numbers 𝛿₅ and 𝜔₁ such that, for any orthogonal matrix 𝑃 with ‖𝑃−𝑄(𝐜^∗)‖<𝛿₅, the skew-symmetric 𝑋 defined by 𝑒^𝑋∶=𝑃^⊤𝑄(𝐜^∗) satisfies ‖𝑋‖≤𝜔₁‖𝑃−𝑄(𝐜^∗)‖.

Based on Lemma 14, by using the similar arguments in the proof of [10, Lemma 5], we can obtain the following lemma. If we write 𝑒^𝑋𝑘∶=𝑃^⊤_𝑘𝑄^∗, then there exists 𝐶>0 such that ‖𝑌_𝑘‖≤𝐶(‖𝐜^𝑘+1−𝐜^∗‖+‖𝑋_𝑘‖).

Lemma 15. Suppose that the given eigenvalues {𝜆^∗_𝑖}^𝑛_𝑖=1 are distinct and the Jacobian matrix 𝐽(𝐜^∗) is invertible. Then, there exist 𝜔₂>0 and 0<𝛿₆<min{𝛿₅,1/(1+𝜔₁)𝐶} such that, for𝑘sufficiently large, if ‖𝐜^𝑘+1−𝐜^∗‖≤𝛿₆ and ‖𝐸_𝑘‖≤𝛿₆, then 𝐸_𝑘+1≤𝜔₂𝐜^𝑘+1−𝐜^∗+𝐸_𝑘², (68) where𝜔₁is determined by Lemma 14.

In order to prove our global convergence result for Algorithm 1, we introduce some notations. Let𝜁₁∶=2(√𝑛𝐿+𝜇₂)/𝜇₃and𝜁₂∶=2𝑀(𝜇₃+𝜁^𝛽₁). Set 𝜏∶=min1,2𝑀𝜇₃+𝜁^𝛽₁^{−(1/(𝛽−1))}×𝜔₂𝜁₂+1^{−(1/(𝛽−1))},𝛿∶=min𝜏,𝛿₂,𝛿₃,𝛿₄,𝛿₆. (69)(70) Our main global convergence result is as follows.

Theorem 16. Assume that {𝐜^𝑘} is generated by Algorithm 1. Suppose that the given eigenvalues {𝜆^∗_𝑖}^𝑛_𝑖=1 are distinct and 𝐜^∗ is an accumulation point of {𝐜^𝑘} such that 𝐽(𝐜^∗) is invertible. Then, 𝝆(𝐜^𝑘)−𝝀^∗→0 and 𝐜^𝑘→𝐜^∗ as 𝑘→∞. Moreover, the convergence is of R-order 𝛽.

Proof. It follows immediately from Corollaries 10 and 11 that 𝝆(𝐜^𝑘)−𝝀^∗→0 and 𝐜^𝑘→𝐜^∗ as 𝑘→∞. For the 𝜏 given in (69), there exists 𝑘₀ sufficiently large such that ‖𝐜^𝑘0−𝐜^∗‖≤𝜏 and ‖𝐸_𝑘0‖≤𝜏. Set 𝛾∶=𝛿/𝜏, where 𝛿 and 𝜏 are given in (70) and (69), respectively. Then, 𝛾≤1. We will show that, for all 𝑘≥𝑘₀ sufficiently large, 𝐜^𝑘−𝐜^∗≤𝜏⋅𝛾^{𝛽0𝑘−𝑘},𝐸_𝑘≤𝜏⋅𝛾^{𝛽0𝑘−𝑘}. (71)(72) Suppose that (71) and (72) hold for some 𝑘=ℓ≥𝑘₀. Consider the case 𝑘=ℓ+1. Thanks to Lemma 3, we have, for all ℓ≥𝑘₀, 𝝆𝐜^ℓ−𝝀^∗≤√𝑛𝐿𝐜^ℓ−𝐜^∗+𝜇₂𝐸_ℓ≤√𝑛𝐿𝜏⋅𝛾^{𝛽0ℓ−𝑘}+𝜇₂𝜏⋅𝛾^{𝛽0ℓ−𝑘}=√𝑛𝐿+𝜇₂𝜏⋅𝛾^{𝛽0ℓ−𝑘}. (73) Then, by using Lemma 13, one has that, for all ℓ≥𝑘₀, 𝐜^ℓ+1−𝐜^∗≤2𝑀𝜇₃𝐸_ℓ²_𝐹+2𝜇₃^𝛽𝝆𝐜^ℓ−𝝀^∗≤2𝑀𝜇₃𝜏²𝛾^{𝛽0ℓ−𝑘2}2+2𝑀𝜇₃^𝛽√𝑛𝐿+𝜇₃^𝛽𝜏^𝛽𝛾^{𝛽0ℓ−𝑘+1}≤2𝑀𝜇₃+2√𝑛𝐿+𝜇₂𝜇₃^𝛽𝜏^𝛽𝛾^{𝛽0ℓ−𝑘+1}≤𝜏𝛾^{𝛽0ℓ−𝑘+1}, (74) where the last inequality follows from the definition of𝜏in (69). By Lemma 15, we have, for all ℓ≥𝑘₀, 𝐸_ℓ+1≤𝜔₂𝐜^ℓ+1−𝐜^∗+𝐸_ℓ²≤2𝑀𝜔₂𝜇₃+𝜁^𝛽₁𝜏^𝛽𝛾^{𝛽0ℓ−𝑘+1}+𝜔₂𝜏²𝛾^{𝛽0ℓ−𝑘2}≤𝜔2𝑀𝜇₃+𝜁^𝛽₁+1𝜏^𝛽⋅𝛾^{𝛽0ℓ−𝑘+1}=𝜔₂𝜁₂+1𝜏^𝛽⋅𝛾^{𝛽0ℓ−𝑘+1}≤𝜏⋅𝛾^{𝛽0ℓ−𝑘+1}. (75) Therefore, we conclude that (71) and (72) hold for all𝑘≥𝑘₀. Moreover, we see from (71) that𝐜^𝑘converges to𝐜^∗with R-order𝛽. This completes the proof.

If{𝐜^𝑘}generated by Algorithm 1 converges to a solution at which the Jacobian matrix is invertible, then the ultimate rate of convergence is governed by the choices of the𝜂_𝑘(𝑘=0,1,…)as in the local theory of [9].

4. Numerical Examples

In this section, we illustrate the effectiveness of Algorithm 1 in solving IEP on three examples. The tests were carried out in MATLAB 7.0 running on a PC Intel Pentium P6200 of 2.13 GHz CPU.

The given parameters used in our algorithm were 𝜂₀=0.5, 𝜂_max=0.9, 𝜉=10⁻⁴, 𝜃_min=0.1, and 𝜃_max=0.9. In the while loop, we choose 𝜃∈[𝜃_min,𝜃_max] to minimize ‖𝝆(𝐜^𝑘+𝜃Δ𝐜^𝑘)−𝝀^∗‖ if 80 iterations of the backtracking loop fail to produce the sufficient decrease in ‖𝝆(𝐜)−𝝀^∗‖.

Linear systems (14) and (19) are all solved iteratively by the QMR method [20] using the MATLAB qmr function. In order to guarantee the orthogonality of 𝑃^new in (14), this system is solved up to machine precision eps (≈2.2×10⁻¹⁶). The inner loop stopping tolerance for (18) is given by (20). The stopping criterion of the outer iteration in our algorithm is 𝑃^⊤_𝑘𝐴(𝐜^𝑘)𝑃_𝑘−Λ_∗𝐹≤10⁻¹⁰. (76)

Example 1. This is an inverse Toeplitz eigenvalue problem (see [2] for more details on this inverse problem) with distinct eigenvalues. The basis matrices {𝐴_𝑖}⁵_𝑖=1 are given as follows: 𝐴₁=1000001000001000001000001,𝐴₂=,𝐴0100010100010100010100010₃=0010000010100010100000100,𝐴₄=,𝐴0001000001000001000001000₅=.0000100000000000000010000 (77) The given real eigenvalues and a solution, respectively, are 𝝀^∗=(−5.2361,−1.5876,−0.7639,−0.5555,18.1431)^⊤,𝐜^∗=(2,3,4,5,6)^⊤. (78) In Table 1, we report our numerical results for various starting points: (a)𝐜⁰=(1,2,3,4,5)^⊤,(b)𝐜⁰=(1,5,10,15,20)^⊤,(c)𝐜⁰=(11,12,13,14,15)^⊤,(d)𝐜⁰=(21,38,46,63,81)^⊤,(e)𝐜⁰=(101,112,123,134,145)^⊤, (79) where 𝐜⁰, errs, ite. and 𝐜^∗ stand for the starting point, the error value of the left-hand side of (76) for the last three iterates of the algorithm, the number of outer iteration, and the accumulation point corresponding to the starting point.

Table 1 

Numerical results for Example 1.

Example 2. This is a Toeplitz-plus-Hankel inverse eigenvalue problem (see [1] for more details on this inverse problem) with distinct eigenvalues. The basis matrices{𝐴_𝑖}⁷_𝑖=1are given as follows: 𝐴₁=,𝐴−1000000010000000100000001000000010000000100000001₂=,𝐴0−100000−101000001010000010100000101000001010000010₃=,𝐴00−100000−201000−10001000100010001000100010000000100₄=,𝐴000−100000−201000−200010−1000001010000000100000001000₅=,𝐴0000−100000−201000−200010−200000−100000001000000010000₆=,𝐴00000−100000−201000−200000−200000−200000−10000000100000₇=.000000−100000−200000−200000−200000−200000−200000−1000000 (80) The given real eigenvalues are 𝝀^∗=()−35.4513,−13.6805,−9.5675,−8.5489,8.7666,11.8220,26596^⊤, (81) and𝐜^∗=(2,3,4,5,6,7,8)^⊤is a solution. In Table 2, we report our numerical results for various starting points: (a)𝐜⁰=(1,2,3,4,5,6,7)^⊤,(b)𝐜⁰=(1,3,5,7,9,11,13)^⊤,(c)𝐜⁰=(11,13,15,17,19,21,23)^⊤,(d)𝐜⁰=(50,52,56,58,62,65,68)^⊤,(e)𝐜⁰=(101,102,103,104,106,108,110)^⊤. (82)

Table 2 

Numerical results for Example 2.

Example 3 (see [10]). Given 𝐵=[𝑏_𝑖𝑗]_8×8∶=𝐈+𝑉𝑉^⊤, where 𝑉=1−1−3−5−611−2−5−171−1−1518111201−12011130−12.50.20.30.50.62−0.20.30.50.8_8×5, (83) the basis matrices {𝐴_𝑘}⁸_𝑘=1 are defined from 𝐵 as follows: 𝐴_𝑘=𝑏_𝑘𝑘𝑒_𝑘𝑒^⊤_𝑘+^𝑘−1_𝑗=1𝑏_𝑘𝑗𝑒_𝑘𝑒^⊤_𝑘+𝑒_𝑗𝑒^⊤_𝑘,𝑘=1,2,…,8, (84) where𝑒_𝑘is the𝑘th column of the identity matrix𝐈. The given real eigenvalues and a solution, respectively, are 𝝀^∗=()×()1.0𝑒+04−1.1672,−0.3139,−0.0140,−0.0008,0.0003,0.0075,0.0327,1.7372^⊤, (85) and𝐜^∗=(2,3,4,5,6,7,8,9)^⊤. Table 3 shows the numerical results for the following various starting points: (a)𝐜⁰=(1,1,1,1,1,1,1,1)^⊤,(b)𝐜⁰=(2.5,2.5,2.5,2.5,2.5,2.5,2.5,2.5)^⊤,(c)𝐜⁰=(1,2,3,4,5,6,7,8)^⊤,(d)𝐜⁰=(2,4,6,8,10,12,14,16)^⊤, (86) where 𝐜^∗₁=(2,3,4,5,6,7,8,9)^⊤,𝐜^∗₂=()−2.7065,−7.0690,8.4130,1.1038,−3.6637,8.8346,5.5078,9.1273^⊤,𝐜^∗₃=()1.2877,2.4677,4.5379,−6.1037,1.7930,1.1093,1.1024,5.3320^⊤ (87) are the three accumulation points of this problem.

Table 3 

Numerical results for Example 3.

We observe from Tables 1–3 that our algorithm is convergent for different starting points. We also see that our algorithm converges to a solution of the IEP, which is not necessarily equal to the original one. An interesting question is to consider the performance of the algorithm when the starting point𝐜⁰is chosen as random vector, which needs future study.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The second author’s work was supported in part by the National Natural Science Foundation of China (Grant no. 61170109).