Newton-Krylov Type Algorithm for Solving Nonlinear Least Squares Problems

Abstract

The minimization of a quadratic function within an ellipsoidal trust region is an important subproblem for many nonlinear programming algorithms. When the number of variables is large, one of the most widely used strategies is to project the original problem into a small dimensional subspace. In this paper, we introduce an algorithm for solving nonlinear least squares problems. This algorithm is based on constructing a basis for the Krylov subspace in conjunction with a model trust region technique to choose the step. The computational step on the small dimensional subspace lies inside the trust region. The Krylov subspace is terminated such that the termination condition allows the gradient to be decreased on it. A convergence theory of this algorithm is presented. It is shown that this algorithm is globally convergent.

1. Introduction

Nonlinear least squares (NLS) problems are unconstrained optimization problems with special structures. These problems arise in many aspects such as the solution of overdetermined systems of nonlinear equations, some scientific experiments, pattern recognition, and maximum likelihood estimation. For more details about these problems [1]. The general formulation of the NLS problem is to determine the solution 𝑥^⋆∈ℜ^𝑛 that minimizes the function𝑓(𝑥)=‖𝑅(𝑥)‖²₂=^𝑚_𝑖=1(𝑅_𝑖(𝑥))², (1.1)where 𝑅(𝑥)=(𝑅₁(𝑥),𝑅₂(𝑥),…,𝑅_𝑚(𝑥))^𝑡,𝑅_𝑖∶ℜ^𝑛→ℜ,1≤𝑖≤𝑚 and m≥n.

There are two general types of algorithms for solving NLS problem (1.1). The first type is most closely related to solving systems of nonlinear equations [2] and it leads to Gauss-Newton and Levenberg-Marquart methods [3]. The second type is closely related to unconstrained minimization techniques [4]. Of course the two types of algorithms are closely related to each other [5].

The presented algorithm is a Newton-Krylov type algorithm. It requires a fixed-size limited storage proportional to the size of the problem and relies only upon matrix vector product. It is based on the implicitly restarted Arnoldi method (IRAM) to construct a basis for the Krylov subspace and to reduce the Jacobian into a Hessenberg matrix in conjunction with a trust region strategy to control the step on that subspace [6].

Trust region methods for unconstrained minimization are blessed with both strong theoretical convergence properties and a good accurate results in practice. The trial computational step in these methods is to find an approximate minimizer of some model of the true objective function within a trust region for which a suitable norm of the correction lies inside a given bound. This restriction is known as the trust region constraint, and the bound on the norm is its radius. The radius is adjusted so that successive model problems minimized the true objective function within the trust region [7].

The trust region subproblem is the problem of finding 𝑠(Δ) so that𝜙(𝑠(Δ))=min{𝜙(𝑠)∶‖𝑠‖≤Δ}, (1.2)where Δ is some positive constant, ‖⋅‖ is the Euclidean norm in ℜ^𝑛, and𝜙(𝑠)=𝑔^𝑡1𝑠+2𝑠^𝑡𝐻𝑠, (1.3)where 𝑔∈ℜ^𝑛 and 𝐻∈ℜ^𝑛×𝑛 are, respectively, the gradient vector and the Hessian matrix or their approximations.

There are two different approaches to solve (1.2). These approaches are based on either solving several linear systems [8] or approximating the curve 𝑠(Δ) by a piecewise linear approximation “dogleg strategies” [9]. In large-scale optimization, solving linear systems is computationally expensive. Moreover, it is not clear how to define the dogleg curves when the matrix 𝐻 is singular [10].

Several authors have studied inexact Newton's methods for solving NLS problems [11]. Xiaofang et al. have introduced stable factorized quassi-Newton methods for solving large-scale NLS [12]. Dennis et al. proposed a convergence theory for structured BFGS secant method with an application for NLS [13]. The Newton-Krylov method is an example of inexact Newton methods. Krylov techniques inside a Newton iteration in the context of system of equations have been proposed in [14]. The recent work of Sorensen, provides an algorithm which is based on recasting the trust region subproblem into a parameterized eigenvalue problem. This algorithm provides a super linearly convergent scheme to adjust the parameter and find the optimal solution from the eigenvector of the parameterized problem, as long as the hard case does not occur [15].

This contribution is organized as follows. In Section 2, we introduce the structure of the NLS problem and a general view about the theory of the trust region strategies and their convergence. The statement of the algorithm and its properties is developed in Section 3. Section 4 states the assumptions and presents the role of restarting mechanism to control the dimension of the Krylov subspace. The global convergence analysis is presented in Section 5. Concluding remarks and future ideas are given in the last section.

2. Structure of the Problem

The subspace technique plays an important role in solving the NLS problem (1.1). We assume that the current iterate 𝑥_𝑘 and a Krylov subspace 𝑆_𝑘. Denote the dimension of 𝑆_𝑘 to be 𝑘_𝑗 and {𝑣₁^(𝑘),𝑣₂^(𝑘),…,𝑣_𝑘^(𝑘)_𝑗} be a set of linearly independent vectors in 𝑆_𝑘. The next iterate 𝑥_𝑘+1 is such that the increment 𝑥_𝑘+1−𝑥_𝑘∈𝑆_𝑘. Thus, we have ‖𝑅_𝑗(𝑥_𝑘+𝑉_𝑘𝑦)‖=0,for𝑗=1,2,…,𝑚 and 𝑦∈ℜ^𝑘𝑗, where the matrix 𝑉_𝑘=[𝑣₁^(𝑘),𝑣₂^(𝑘),…,𝑣_𝑘^(𝑘)_𝑗]. The second-order Taylor expansion of ‖𝑅(𝑥)‖²₂ at 𝑥_𝑘 is‖𝑅(𝑥_𝑘)+𝐽_𝑘𝑠‖²₂+𝑠^𝑡𝑄_𝑘𝑠, (2.1)where 𝑄_𝑘∈ℜ^𝑛×𝑛 is defined by 𝑄_𝑘=∑^𝑛_𝑗=1𝑅_𝑗(𝑥_𝑘)∇²𝑅_𝑗(𝑥_𝑘) and the jacobian matrix is 𝐽_𝑘=𝐽(𝑥_𝑘)=[∇𝑅₁(𝑥_𝑘),∇𝑅₂(𝑥_𝑘),…,∇𝑅_𝑚(𝑥_𝑘)]^𝑡. If we consider 𝑠∈𝑆_𝑘, we get the quadratic model𝑞_𝑘(𝑦)=‖𝑅(𝑥_𝑘)+𝐽_𝑘𝑉_𝑘𝑦‖²₂+𝑦^𝑡𝐴_𝑘𝑦, (2.2)where 𝐴_𝑘∈ℜ^{𝑘𝑗,𝑘𝑗} approximates the reduced matrix 𝑉^𝑡_𝑘𝑄_𝑘𝑉_𝑘=∑^𝑚_𝑗=1𝑅_𝑗(𝑥_𝑘)𝑉^𝑡_𝑘∇²𝑅_𝑗(𝑥_𝑘)𝑉_𝑘.

