# Newton-Krylov Type Algorithm for Solving Nonlinear Least Squares Problems

Authors: Mohammedi R. Abdel-Aziz; Mahmoud M. El-Alem
Journal: International Journal of Mathematics and Mathematical Sciences
 (2009)

---

## 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.


---

## Body

## 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 x⋆∈ℜn that minimizes
the function(1.1)f(x)=∥R(x)∥22=∑i=1m(Ri(x))2,where R(x)=(R1(x),R2(x),…,Rm(x))t,Ri:ℜn→ℜ,1≤i≤m 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 findings(Δ) so that(1.2)ϕ(s(Δ))=min{ϕ(s):∥s∥≤Δ},where Δ is some
positive constant, ∥·∥ is the
Euclidean norm in ℜn, and(1.3)ϕ(s)=gts+12stHs,where g∈ℜn and H∈ℜn×n 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 s(Δ) 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 H 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 Section2, 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 xk and a Krylov
subspace Sk. Denote the dimension of Sk to be kj and {v1(k),v2(k),…,vkj(k)} be a set of
linearly independent vectors in Sk. The next iterate xk+1 is such that
the increment xk+1−xk∈Sk. Thus, we have ∥Rj(xk+Vky)∥=0,forj=1,2,…,m and y∈ℜkj, where the matrix Vk=[v1(k),v2(k),…,vkj(k)]. The second-order Taylor expansion of ∥R(x)∥22 at xk is(2.1)∥R(xk)+Jks∥22+stQks,where Qk∈ℜn×n is defined by Qk=∑j=1nRj(xk)∇2Rj(xk) and the
jacobian matrix is Jk=J(xk)=[∇R1(xk),∇R2(xk),…,∇Rm(xk)]t. If we consider s∈Sk, we get the quadratic model(2.2)qk(y)=∥R(xk)+JkVky∥22+ytAky,where Ak∈ℜkj,kj approximates
the reduced matrix VktQkVk=∑j=1mRj(xk)Vkt∇2Rj(xk)Vk.The first order necessary conditions of (2.1)
is(2.3)(JktJk+Qk)s+RktJk=0.Thus, a solution of the
quadratic model is a solution to an equation of the form(2.4)(JktJk+Qk)s=−RktJk.The model trust region algorithm generates a sequence
of pointsxk, and at the kth stage of the
iteration the quadratic model of problem (1.2) has the form(2.5)ϕk(s)=JktRks+12st(JktJk+Qk)s.At this stage, an initial value
for the trust region radius Δk is also
available. An inner iteration is then performed which consists of using the
current trust region radius, Δk, and the information contained in the quadratic model
to compute a step, s(Δk). Then a comparison of the actual reduction of the
objective function(2.6)aredk(Δk)=∥R(xk)∥22−∥R(xk+sk(Δk))∥22,and the reduction predicted by
the quadratic model(2.7)predk(Δk)=∥R(xk)∥22−qk(y)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 s(Δ), and how the trust region radius Δ is changed. We
also leave unspecified the selection of (JktJk+Qk) except to
restrict it to be symmetric. Details on these points will be addressed in our
forth coming paper.The solutions of the
quadratic model (2.2) is a solution to an equation of the form(2.8)[(JtJ+Q)+μI]s=−RtJ,with μ≥0, μ(sts−Δ2)=0 and [(JtJ+Q)+μI] is positive
semidefinite. This represents the first-order necessary conditions concerning
the pair μ and s, where s is a solution
to model (2.2) and μ is the Lagrange
multiplier associated with the constraint sts≤Δ2. The sufficient conditions that will ensure s 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 s∈ℜn satisfy [(JtJ+Q)+μI]s=−RtJ with [(JtJ+Q)+μI] positive
semidefinite, if(2.9)μ=0,∥s∥≤Δthenssolves(2.2),∥s∥=Δ,thenssolvesmin{ϕ(s):∥s∥=Δ},μ≥0with∥s∥=Δ,thenssolves(2.2).
Of course, if [(JtJ+Q)+μI] is positive
definite then s is a unique
solution.The main result which is used to prove that the
sequence of gradients tends to zero for modified Newton methods is(2.10)f−ϕ(s)≥12∥RtJ∥min{Δ,∥RtJ∥∥JtJ+Q∥},where s is a solution
to (2.2). The geometric interpretation of this inequality is that for a
quadratic function, any solution s to (2.2)
produces a decrease f−ϕ(s) that is at
least as much as the decrease along the steepest descent direction −RtJ. A proof of this result may be found in [17].