The first order necessary conditions of (2.1) is(𝐽^𝑡_𝑘𝐽_𝑘+𝑄_𝑘)𝑠+𝑅^𝑡_𝑘𝐽_𝑘=0. (2.3)Thus, a solution of the quadratic model is a solution to an equation of the form(𝐽^𝑡_𝑘𝐽_𝑘+𝑄_𝑘)𝑠=−𝑅^𝑡_𝑘𝐽_𝑘. (2.4)

The model trust region algorithm generates a sequence of points 𝑥_𝑘, and at the 𝑘th stage of the iteration the quadratic model of problem (1.2) has the form𝜙_𝑘(𝑠)=𝐽^𝑡_𝑘𝑅_𝑘1𝑠+2𝑠^𝑡(𝐽^𝑡_𝑘𝐽_𝑘+𝑄_𝑘)𝑠. (2.5)At this stage, an initial value for the trust region radius Δ_𝑘 is also available. An inner iteration is then performed which consists of using the current trust region radius, Δ_𝑘, and the information contained in the quadratic model to compute a step, 𝑠(Δ_𝑘). Then a comparison of the actual reduction of the objective functionared_𝑘(Δ_𝑘)=‖𝑅(𝑥_𝑘)‖²₂−‖𝑅(𝑥_𝑘+𝑠_𝑘(Δ_𝑘))‖²₂, (2.6)and the reduction predicted by the quadratic modelpred_𝑘(Δ_𝑘)=‖𝑅(𝑥_𝑘)‖²₂−𝑞_𝑘(𝑦) (2.7)is performed. If there is a satisfactory reduction, then the step can be taken, or a possibly larger trust region used. If not, then the trust region is reduced and the inner iteration is repeated. For now, we leave unspecified what algorithm is used to form the step 𝑠(Δ), and how the trust region radius Δ is changed. We also leave unspecified the selection of (𝐽^𝑡_𝑘𝐽_𝑘+𝑄_𝑘) except to restrict it to be symmetric. Details on these points will be addressed in our forth coming paper.

The solution 𝑠 of the quadratic model (2.2) is a solution to an equation of the form[(𝐽^𝑡𝐽+𝑄)+𝜇𝐼]𝑠=−𝑅^𝑡𝐽, (2.8)with 𝜇≥0, 𝜇(𝑠^𝑡𝑠−Δ²)=0 and [(𝐽^𝑡𝐽+𝑄)+𝜇𝐼] is positive semidefinite. This represents the first-order necessary conditions concerning the pair 𝜇 and 𝑠, where 𝑠 is a solution to model (2.2) and 𝜇 is the Lagrange multiplier associated with the constraint 𝑠^𝑡𝑠≤Δ². The sufficient conditions that will ensure 𝑠 to be a solution to (2.2) can be established in the following lemma which has a proof in [16]. Lemma 2.1. Let 𝜇∈ℜ and 𝑠∈ℜ^𝑛 satisfy [(𝐽^𝑡𝐽+𝑄)+𝜇𝐼]𝑠=−𝑅^𝑡𝐽 with [(𝐽^𝑡𝐽+𝑄)+𝜇𝐼] positive semidefinite, if𝜇=0,‖𝑠‖≤Δthen𝑠solves(2.2),‖𝑠‖=Δ,then𝑠solvesmin{𝜙(𝑠)∶‖𝑠‖=Δ},𝜇≥0with‖𝑠‖=Δ,then𝑠solves(2.2). (2.9) Of course, if [(𝐽^𝑡𝐽+𝑄)+𝜇𝐼] is positive definite then 𝑠 is a unique solution.

The main result which is used to prove that the sequence of gradients tends to zero for modified Newton methods is1𝑓−𝜙(𝑠)≥2‖𝑅^𝑡𝐽‖min{Δ,‖𝑅^𝑡𝐽‖‖𝐽^𝑡𝐽+𝑄‖}, (2.10)where 𝑠 is a solution to (2.2). The geometric interpretation of this inequality is that for a quadratic function, any solution 𝑠 to (2.2) produces a decrease 𝑓−𝜙(𝑠) that is at least as much as the decrease along the steepest descent direction −𝑅^𝑡𝐽. A proof of this result may be found in [17].

3. Algorithmic Framework

The algorithm we will discuss here requires that 𝑓 is twice differentiable at any point 𝑥 in the domain of 𝑓. This algorithm involves two levels. In the first level we use IRA method to reduce the Hessian to a tridiagonal matrix and construct an orthonormal basis for the invariant subspace of the Hessian matrix. The second level is used to compute the step and update the trust region radius for the reduced local model (a model defined on the subspace). The desired properties of this algorithm are:

(1)the algorithm should be well defined for a sufficiently general class of functions and it is globally convergent;(2)the algorithm should be invariant under linear affine scalings of the variables, that is, if we replace 𝑓(𝑥) by 𝑓(𝑦)=𝑓(𝑠+𝑉𝑦), where 𝑉∈ℜ^𝑛×𝑘 is an orthonormal matrix, 𝑥∈ℜ^𝑛 and 𝑦∈ℜ^𝑘, then applying the iteration to 𝑓 with initial guess 𝑦₀ satisfying 𝑥₀=𝑠+𝑉𝑦₀ should produce a sequence {𝑦_𝑙} related to the sequence {𝑥_𝑙} by 𝑥_𝑙=𝑉𝑦_𝑙+𝑠, where 𝑥_𝑙 is produced by applying the algorithm to 𝑓 with initial guess 𝑥₀;(3) the algorithm should provide a decrease that is at least as large as a given multiple of the minimum decrease that would be provided by a quadratic search along the steepest descent direction;(4)the algorithm should give as good a decrease of the quadratic model as a direction of the negative gradient when the Hessian, (𝐽^𝑡𝐽+𝑄), is indefinite and should force the direction to be equal to the Newton direction when (𝐽^𝑡𝐽+𝑄) is symmetric positive definite.

The following describes a full iteration of a truncated Newton type method. Some of the previous characteristics will be obvious and the other ones will be proved in the next section.

Algorithm 3.1. (1)Step 0 (initialization).Given 𝑥₀∈ℜ^𝑛, compute (𝑅^𝑡𝐽)(𝑥₀), (𝐽^𝑡𝐽+𝑄)(𝑥₀).Choose 𝜉₁<𝜉₂, 𝜂₁∈(0,1),𝜂₂>1, Δ₁>0,𝑘=1.(2)Step 1 (construction a basis for the Krylov subspace). (a)Choose 𝜖∈(0,1), an initial guess 𝑠₀ which can be chosen to be zero, form 𝑟₀=−𝑅^𝑡𝐽−(𝐽^𝑡𝐽+𝑄)𝑠₀, where 𝑅^𝑡𝐽=(𝑅^𝑡𝐽)(𝑥_𝑛),(𝐽^𝑡𝐽+𝑄)=(𝐽^𝑡𝐽+𝑄)(𝑥_𝑛), and compute 𝜏=‖𝑟₀‖ and 𝑣₁=𝑟₀/𝜏.(b)For 𝑗=1,2,… until convergence, do.(i)Form (𝐽^𝑡𝐽+𝑄)𝑣_𝑗 and orthogonalize it against the previous 𝑣₁,𝑣₂,…,𝑣_𝑗 thenℎ(𝑖,𝑗)=((𝐽^𝑡𝐽+𝑄)𝑣_𝑗,𝑣_𝑖),𝑖=1,2,…,𝑗,𝑣_𝑗+1=(𝐽^𝑡𝐽+𝑄)𝑣_𝑗−∑^𝑗_𝑖=1ℎ(𝑖,𝑗)𝑣_𝑖,ℎ(𝑖+1,𝑗)=‖𝑣_𝑗+1‖.𝑣_𝑗+1=𝑣_𝑗+1/ℎ(𝑖+1,𝑗).(ii)Compute the residual norm 𝛽_𝑗=‖𝑅^𝑡𝐽+(𝐽^𝑡𝐽+𝑄)𝑠_𝑗‖ of the solution 𝑠_𝑗 that would be obtained if we stopped.(iii)If 𝛽_𝑗≤𝜖_𝑗‖𝑅^𝑡𝐽‖, where 𝜖_𝑗≤(1/𝒟_𝑗)∈(0,𝜖], and 𝒟_𝑗 is an approximation to the condition number of (𝐽^𝑡𝐽+𝑄)_𝑗, then set 𝑘=𝑗 and go to 3, else continue.(3)Step 2 (compute the trial step). (i)Construct the local model,𝜓_𝑘(𝑦)=𝑓_𝑘+ℎ^𝑡_𝑘1𝑦+2𝑦^𝑡𝐻_𝑘𝑦, (3.1)where ℎ_𝑘=𝑉^𝑡_𝑘(𝑅^𝑡𝐽)_𝑘∈ℜ^𝑘𝑗 and 𝐻_𝑘=𝑉^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑘𝑉_𝑘∈ℜ^{𝑘𝑗×𝑘𝑗}.(ii)Compute the solution 𝑦 to the problem min{𝜓_𝑘(𝑦)∶‖𝑦‖≤Δ_𝑘} and set 𝑠_𝑘=𝑉_𝑘𝑦.(4)Step 3 (test the step and update Δ_𝑘). (i)Evaluate 𝑓_𝑘+1=𝑓(𝑥_𝑘+𝑉_𝑘𝑦), ared_𝑘=𝑓_𝑘−𝑓(𝑥_𝑘+𝑉_𝑘𝑦) and pred_𝑘=𝑓_𝑘−𝜓(𝑉_𝑘𝑦).(ii)If (ared_𝑘/pred_𝑘)<𝜉₁ thenbegin Δ_𝑘=𝜂₁Δ_𝑘 go to 3(i) Comment: Do not accept the step “‖𝑦‖”, reduce the trust region radius “Δ_𝑘” and compute a new stepElse if 𝜉₁≤(ared_𝑘/pred_𝑘)≤𝜉₂ then𝑥_𝑘+1=𝑥_𝑘+𝑉_𝑘𝑦Comment: Accept the step and keep the previous trust region radius the sameElse (ared_𝑘/pred_𝑘)>𝜉₂ thenΔ_𝑘=𝜂₂Δ_𝑘Comment: Accept the step and increase the trust region radiusend if(5)Step 4.Set 𝑘+1←𝑘 and go to Step 1.

4. The Restarting Mechanism

In this section, we discuss the important role of the restarting mechanism to control the possibility of the failure of nonsingularity of the Hessian matrix, and introduce the assumptions under which we prove the global convergence.

Let the sequence of the iterates generated by Algorithm 3.1 be {𝑥_𝑘}; for such sequence we assume

(G1)for all 𝑘,𝑥_𝑘 and 𝑥_𝑘+𝑠_𝑘∈Ω where Ω⊂ℜ^𝑛 is a convex set;(G2)𝑓∈𝐶²(Ω) and (𝑅^𝑡𝐽)=∇𝑓(𝑥)≠0 for all 𝑥∈Ω;(G3)𝑓(𝑥) is bounded in norm in Ω and (𝐽^𝑡𝐽+𝑄) is nonsingular for all 𝑥∈Ω;(G4)for all 𝑠_𝑘 such that 𝑥_𝑘+𝑠_𝑘∈Ω, the termination condition, ‖(𝑅^𝑡𝐽)_𝑘+(𝐽^𝑇𝐽+𝑄)𝑠_𝑘‖≤𝜖_𝑗‖(𝑅^𝑡𝐽)_𝑘‖, 𝜖_𝑘∈(0,1) is satisfied;

An immediate consequence of the global assumptions is that the Hessian matrix (𝐽^𝑡𝐽+𝑄)(𝑥) is bounded, that is, there exists a constant 𝛽₁>0 such that ‖(𝐽^𝑡𝐽+𝑄)(𝑥)‖<𝛽₁. Therefore, the condition numbers of the Hessian, 𝒟_𝑘=cond₂((𝐽^𝑡𝐽+𝑄)_𝑘) are bounded.

Assumption G3 that (𝐽^𝑡𝐽+𝑄) is nonsingular is necessary since the IRA iteration is only guaranteed to converge for nonsingular systems. We first discuss the two ways in which the IRA method can breakdown. This can happen if either 𝑣_𝑘+1=0 so that 𝑣_𝑘+1 cannot be formed, or if 𝐻_𝑘max is singular which means that the maximum number of IRA steps has been taken, but the final iterate cannot be formed. The first situation has often been referred to as a soft breakdown, since 𝑣_𝑘+1=0 implies (𝐽^𝑡𝐽+𝑄)_𝑘 is nonsingular and 𝑥_𝑘 is the solution [18]. The second case is more serious in that it causes a convergence failure of the Arnoldi process. One possible recourse is to hope that (𝐽^𝑡𝐽+𝑄)_𝑘 is nonsingular for some 𝑘 among 1,2,…,𝑘_max−1. If such 𝑘 exists, we can compute 𝑥_𝑗 and then restart the algorithm using 𝑥_𝑘 as the new initial guess 𝑥₀. It may not be possible to do this.