## 3. Algorithmic Framework

The algorithm
we will discuss here requires thatf is twice
differentiable at any point x in the domain
of f. 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 replacef(x) by f^(y)=f(s+Vy), where V∈ℜn×k is an
orthonormal matrix, x∈ℜn and y∈ℜk, then applying the iteration to f^ with initial
guess y0 satisfying x0=s+Vy0 should produce
a sequence {yl} related to the
sequence {xl} by xl=Vyl+s, where xl is produced by
applying the algorithm to f with initial
guess x0;(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,(JtJ+Q), is indefinite and should force the direction to be
equal to the Newton direction when (JtJ+Q) 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).
Givenx0∈ℜn, compute (RtJ)(x0), (JtJ+Q)(x0).
Chooseξ1<ξ2, η1∈(0,1),η2>1, Δ1>0,k=1.(2)
Step 1 (construction a basis for the Krylov
subspace).(a)
Chooseϵ∈(0,1), an initial guess s0 which can be
chosen to be zero, form r0=−RtJ−(JtJ+Q)s0, where RtJ=(RtJ)(xn),(JtJ+Q)=(JtJ+Q)(xn), and compute τ=∥r0∥ and v1=r0/τ.(b)
Forj=1,2,… until
convergence, do.(i)
Form(JtJ+Q)vj and
orthogonalize it against the previous v1,v2,…,vj then
h(i,j)=((JtJ+Q)vj,vi),i=1,2,…,j,
v^j+1=(JtJ+Q)vj−∑i=1jh(i,j)vi,
h(i+1,j)=∥v^j+1∥.
vj+1=v^j+1/h(i+1,j).(ii)
Compute the residual normβj=∥RtJ+(JtJ+Q)sj∥ of the solution sj that would be
obtained if we stopped.(iii)
Ifβj≤ϵj∥RtJ∥, where ϵj≤(1/𝒟j)∈(0,ϵ], and 𝒟j is an
approximation to the condition number of (JtJ+Q)j, then set k=j and go to 3, else
continue.(3)
Step 2 (compute the trial step).(i)
Construct the local model,(3.1)ψk(y)=fk+hkty+12ytHky,where hk=Vkt(RtJ)k∈ℜkj and Hk=Vkt(JtJ+Q)kVk∈ℜkj×kj.(ii)
Compute the solutiony to the problem min{ψk(y):∥y∥≤Δk} and set sk=Vky.(4)
Step 3 (test the step and updateΔk).
(i)
Evaluatefk+1=f(xk+Vky), aredk=fk−f(xk+Vky) and predk=fk−ψ(Vky).(ii)
If(aredk/predk)<ξ1 then
beginΔk=η1Δk go to 3(i)
Comment: Do not accept the step “∥y∥”, reduce the
trust region radius “Δk” and compute
a new step
Else ifξ1≤(aredk/predk)≤ξ2 then
xk+1=xk+Vky
Comment: Accept the step and keep the previous trust
region radius the same
Else(aredk/predk)>ξ2 then
Δk=η2Δk
Comment: Accept the step and increase the trust region
radius
end if(5)
Step 4.
Setk+1←k 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
Algorithm3.1 be {xk}; for such sequence we assume(G1)
for allk,xk and xk+sk∈Ω where Ω⊂ℜn is a convex set;(G2)
f∈C2(Ω) and (RtJ)=∇f(x)≠0 for all x∈Ω;(G3)
f(x) is bounded in
norm in Ω and (JtJ+Q) is nonsingular
for all x∈Ω;(G4)
for allsk such that xk+sk∈Ω, the termination condition, ∥(RtJ)k+(JTJ+Q)sk∥≤ϵj∥(RtJ)k∥, ϵk∈(0,1) is satisfied;An immediate consequence of the global assumptions is
that the Hessian matrix(JtJ+Q)(x) is bounded,
that is, there exists a constant β1>0 such that ∥(JtJ+Q)(x)∥<β1. Therefore, the condition numbers of the Hessian, 𝒟k=cond2((JtJ+Q)k) are bounded.Assumption G3 that(JtJ+Q) 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 v^k+1=0 so that vk+1 cannot be
formed, or if Hkmax 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 v^k+1=0 implies (JtJ+Q)k is nonsingular
and xk 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 (JtJ+Q)k is nonsingular
for some k among 1,2,…,kmax−1. If such k exists, we can
compute xj and then
restart the algorithm using xk as the new
initial guess x0. It may not be possible to do this.To handle the problem of singularHk, we can use a technique similar to that done in the
full dimensional space. If Hk is singular, or
its condition number is greater than 1/ν, where ν is the machine
epsilon, then we perturb the quadratic model ψ(y) to(4.1)ψ^(y)=ψ(y)+12γyty=f+hty+12yt(Hk+γI)y,where γ=nν∥Hk∥1. The condition number of (Hk+γI) 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,𝒮k. Indeed in Algorithm 3.1, the Krylov subspace will
contain the steepest descent direction because v1=−RtJ/∥RtJ∥ and the Hessian
is nonsingular;(2)
there is no
difficulty, if one required the Hessian matrix,(JtJ+Q)(x) to be positive
definite for all x∈Ω, 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 (JtJ+Q)(x) 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,kj can be chosen
as large as the termination condition will eventually be satisfied, but, when kj is too large
the computational cost and the storage become too high. An alternative to that
is to restart the algorithm keeping (JtJ+Q) nonsingular and kj<n. 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 Algorithm3.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 anys∈Ω, one has(5.1)|(RtJ)kts|∥s∥2≥𝒟k−ϵk∥(JtJ+Q)k∥(1+ϵk)𝒟k∥(RtJ)k∥>0.Proof.
Supposerk be the residual
associated with s so that rk=(RtJ)k+(JtJ+Q)ks and (RtJ)k≠0. Then ∥rk∥≤ϵk∥(RtJ)k∥, and s=−(JtJ+Q)k−1((RtJ)k−rk). So,(5.2)(RtJ)kts=(RtJ)kt[−(JtJ+Q)k−1((RtJ)k−rk)]=−(RtJ)kt(JtJ+Q)k−1(RtJ)k+(RtJ)ktrk.Hence,(5.3)|(RtJ)kts|∥s∥=|(RtJ)kt(JtJ+Q)k−1(RtJ)k−(RtJ)ktrk|∥(JtJ+Q)k−1((RtJ)k−rk)∥≥(RtJ)kt(JtJ+Q)k−1(RtJ)k−|(RtJ)ktrk|∥(JtJ+Q)k−1((RtJ)k−rk)∥.Next, ∥rk∥≤ϵ∥(RtJ)k∥ implies |(RtJ)ktrkj|≤ϵk∥(RtJ)k∥2 which gives(5.4)(RtJ)kt(JtJ+Q)k−1(RtJ)k−|(RtJ)ktrk|≥(∥(JtJ+Q)k−1∥−ϵk)∥(RtJ)k∥2,∥(JtJ+Q)k−1((RtJ)k−rk)∥≤∥(JtJ+Q)k−1∥∥(RtJ)k∥+∥(JtJ+Q)k−1rk∥≤(1+ϵk)∥(JtJ+Q)k−1∥∥(RtJ)k∥.Introduce (5.4)
into
(5.3), we obtain the following:(5.5)|(RtJ)kts|∥s∥≥(∥(JtJ+Q)k−1∥−ϵk)∥(RtJ)k∥2(1+ϵk)∥(JtJ+Q)k−1∥∥(RtJ)k∥≥(∥(JtJ+Q)k−1∥−ϵk)(1+ϵk)∥(JtJ+Q)k−1∥∥(RtJ)k∥=(𝒟k−ϵk∥(JtJ+Q)k∥(1+ϵk)𝒟k)∥(RtJ)k∥.Inequality (5.1) can be written
as(5.6)|(RtJ)kts|∥s∥∥(RtJ)k∥≥𝒟k−ϵk∥(JtJ+Q)k∥(1+ϵk)𝒟k.This condition (5.6) tells us at
each iteration of Algorithm 3.1, we require that the termination condition
holds. Since the condition numbers, 𝒟k, are bounded from above, condition (5.6) gives(5.7)cosθ≥𝒟−ϵ∥(JtJ+Q)k∥(1+ϵ)𝒟,where θ is the angle
between (RtJ)k and s and 𝒟 are the upper
bound of the condition numbers 𝒟k. 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(5.8)∥Vkt(RtJ)k∥≥𝒟k−ϵk∥(JtJ+Q)k∥(1+ϵk)𝒟k∥(RtJ)k∥.Proof.
Supposes=Vky be a vector of 𝒮k satisfying the
termination condition, ∥(RtJ)k+(JtJ+Q)ks∥≤ϵk∥(RtJ)k∥,ϵk∈(0,1). From Lemma 5.1 and the fact that ∥s∥=∥Vy∥=∥y∥, we obtain(5.9)|(RtJ)ktVky|∥y∥≥𝒟k−ϵk∥(JtJ+Q)k∥(1+ϵk)𝒟k∥(RtJ)k∥.From Cauchy-Schwartz
inequality, we obtain(5.10)|(RtJ)ktVky|≤∥Vkt(RtJ)k∥∥y∥.Combining formulas (5.9) and
(5.10), we obtain the result.The following Lemma establishes that Algorithm3.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 lety be a solution
of the minimization problem, min∥y∥≤Δψ(y). Then(5.11)∥hk∥min{Δk,∥hk∥∥Hk∥}≥αk∥(RtJ)k∥min{Δk,αk∥(RtJ)k∥∥(JtJ+Q)k∥},where αk=(𝒟k−ϵk∥(JtJ+Q)k∥)/(1+ϵk)𝒟k.Proof.
Since,ψ(y)=ϕ(Vky)=fk+(RtJ)ktVky+(1/2)(Vky)t(JtJ+Q)(Vky). Thus ∇ψ(0)=Vkt(RtJ)k. From Step 3 of Algorithm 3.1, we obtain(5.12)fk−fk+1≥ξ1(fk−ϕk(y))(5.13)≥ξ12∥Vkt(RtJ)k∥min{Δk,∥Vkt(RtJ)k∥∥Vkt(JtJ+Q)kVk∥}(5.14)=ξ12∥hk∥min{Δk,∥hk∥∥Hk∥}.We have to convert the lower
bound of this inequality to one involving ∥(RtJ)k∥ rather than ∥Vkt(RtJ)k∥. Using Lemma 5.2, we obtain ∥hk∥≥αk∥(RtJ)k∥, but, ∥Hk∥=∥Vkt(JtJ+Q)kVk∥≤∥(JtJ+Q)k∥. Substituting in (5.14), we obtain(5.15)∥hk∥2min{Δk,∥hk∥∥Hk∥}≥αk2∥(RtJ)k∥min{Δk,αk∥(RtJ)k∥∥(JtJ+Q)k∥},which completes the proof.The following two facts will be used in the remainder
of the proof.Fact 1.
By Taylor's
theorem for anyk and Δ>0, we have(5.16)|aredk(Δ)−predk(Δ)|=|(fk−fk+1)−(−hktyk(Δ)−12ykt(Δ)Vkt(JtJ+Q)kVkyk(Δ))|=|12ykt(Δ)Hkyk(Δ)−∫01ykt(Δ)Hk(η)yk(Δ)(1−η)dη|≤∥yk(Δ)∥2∫01∥Hk−Hk(η)∥(1−η)dη.where, Hk(η)=Vkt(JtJ+Q)(xk+ηsk)Vk with sk=xk+1−xk=Vyk. Thus,(5.17)|aredk(Δ)predk(Δ)−1|≤∥yk(Δ)∥2|predk(Δ)|∫01∥Hk−Hk(η)∥(1−η)dη.Fact 2.
For any
sequence{xk} generated by an
algorithm satisfying the global assumptions, the related sequence {fk} is
monotonically decreasing and bounded from below, that is, fk→f⋆ as k→∞.The next result establishes that every limit point of
the sequence{xk} satisfies the
first-order necessary conditions for a minimum.Lemma 5.4.
Let the global assumptions hold and Algorithm3.1 is applied to f(x), generating a sequence {xk},xk∈Ω,andk=1,2,…, then hk→0.Proof.
Sinceϵk≤ϵ<1 for all k, we have(5.18)αk≥𝒟−ϵ∥(JtJ+Q)k∥(1+ϵ)𝒟=α>0∀k.Consider any xl with ∥(RtJ)l∥≠0, ∥(RtJ)(x)−(RtJ)l∥≤β1∥x−xl∥. Thus, if ∥x−xk∥≤(∥(RtJ)k∥/2β1), then ∥(RtJ)(x)∥≥∥(RtJ)k(RtJ)l∥−∥(RtJ)(x)−(RtJ)k∥≥(∥(RtJ)k∥/2). Let, ρ=∥(RtJ)k∥/2 and ℬρ={x:∥x−xk∥<ρ}. There are two cases, either for all k≥l,xl∈ℬρ or eventually {xl} leaves the ball ℬρ. Suppose xl∈ℬρ for all k≥l, then ∥gl∥≥∥gk∥/2=σ¯ for all l>k. Thus by Lemma 5.3, we have(5.19)predk(Δ)≥ξ1ασ¯2min{Δk,ασ¯∥Hk∥}∀k≥l.Put ασ¯=σ, introduce inequality (5.19) into (5.17) and since ∥Hk∥≤∥(JtJ+Q)k∥≤β1, we
obtain(5.20)|aredk(Δ)predk(Δ)−1|≤∥yk∥2∫01∥Hk−Hk(η)∥(1−η)dηξ1σmin{Δk,σ/β1}≤Δk2(2β1)ξ1σmin{Δk,σ/β1}≤2Δkβ1ξ1σ,Δk≤σβ1∀k≥l.This gives for Δk sufficiently
small and k≥l that(5.21)aredk(Δ)predk(Δ)>ξ2.In addition, we have(5.22)∥(JtJ+Q)k−1hk∥≥∥hk∥∥Hk∥≥α∥(RtJ)k∥β1≥σβ1.Inequalities (5.21) and (5.22)
tell us that, for Δk sufficiently
small, we cannot get a decrease in Δk in Step 3(ii)
of Algorithm 3.1. It follows that Δk is bounded away
from 0, but, since(5.23)fk−fk+1=aredk(Δ)≥ξ1predk(Δ)≥ξ1σmin{Δk,σβ1},and since f is bounded from
below, we must have Δk→0 which is a
contradiction. Therefore, {xk} must be outside ℬρ for some k>l.
Letxj+1 be the first
term after xl that does not
belong to ℬρ. From inequality (5.23), we get(5.24)f(xl)−f(xj+1)=∑k=lj(fk−fk+1)≥∑j+1lξ1predk(Δ)≥∑k=ljξ1σmin{Δk,σβ1}.If Δk≤σ/β1 for all l≤k≤j, we have(5.25)f(xl)−f(xj+1)≥ξ1σ∑k=ljΔk≥ξ1σρ.If there exists at least one k with Δk>σ/β1, we get(5.26)f(xm)−f(xm+1)≥ξ1σ2β1.In either case, we
have(5.27)f(xm)−f(xj+1)≥ξ1σmin{ρ,σβ1}=ξ1α∥(RtJ)k∥2min{α∥(RtJ)k∥2β1,α∥(RtJ)k∥2β1}≥ξα24β1∥(RtJ)k∥2.Since, f is bounded
below, {fk} is a
monotonically decreasing sequence (i.e.,fk→f⋆). Hence,(5.28)∥(RtJ)k∥2≤4β1ξ1α2(fk−f⋆),(i.e.,(RtJ)k→0 as k→∞). 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 Algorithm3.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,Hk=Vkt(JtJ+Q)(xk)Vk for all k, (JtJ+Q)(x) is Lipschitz
continuous with L, and x* is a limit
point of the sequence {xk} with (JtJ+Q)(x*) be positive
definite. Then {xk} converges q-Superlinearly
to x⋆.Proof.
Let{xkm} be a
subsequence of the sequence {xk} that converges
to x⋆, where, x⋆ is a limit of {xk}. From Lemma 5.4, (RtJ)(x⋆)=0. Since (JtJ+Q)(x⋆) is positive
definite and (JtJ+Q)(x) is continuous
there exists δ1>0 such that if ∥x−x⋆∥<δ1, then (JtJ+Q)(x) is positive
definite, and if x≠x⋆ then (RtJ)(x⋆)≠0. Let, B1={x:∥x−x⋆∥<δ1}.
Since(RtJ)(x⋆)=0, we can find δ2>0 with δ2<δ1/4 and ∥(JtJ+Q)−1∥∥(RtJ)∥≤δ1/2 for all x∈B2={x:∥x−x⋆∥<δ2}.
Findlo such that f(xlo)<inf{f(x):x∈B1−B2} and xlo∈B2. Consider any xi with i≥kl,xi∈B2.
We claim thatxi+1∈B2, which implies that the entire sequence beyond xklo∈B2. Suppose xi+1 does not belong
to B2. Since fi+1<fklo, xi+1 does not belong
to B1 either.
So(5.29)Δi≥∥xi+1−xi∥≥∥xi+1−x⋆∥−∥xi−x⋆∥≥δ1−δ14=3δ14>δ12≥∥(JtJ+Q)i−1∥∥(RtJ)i∥≥∥Hi−1∥∥hi∥.
So,
the truncated-Newton step is within the trust region, we obtain
(5.30)yi(Δi)=−(Vit(JtJ+Q)iVi)−1Vit(RtJ)i.Since ∥yi(Δ)∥<δ1/2, it follows that xi+1∈B1, which is a contradiction. Hence, for all k≥klo,xk∈B2, and so since f(xk) is strictly
decreasing sequence and x⋆ is the unique
minimizer of f in B2, we have that {xk}→x⋆.The following lemma establishes the rate of
convergence of the sequence of iterates generated by Algorithm3.1 to be q-super
linear.Lemma 5.6.
Let the assumptions of Lemma5.5
hold. Then the rate of convergence of the sequence {xk} is q-Super
linear.Proof.
To show that the convergence rate is super linear, we
will show eventually that∥Hk−1hk∥ will always be
less than Δj and hence the
truncated-Newton step will always be taken. Since (JtJ+Q)(x⋆) is positive
definite, it follows that the convergence rate of {xk} to x⋆ 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
onpredk(Δ). By assumptions, for all k large enough, Hk is positive
definite. Hence, either the truncated-Newton step is longer than the trust
region radius, or yk(Δ) is the
truncated-Newton step. In either case(5.31)∥yk(Δ)∥≤∥(Vkt(JtJ+Q)jVj)−1Vkt(RtJ)k∥≤∥(Vkt(JtJ+Q)kVk)−1∥∥Vkt(RtJ)k∥,and so it follows that ∥Vkt(RtJ)k∥≥∥sk(Δ)∥/∥(Vkt(JtJ+Q)kVk)−1∥. By Lemma 5.3, for all k large enough,
we have(5.32)predk(Δ)≥ξ1α2∥yk(Δ)∥∥(JtJ+Q)k−1∥min{Δk,α∥(RtJ)k∥∥(JtJ+Q)k∥}.Thus,(5.33)predk(Δ)≥ξ1α2∥yk(Δ)∥∥(JtJ+Q)k−1∥min{∥yk(Δ)∥,α∥yk(Δ)∥∥(JtJ+Q)k∥∥(JtJ+Q)k−1∥}≥ξ1α22∥yk(Δ)∥2𝒟∥(JtJ+Q)k−1∥.So, by the continuity of (JtJ+Q)(x), for all k large enough,
we obtain(5.34)predk(Δ)≥ξ1α24∥yk(Δ)∥2∥(JtJ+Q)−1(x⋆)∥.Finally, by the argument leading
up to Fact 1 and the Lipschitz continuity assumption, we have(5.35)|aredk(Δ)−predk(Δ)|≤∥yk(Δ)∥3L2.Thus, for any Δk>0 and k large enough,
we obtain(5.36)|aredk(Δ)predk(Δ)−1|≤∥sk(Δ)∥3L2.4∥(JtJ+Q)−1(x⋆)∥ξ1α2∥yk(Δ)∥2=2Lξ1α2∥(JtJ+Q)−1(x⋆)∥∥yk(Δ)∥≤2L∥(JtJ+Q)−1(x⋆)∥ξ1α2Δk.Thus by Step 3(ii) of Algorithm
3.1, there is a Δ¯ such that if Δk−1<Δ¯, then Δk will be less
than Δk−1 only if Δk≥∥(Vk−1t(JtJ+Q)k−1Vk−1)−1Vk−1t(RtJ)k−1∥. It follows from the superlinear convergence of the
truncated-Newton method that for xk−1 close enough to x⋆ and k large
enough,(5.37)∥(Vkt(JtJ+Q)kVk)−1Vkt(RtJ)k∥<∥(Vk−1t(JtJ+Q)k−1Vk−1)−1Vk−1t(RtJ)k−1∥.Now, if Δk is bounded away
from 0 for all large k, then we are done. Otherwise, if for an arbitrarily
large k, Δk is reduced,
that is, Δk<Δk−1, then we have(5.38)Δk≥∥(Vk−1t(JtJ+Q)k−1Vk−1)−1Vk−1t(RtJ)k−1∥>∥(Vkt(JtJ+Q)kVk)−1Vkt(RtJ)k∥,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{xk} satisfies the
second-order necessary conditions for a minimum.Lemma 5.7.
Let the global assumptions hold, and
assumepredk(Δ)≥α¯2(−λ1(H))Δ2 for all the
symmetric matrix H∈ℜkj×kj, h∈ℜkj, Δ>0 and some α¯2>0. If {xk}→x⋆, then H(x⋆) is positive
semidefinite, where λ1(H) is the smallest
eigenvalue of the matrix H.Proof.
Suppose to the contrary thatλ1(H(x⋆))<0. There exists K such that if k≥K, then λ1(Hk)<(1/2)λ1(H(x⋆)). For all k≥K and for all Δk>0, we obtain(5.39)predk(Δ)≥α¯2(−λ1(Hk))Δk2≥α¯22(−λ1(H(x⋆)))Δk2.Using inequality (5.17), we
obtain(5.40)|aredk(Δ)predk(Δ)−1|≤∥yk(Δ)∥2∥predk(Δ)∥∫01∥Hk(η)−Hk∥(1−η)dη≤2∥yk(Δ)∥2Δk2∫01∥Hk(η)−Hk∥(1−η)dηα¯2(−λ1(H(x⋆))).Since the last quantity goes to
zero as Δk→0 and since a
Newton step is never taken for k>K, it follows from Step 3(ii) of Algorithm 3.1 that
for k≥KΔk sufficiently
small. Δk cannot be
decreased further. Thus, Δ is bounded
below, but(5.41)aredk(Δ)≥ξpredk(Δ)≥ξα¯22(−λ1(H(x⋆)))Δk2.Since f is bounded
below, aredk(Δ)→0, so Δk→0 which is a
contradiction. Hence λ1(H(x⋆))≥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.



---

*Source: 435851-2009-07-22.xml*