To handle the problem of singular 𝐻_𝑘, we can use a technique similar to that done in the full dimensional space. If 𝐻_𝑘 is singular, or its condition number is greater than 1/√𝜈, where 𝜈 is the machine epsilon, then we perturb the quadratic model 𝜓(𝑦) to1𝜓(𝑦)=𝜓(𝑦)+2𝛾𝑦^𝑡𝑦=𝑓+ℎ^𝑡1𝑦+2𝑦^𝑡(𝐻_𝑘+𝛾𝐼)𝑦, (4.1)where 𝛾=√𝑛𝜈‖𝐻_𝑘‖₁. The condition number of (𝐻_𝑘+𝛾𝐼) is roughly 1/√𝜈. For more details about this technique see reference [4].

We have discussed the possibility of stagnation in the linear IRA method which results in a break down in the nonlinear iteration. Sufficient conditions under which stagnation of the linear iteration never occurs are

(1)we have to ensure that the steepest descent direction belongs to the subspace, 𝒮_𝑘. Indeed in Algorithm 3.1, the Krylov subspace will contain the steepest descent direction because 𝑣₁=−𝑅^𝑡𝐽/‖𝑅^𝑡𝐽‖ and the Hessian is nonsingular;(2)there is no difficulty, if one required the Hessian matrix, (𝐽^𝑡𝐽+𝑄)(𝑥) to be positive definite for all 𝑥∈Ω, then the termination condition can be satisfied, and so convergence of the sequence of iterates will follow. This will be the case for some problems, but, requiring (𝐽^𝑡𝐽+𝑄)(𝑥) to be positive definite every where is very restrictive.

One of the main restrictions of most of the Newton-Krylov schemes is that the subspace onto which a given Newton step is projected must solve the Newton equations with a certain accuracy which is monitored by the termination condition (assumption G4). This condition is enough to essentially guarantee convergence of the trust region algorithm. Practically, the main difficulty is that one does not know in advance if the subspace chosen for projection will be good enough to guarantee this condition. Thus, 𝑘_𝑗 can be chosen as large as the termination condition will eventually be satisfied, but, when 𝑘_𝑗 is too large the computational cost and the storage become too high. An alternative to that is to restart the algorithm keeping (𝐽^𝑡𝐽+𝑄) nonsingular and 𝑘_𝑗<𝑛. Moreover, preconditioning and scaling are essential for the successful application of these schemes [19].

5. Global Convergence Analysis

In this section, we are going to establish some convergence properties which are possessed by Algorithm 3.1. The major differences between the proposed results and the preexisting ones arise from the fact that a lower dimensional quadratic model is used, rather than the full dimension model that is used in the preexisting results [10].Lemma 5.1. Let the global assumptions hold. Then for any 𝑠∈Ω, one has|(𝑅^𝑡𝐽)^𝑡_𝑘𝑠|‖𝑠‖₂≥𝒟_𝑘−𝜖_𝑘‖(𝐽^𝑡𝐽+𝑄)_𝑘‖(1+𝜖_𝑘)𝒟_𝑘‖(𝑅^𝑡𝐽)_𝑘‖>0. (5.1)Proof. Suppose 𝑟_𝑘 be the residual associated with 𝑠 so that 𝑟_𝑘=(𝑅^𝑡𝐽)_𝑘+(𝐽^𝑡𝐽+𝑄)_𝑘𝑠 and (𝑅^𝑡𝐽)_𝑘≠0. Then ‖𝑟_𝑘‖≤𝜖_𝑘‖(𝑅^𝑡𝐽)_𝑘‖, and 𝑠=−(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹((𝑅^𝑡𝐽)_𝑘−𝑟_𝑘). So,(𝑅^𝑡𝐽)^𝑡_𝑘𝑠=(𝑅^𝑡𝐽)^𝑡_𝑘[−(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹((𝑅^𝑡𝐽)_𝑘−𝑟_𝑘)]=−(𝑅^𝑡𝐽)^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹(𝑅^𝑡𝐽)_𝑘+(𝑅^𝑡𝐽)^𝑡_𝑘𝑟_𝑘. (5.2)Hence,|(𝑅^𝑡𝐽)^𝑡_𝑘𝑠|=‖𝑠‖|(𝑅^𝑡𝐽)^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹(𝑅^𝑡𝐽)_𝑘−(𝑅^𝑡𝐽)^𝑡_𝑘𝑟_𝑘|‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹((𝑅^𝑡𝐽)_𝑘−𝑟_𝑘≥)‖(𝑅^𝑡𝐽)^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹(𝑅^𝑡𝐽)_𝑘−|(𝑅^𝑡𝐽)^𝑡_𝑘𝑟_𝑘|‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹((𝑅^𝑡𝐽)_𝑘−𝑟_𝑘.)‖ (5.3)Next, ‖𝑟_𝑘‖≤𝜖‖(𝑅^𝑡𝐽)_𝑘‖ implies |(𝑅^𝑡𝐽)^𝑡_𝑘𝑟_𝑘𝑗|≤𝜖_𝑘‖(𝑅^𝑡𝐽)_𝑘‖² which gives(𝑅^𝑡𝐽)^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹(𝑅^𝑡𝐽)_𝑘−|(𝑅^𝑡𝐽)^𝑡_𝑘𝑟_𝑘|≥(‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹‖−𝜖_𝑘)‖(𝑅^𝑡𝐽)_𝑘‖²,‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹((𝑅^𝑡𝐽)_𝑘−𝑟_𝑘)‖≤‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹‖‖(𝑅^𝑡𝐽)_𝑘‖+‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹𝑟_𝑘‖≤(1+𝜖_𝑘)‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹‖‖(𝑅^𝑡𝐽)_𝑘‖. (5.4)Introduce (5.4) into (5.3), we obtain the following:|(𝑅^𝑡𝐽)^𝑡_𝑘𝑠|≥‖𝑠‖(‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹‖−𝜖_𝑘)‖(𝑅^𝑡𝐽)_𝑘‖²(1+𝜖_𝑘)‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹‖‖(𝑅^𝑡𝐽)_𝑘‖≥(‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹‖−𝜖_𝑘)(1+𝜖_𝑘)‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹‖‖(𝑅^𝑡𝐽)_𝑘‖𝒟=(^𝑘−𝜖_𝑘‖(𝐽^𝑡𝐽+𝑄)_𝑘‖(1+𝜖_𝑘)𝒟_𝑘)‖(𝑅^𝑡𝐽)_𝑘‖. (5.5)Inequality (5.1) can be written as|(𝑅^𝑡𝐽)^𝑡_𝑘𝑠|‖𝑠‖‖(𝑅^𝑡𝐽)_𝑘‖≥𝒟_𝑘−𝜖_𝑘‖(𝐽^𝑡𝐽+𝑄)_𝑘‖(1+𝜖_𝑘)𝒟_𝑘. (5.6)This condition (5.6) tells us at each iteration of Algorithm 3.1, we require that the termination condition holds. Since the condition numbers, 𝒟_𝑘, are bounded from above, condition (5.6) givescos𝜃≥𝒟−𝜖‖(𝐽^𝑡𝐽+𝑄)_𝑘‖(1+𝜖)𝒟, (5.7)where 𝜃 is the angle between (𝑅^𝑡𝐽)_𝑘 and 𝑠 and 𝒟 are the upper bound of the condition numbers 𝒟_𝑘. Inequality (5.7) shows that the acute angle 𝜃 is bounded away from 𝜋/2.

The following lemma shows that the termination norm assumption G4 implies that the cosine of the angle between the gradient and the Krylov subspace is bounded below.Lemma 5.2. Let the global assumptions hold. Then‖𝑉^𝑡_𝑘(𝑅^𝑡𝐽)_𝑘𝒟‖≥^𝑘−𝜖_𝑘‖(𝐽^𝑡𝐽+𝑄)_𝑘‖(1+𝜖_𝑘)𝒟_𝑘‖(𝑅^𝑡𝐽)_𝑘‖. (5.8)Proof. Suppose 𝑠=𝑉_𝑘𝑦 be a vector of 𝒮_𝑘 satisfying the termination condition, ‖(𝑅^𝑡𝐽)_𝑘+(𝐽^𝑡𝐽+𝑄)_𝑘𝑠‖≤𝜖_𝑘‖(𝑅^𝑡𝐽)_𝑘‖,𝜖_𝑘∈(0,1). From Lemma 5.1 and the fact that ‖𝑠‖=‖𝑉𝑦‖=‖𝑦‖, we obtain|(𝑅^𝑡𝐽)^𝑡_𝑘𝑉_𝑘𝑦|≥𝒟‖𝑦‖^𝑘−𝜖_𝑘‖(𝐽^𝑡𝐽+𝑄)_𝑘‖(1+𝜖_𝑘)𝒟_𝑘‖(𝑅^𝑡𝐽)_𝑘‖. (5.9)From Cauchy-Schwartz inequality, we obtain|(𝑅^𝑡𝐽)^𝑡_𝑘𝑉_𝑘𝑦|≤‖𝑉^𝑡_𝑘(𝑅^𝑡𝐽)_𝑘‖‖𝑦‖. (5.10)Combining formulas (5.9) and (5.10), we obtain the result.

The following Lemma establishes that Algorithm 3.1 converges with a decrease in the quadratic model on the lower dimensional space at least equal to the decrease in the quadratic model on the full dimensional space.Lemma 5.3. Let the global assumptions hold and let 𝑦 be a solution of the minimization problem, min_{‖𝑦‖≤Δ}𝜓(𝑦). Then‖ℎ_𝑘‖min{Δ_𝑘,‖ℎ_𝑘‖‖𝐻_𝑘‖}≥𝛼_𝑘‖(𝑅^𝑡𝐽)_𝑘‖min{Δ_𝑘,𝛼_𝑘‖(𝑅^𝑡𝐽)_𝑘‖‖(𝐽^𝑡𝐽+𝑄)_𝑘‖}, (5.11)where 𝛼_𝑘=(𝒟_𝑘−𝜖_𝑘‖(𝐽^𝑡𝐽+𝑄)_𝑘‖)/(1+𝜖_𝑘)𝒟_𝑘.Proof. Since, 𝜓(𝑦)=𝜙(𝑉_𝑘𝑦)=𝑓_𝑘+(𝑅^𝑡𝐽)^𝑡_𝑘𝑉_𝑘𝑦+(1/2)(𝑉_𝑘𝑦)^𝑡(𝐽^𝑡𝐽+𝑄)(𝑉_𝑘𝑦). Thus ∇𝜓(0)=𝑉^𝑡_𝑘(𝑅^𝑡𝐽)_𝑘. From Step 3 of Algorithm 3.1, we obtain𝑓_𝑘−𝑓_𝑘+1≥𝜉₁(𝑓_𝑘−𝜙_𝑘≥𝜉(𝑦))₁2‖𝑉^𝑡_𝑘(𝑅^𝑡𝐽)_𝑘‖min{Δ_𝑘,‖𝑉^𝑡_𝑘(𝑅^𝑡𝐽)_𝑘‖‖𝑉^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑘𝑉_𝑘‖}=𝜉₁2‖ℎ_𝑘‖min{Δ_𝑘,‖ℎ_𝑘‖‖𝐻_𝑘‖}. (5.12)(5.13)(5.14)We have to convert the lower bound of this inequality to one involving ‖(𝑅^𝑡𝐽)_𝑘‖ rather than ‖𝑉^𝑡_𝑘(𝑅^𝑡𝐽)_𝑘‖. Using Lemma 5.2, we obtain ‖ℎ_𝑘‖≥𝛼_𝑘‖(𝑅^𝑡𝐽)_𝑘‖, but, ‖𝐻_𝑘‖=‖𝑉^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑘𝑉_𝑘‖≤‖(𝐽^𝑡𝐽+𝑄)_𝑘‖. Substituting in (5.14), we obtain‖ℎ_𝑘‖2min{Δ_𝑘,‖ℎ_𝑘‖‖𝐻_𝑘‖𝛼}≥^𝑘2‖(𝑅^𝑡𝐽)_𝑘‖min{Δ_𝑘,𝛼_𝑘‖(𝑅^𝑡𝐽)_𝑘‖‖(𝐽^𝑡𝐽+𝑄)_𝑘‖}, (5.15)which completes the proof.

The following two facts will be used in the remainder of the proof.

Fact 1.
By Taylor's theorem for any 𝑘 and Δ>0, we have|ared_𝑘(Δ)−pred_𝑘(Δ)|=|(𝑓_𝑘−𝑓_𝑘+1)−(−ℎ^𝑡_𝑘𝑦_𝑘1(Δ)−2𝑦^𝑡_𝑘(Δ)𝑉^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑘𝑉_𝑘𝑦_𝑘1(Δ))|=|2𝑦^𝑡_𝑘(Δ)𝐻_𝑘𝑦_𝑘(Δ)−¹₀𝑦^𝑡_𝑘(Δ)𝐻_𝑘(𝜂)𝑦_𝑘(Δ)(1−𝜂)𝑑𝜂|≤‖𝑦_𝑘(Δ)‖²¹₀‖𝐻_𝑘−𝐻_𝑘(𝜂)‖(1−𝜂)𝑑𝜂. (5.16)where, 𝐻_𝑘(𝜂)=𝑉^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)(𝑥_𝑘+𝜂𝑠_𝑘)𝑉_𝑘 with 𝑠_𝑘=𝑥_𝑘+1−𝑥_𝑘=𝑉𝑦_𝑘. Thus,|ared_𝑘(Δ)pred_𝑘(Δ)−1|≤‖𝑦_𝑘(Δ)‖²|pred_𝑘(Δ)|¹₀‖𝐻_𝑘−𝐻_𝑘(𝜂)‖(1−𝜂)𝑑𝜂. (5.17)

Fact 2.
For any sequence {𝑥_𝑘} generated by an algorithm satisfying the global assumptions, the related sequence {𝑓_𝑘} is monotonically decreasing and bounded from below, that is, 𝑓_𝑘→𝑓_⋆ as 𝑘→∞.

The next result establishes that every limit point of the sequence {𝑥_𝑘} satisfies the first-order necessary conditions for a minimum.Lemma 5.4. Let the global assumptions hold and Algorithm 3.1 is applied to 𝑓(𝑥), generating a sequence {𝑥_𝑘},𝑥_𝑘∈Ω,and𝑘=1,2,…, then ℎ_𝑘→0.Proof. Since 𝜖_𝑘≤𝜖<1 for all 𝑘, we have𝛼_𝑘≥𝒟−𝜖‖(𝐽^𝑡𝐽+𝑄)_𝑘‖(1+𝜖)𝒟=𝛼>0∀𝑘. (5.18)Consider any 𝑥_𝑙 with ‖(𝑅^𝑡𝐽)_𝑙‖≠0, ‖(𝑅^𝑡𝐽)(𝑥)−(𝑅^𝑡𝐽)_𝑙‖≤𝛽₁‖𝑥−𝑥_𝑙‖. Thus, if ‖𝑥−𝑥_𝑘‖≤(‖(𝑅^𝑡𝐽)_𝑘‖/2𝛽₁), then ‖(𝑅^𝑡𝐽)(𝑥)‖≥‖(𝑅^𝑡𝐽)_𝑘(𝑅^𝑡𝐽)_𝑙‖−‖(𝑅^𝑡𝐽)(𝑥)−(𝑅^𝑡𝐽)_𝑘‖≥(‖(𝑅^𝑡𝐽)_𝑘‖/2). Let, 𝜌=‖(𝑅^𝑡𝐽)_𝑘‖/2 and ℬ_𝜌={𝑥∶‖𝑥−𝑥_𝑘‖<𝜌}. There are two cases, either for all 𝑘≥𝑙,𝑥_𝑙∈ℬ_𝜌 or eventually {𝑥_𝑙} leaves the ball ℬ_𝜌. Suppose 𝑥_𝑙∈ℬ_𝜌 for all 𝑘≥𝑙, then ‖𝑔_𝑙‖≥‖𝑔_𝑘‖/2=𝜎 for all 𝑙>𝑘. Thus by Lemma 5.3, we havepred_𝑘𝜉(Δ)≥¹𝛼𝜎2min{Δ_𝑘,𝛼𝜎‖𝐻_𝑘‖}∀𝑘≥𝑙. (5.19)Put 𝛼𝜎=𝜎, introduce inequality (5.19) into (5.17) and since ‖𝐻_𝑘‖≤‖(𝐽^𝑡𝐽+𝑄)_𝑘‖≤𝛽₁, we obtain|ared_𝑘(Δ)pred_𝑘(Δ)−1|≤‖𝑦_𝑘‖²∫¹₀‖𝐻_𝑘−𝐻_𝑘(𝜂)‖(1−𝜂)𝑑𝜂𝜉₁𝜎min{Δ_𝑘,𝜎/𝛽₁}≤Δ²_𝑘(2𝛽₁)𝜉₁𝜎min{Δ_𝑘,𝜎/𝛽₁}≤2Δ_𝑘𝛽₁𝜉₁𝜎,Δ_𝑘≤𝜎𝛽₁∀𝑘≥𝑙. (5.20)This gives for Δ_𝑘 sufficiently small and 𝑘≥𝑙 thatared_𝑘(Δ)pred_𝑘(Δ)>𝜉₂. (5.21)In addition, we have‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹ℎ_𝑘‖≥‖ℎ_𝑘‖‖𝐻_𝑘‖≥𝛼‖(𝑅^𝑡𝐽)_𝑘‖𝛽₁≥𝜎𝛽₁. (5.22)Inequalities (5.21) and (5.22) tell us that, for Δ_𝑘 sufficiently small, we cannot get a decrease in Δ_𝑘 in Step 3(ii) of Algorithm 3.1. It follows that Δ_𝑘 is bounded away from 0, but, since𝑓_𝑘−𝑓_𝑘+1=ared_𝑘(Δ)≥𝜉₁pred_𝑘(Δ)≥𝜉₁𝜎min{Δ_𝑘,𝜎𝛽₁}, (5.23)and since 𝑓 is bounded from below, we must have Δ_𝑘→0 which is a contradiction. Therefore, {𝑥_𝑘} must be outside ℬ_𝜌 for some 𝑘>𝑙.
Let 𝑥_𝑗+1 be the first term after 𝑥_𝑙 that does not belong to ℬ_𝜌. From inequality (5.23), we get𝑓(𝑥_𝑙)−𝑓(𝑥_𝑗+1)=^𝑗_𝑘=𝑙(𝑓_𝑘−𝑓_𝑘+1)≥^𝑙_𝑗+1𝜉₁pred_𝑘≥(Δ)_{𝑗𝑘=𝑙}𝜉₁𝜎min{Δ_𝑘,𝜎𝛽₁}. (5.24)If Δ_𝑘≤𝜎/𝛽₁ for all 𝑙≤𝑘≤𝑗, we have𝑓(𝑥_𝑙)−𝑓(𝑥_𝑗+1)≥𝜉₁𝜎^𝑗_𝑘=𝑙Δ_𝑘≥𝜉₁𝜎𝜌. (5.25)If there exists at least one 𝑘 with Δ_𝑘>𝜎/𝛽₁, we get𝑓(𝑥_𝑚)−𝑓(𝑥_𝑚+1)≥𝜉₁𝜎²𝛽₁. (5.26)In either case, we have𝑓(𝑥_𝑚)−𝑓(𝑥_𝑗+1)≥𝜉₁𝜎𝜎min{𝜌,𝛽₁}=𝜉₁𝛼‖(𝑅^𝑡𝐽)_𝑘‖2min{𝛼‖(𝑅^𝑡𝐽)_𝑘‖2𝛽₁,𝛼‖(𝑅^𝑡𝐽)_𝑘‖2𝛽₁}≥𝜉𝛼²4𝛽₁‖(𝑅^𝑡𝐽)_𝑘‖². (5.27)Since, 𝑓 is bounded below, {𝑓_𝑘} is a monotonically decreasing sequence (i.e.,𝑓_𝑘→𝑓_⋆). Hence,‖(𝑅^𝑡𝐽)_𝑘‖²≤4𝛽₁𝜉₁𝛼²(𝑓_𝑘−𝑓_⋆), (5.28)(i.e.,(𝑅^𝑡𝐽)_𝑘→0 as 𝑘→∞). The proof is completed by using Lemma 5.2.

The following lemma proves that under the global assumptions, if each member of the sequence of iterates generated by Algorithm 3.1 satisfies the termination condition of the algorithm, then this sequence converges to one of its limit points.Lemma 5.5. Let the global assumptions hold, 𝐻_𝑘=𝑉^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)(𝑥_𝑘)𝑉_𝑘 for all 𝑘, (𝐽^𝑡𝐽+𝑄)(𝑥) is Lipschitz continuous with 𝐿, and 𝑥_∗ is a limit point of the sequence {𝑥_𝑘} with (𝐽^𝑡𝐽+𝑄)(𝑥_∗) be positive definite. Then {𝑥_𝑘} converges q-Superlinearly to 𝑥_⋆.

Proof. Let {𝑥_𝑘𝑚} be a subsequence of the sequence {𝑥_𝑘} that converges to 𝑥_⋆, where, 𝑥_⋆ is a limit of {𝑥_𝑘}. From Lemma 5.4, (𝑅^𝑡𝐽)(𝑥_⋆)=0. Since (𝐽^𝑡𝐽+𝑄)(𝑥_⋆) is positive definite and (𝐽^𝑡𝐽+𝑄)(𝑥) is continuous there exists 𝛿₁>0 such that if ‖𝑥−𝑥_⋆‖<𝛿₁, then (𝐽^𝑡𝐽+𝑄)(𝑥) is positive definite, and if 𝑥≠𝑥_⋆ then (𝑅^𝑡𝐽)(𝑥_⋆)≠0. Let, 𝐵₁={𝑥∶‖𝑥−𝑥_⋆‖<𝛿₁}.
Since (𝑅^𝑡𝐽)(𝑥_⋆)=0, we can find 𝛿₂>0 with 𝛿₂<𝛿₁/4 and ‖(𝐽^𝑡𝐽+𝑄)⁻¹‖‖(𝑅^𝑡𝐽)‖≤𝛿₁/2 for all 𝑥∈𝐵₂={𝑥∶‖𝑥−𝑥_⋆‖<𝛿₂}.
Find 𝑙_𝑜 such that 𝑓(𝑥_𝑙𝑜)<inf{𝑓(𝑥)∶𝑥∈𝐵₁−𝐵₂} and 𝑥_𝑙𝑜∈𝐵₂. Consider any 𝑥_𝑖 with 𝑖≥𝑘_𝑙,𝑥_𝑖∈𝐵₂.
We claim that 𝑥_𝑖+1∈𝐵₂, which implies that the entire sequence beyond 𝑥_𝑘𝑙𝑜∈𝐵₂. Suppose 𝑥_𝑖+1 does not belong to 𝐵₂. Since 𝑓_𝑖+1<𝑓_𝑘𝑙𝑜, 𝑥_𝑖+1 does not belong to 𝐵₁ either. SoΔ_𝑖≥‖𝑥_𝑖+1−𝑥_𝑖‖≥‖𝑥_𝑖+1−𝑥_⋆‖−‖𝑥_𝑖−𝑥_⋆‖≥𝛿₁−𝛿₁4=3𝛿₁4>𝛿₁2≥‖(𝐽^𝑡𝐽+𝑄)_𝑖⁻¹‖‖(𝑅^𝑡𝐽)_𝑖‖≥‖𝐻_𝑖⁻¹‖‖ℎ_𝑖‖. (5.29) So, the truncated-Newton step is within the trust region, we obtain 𝑦_𝑖(Δ_𝑖)=−(𝑉^𝑡_𝑖(𝐽^𝑡𝐽+𝑄)_𝑖𝑉_𝑖)⁻¹𝑉^𝑡_𝑖(𝑅^𝑡𝐽)_𝑖. (5.30)Since ‖𝑦_𝑖(Δ)‖<𝛿₁/2, it follows that 𝑥_𝑖+1∈𝐵₁, which is a contradiction. Hence, for all 𝑘≥𝑘_𝑙𝑜,𝑥_𝑘∈𝐵₂, and so since 𝑓(𝑥_𝑘) is strictly decreasing sequence and 𝑥_⋆ is the unique minimizer of 𝑓 in 𝐵₂, we have that {𝑥_𝑘}→𝑥_⋆.

The following lemma establishes the rate of convergence of the sequence of iterates generated by Algorithm 3.1 to be 𝑞-super linear.Lemma 5.6. Let the assumptions of Lemma 5.5 hold. Then the rate of convergence of the sequence {𝑥_𝑘} is 𝑞-Super linear.

Proof. To show that the convergence rate is super linear, we will show eventually that ‖𝐻_𝑘⁻¹ℎ_𝑘‖ will always be less than Δ_𝑗 and hence the truncated-Newton step will always be taken. Since (𝐽^𝑡𝐽+𝑄)(𝑥_⋆) is positive definite, it follows that the convergence rate of {𝑥_𝑘} to 𝑥_⋆ is super linear.
To show that eventually the truncated-Newton step is always shorter than the trust region radius, we need a particular lower bound on pred_𝑘(Δ). By assumptions, for all 𝑘 large enough, 𝐻_𝑘 is positive definite. Hence, either the truncated-Newton step is longer than the trust region radius, or 𝑦_𝑘(Δ) is the truncated-Newton step. In either case‖𝑦_𝑘(Δ)‖≤‖(𝑉^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑗𝑉_𝑗)⁻¹𝑉^𝑡_𝑘(𝑅^𝑡𝐽)_𝑘‖≤‖(𝑉^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑘𝑉_𝑘)⁻¹‖‖𝑉^𝑡_𝑘(𝑅^𝑡𝐽)_𝑘‖, (5.31)and so it follows that ‖𝑉^𝑡_𝑘(𝑅^𝑡𝐽)_𝑘‖≥‖𝑠_𝑘(Δ)‖/‖(𝑉^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑘𝑉_𝑘)⁻¹‖. By Lemma 5.3, for all 𝑘 large enough, we havepred_𝑘𝜉(Δ)≥¹𝛼2‖𝑦_𝑘(Δ)‖‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹‖min{Δ_𝑘,𝛼‖(𝑅^𝑡𝐽)_𝑘‖‖(𝐽^𝑡𝐽+𝑄)_𝑘‖}. (5.32)Thus,pred_𝑘𝜉(Δ)≥¹𝛼2‖𝑦_𝑘(Δ)‖‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹‖min{‖𝑦_𝑘(Δ)‖,𝛼‖𝑦_𝑘(Δ)‖‖(𝐽^𝑡𝐽+𝑄)_𝑘‖‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹‖}≥𝜉₁𝛼²2‖𝑦_𝑘(Δ)‖²𝒟‖(𝐽^𝑡𝐽+𝑄)_𝑘⁻¹‖. (5.33)So, by the continuity of (𝐽^𝑡𝐽+𝑄)(𝑥), for all 𝑘 large enough, we obtainpred_𝑘𝜉(Δ)≥¹𝛼²4‖𝑦_𝑘(Δ)‖²‖(𝐽^𝑡𝐽+𝑄)⁻¹(𝑥_⋆)‖. (5.34)Finally, by the argument leading up to Fact 1 and the Lipschitz continuity assumption, we have|ared_𝑘(Δ)−pred_𝑘(Δ)|≤‖𝑦_𝑘(Δ)‖³𝐿2. (5.35)Thus, for any Δ_𝑘>0 and 𝑘 large enough, we obtain|ared_𝑘(Δ)pred_𝑘(Δ)−1|≤‖𝑠_𝑘(Δ)‖³𝐿2.4‖(𝐽^𝑡𝐽+𝑄)⁻¹(𝑥_⋆)‖𝜉₁𝛼²‖𝑦_𝑘(Δ)‖²=2𝐿𝜉₁𝛼²‖(𝐽^𝑡𝐽+𝑄)⁻¹(𝑥_⋆)‖‖𝑦_𝑘≤(Δ)‖2𝐿‖(𝐽^𝑡𝐽+𝑄)⁻¹(𝑥_⋆)‖𝜉₁𝛼²Δ_𝑘. (5.36)Thus by Step 3(ii) of Algorithm 3.1, there is a Δ such that if Δ_𝑘−1<Δ, then Δ_𝑘 will be less than Δ_𝑘−1 only if Δ_𝑘≥‖(𝑉^𝑡_𝑘−1(𝐽^𝑡𝐽+𝑄)_𝑘−1𝑉_𝑘−1)⁻¹𝑉^𝑡_𝑘−1(𝑅^𝑡𝐽)_𝑘−1‖. It follows from the superlinear convergence of the truncated-Newton method that for 𝑥_𝑘−1 close enough to 𝑥_⋆ and 𝑘 large enough,‖(𝑉^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑘𝑉_𝑘)⁻¹𝑉^𝑡_𝑘(𝑅^𝑡𝐽)_𝑘‖<‖(𝑉^𝑡_𝑘−1(𝐽^𝑡𝐽+𝑄)_𝑘−1𝑉_𝑘−1)⁻¹𝑉^𝑡_𝑘−1(𝑅^𝑡𝐽)_𝑘−1‖. (5.37)Now, if Δ_𝑘 is bounded away from 0 for all large 𝑘, then we are done. Otherwise, if for an arbitrarily large 𝑘, Δ_𝑘 is reduced, that is, Δ_𝑘<Δ_𝑘−1, then we haveΔ_𝑘≥‖(𝑉^𝑡_𝑘−1(𝐽^𝑡𝐽+𝑄)_𝑘−1𝑉_𝑘−1)⁻¹𝑉^𝑡_𝑘−1(𝑅^𝑡𝐽)_𝑘−1‖>‖(𝑉^𝑡_𝑘(𝐽^𝑡𝐽+𝑄)_𝑘𝑉_𝑘)⁻¹𝑉^𝑡_𝑘(𝑅^𝑡𝐽)_𝑘‖, (5.38)and so the full truncated-Newton step is taken. Inductively, this occurs for all subsequence iterates and super linear convergence follows.

The next result shows that under the global assumptions every limit point of the sequence {𝑥_𝑘} satisfies the second-order necessary conditions for a minimum.Lemma 5.7. Let the global assumptions hold, and assume pred_𝑘(Δ)≥𝛼₂(−𝜆₁(𝐻))Δ² for all the symmetric matrix 𝐻∈ℜ^{𝑘𝑗×𝑘𝑗}, ℎ∈ℜ^𝑘𝑗, Δ>0 and some 𝛼₂>0. If {𝑥_𝑘}→𝑥_⋆, then 𝐻(𝑥_⋆) is positive semidefinite, where 𝜆₁(𝐻) is the smallest eigenvalue of the matrix 𝐻.Proof. Suppose to the contrary that 𝜆₁(𝐻(𝑥_⋆))<0. There exists 𝐾 such that if 𝑘≥𝐾, then 𝜆₁(𝐻_𝑘)<(1/2)𝜆₁(𝐻(𝑥_⋆)). For all 𝑘≥𝐾 and for all Δ_𝑘>0, we obtainpred_𝑘(Δ)≥𝛼₂(−𝜆₁(𝐻_𝑘))Δ²_𝑘≥𝛼₂2(−𝜆₁(𝐻(𝑥_⋆)))Δ²_𝑘. (5.39)Using inequality (5.17), we obtain|ared_𝑘(Δ)pred_𝑘(Δ)−1|≤‖𝑦_𝑘(Δ)‖²‖pred_𝑘(Δ)‖¹₀‖𝐻_𝑘(𝜂)−𝐻_𝑘≤‖(1−𝜂)𝑑𝜂2‖𝑦_𝑘(Δ)‖²Δ²_𝑘¹₀‖𝐻_𝑘(𝜂)−𝐻_𝑘‖(1−𝜂)𝑑𝜂𝛼₂(−𝜆₁(𝐻(𝑥_⋆.))) (5.40)Since the last quantity goes to zero as Δ_𝑘→0 and since a Newton step is never taken for 𝑘>𝐾, it follows from Step 3(ii) of Algorithm 3.1 that for 𝑘≥𝐾Δ_𝑘 sufficiently small. Δ_𝑘 cannot be decreased further. Thus, Δ is bounded below, butared_𝑘(Δ)≥𝜉pred_𝑘𝜉(Δ)≥𝛼₂2(−𝜆₁(𝐻(𝑥_⋆)))Δ²_𝑘. (5.41)Since 𝑓 is bounded below, ared_𝑘(Δ)→0, so Δ_𝑘→0 which is a contradiction. Hence 𝜆₁(𝐻(𝑥_⋆))≥0. This completes the proof of the lemma.

6. Concluding Remarks

In this paper, we have shown that the implicitly restarted Arnoldi method can be combined with the Newton iteration and trust region strategy to obtain a globally convergent algorithm for solving large-scale NLS problems.

The main restriction of this scheme is that the subspace onto which a given Newton step is projected must solve the Newton equations with a certain accuracy which is monitored by the termination condition. This is enough to guarantee convergence.

This theory is sufficiently general that is hold for any algorithm that projects the problem on a lower dimensional subspace. The convergence results indicate promise for this research to solve large-scale NLS problems. Our next step is to investigate the performance of this algorithm on some NLS problems. The results will be reported in our forthcoming paper.

Acknowledgment

The authors would like to thank the editor-in-chief and the anonymous referees for their comments and useful suggestions in the improvement of the final form of this paper.