# A Globally Convergent Matrix-Free Method for Constrained Equations and Its Linear Convergence Rate

Authors: Min Sun; Jing Liu
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A matrix-free method for constrained equations is proposed, which is a combination of the well-known PRP (Polak-Ribière-Polyak) conjugate gradient method and the famous hyperplane projection method. The new method is not only derivative-free, but also completely matrix-free, and consequently, it can be applied to solve large-scale constrained equations. We obtain global convergence of the new method without any differentiability requirement on the constrained equations. Compared with the existing gradient methods for solving such problem, the new method possesses linear convergence rate under standard conditions, and a relax factorγ is attached in the update step to accelerate convergence. Preliminary numerical results show that it is promising in practice.


---

## Body

## 1. Introduction

LetF
:
R
n
→
R
n be a continuous nonlinear mapping and C a nonempty closed convex set of R
n. In this paper, we consider the problem of finding x
∈
C such that
(1)
F
(
x
)
=
0
.
Nonlinear constrained equations (1), denoted by CES (F
,
C), arise in various applications, for instance, ballistic trajectory computation and vibration systems [1], the power flow equations [2], chemical equilibrium systems [3], and so forth.In recent years, many numerical methods have been proposed to find a solution of nonsmooth CES (F
,
C), which include the trust region methods [4, 5], the Levenberg-Marquardt method [6], and the projection methods [7–9]. Compared with the trust region method and the Levenberg-Marquardt method, the projection method is more efficient for solving large-scale CES (F
,
C). Noting this, Wang et al. [7] proposed a projection method for solving CES (F
,
C), which possesses global convergence property without the differentiability. A drawback of this method is that it needs to solve a linear equation inexactly at each iteration, and its variants [8, 10] also have this drawback.It is well-known that the spectral gradient method and the conjugate gradient method are two efficient methods for solving large-scale unconstrained optimization problems due to their simplicity and low storage. Recently, La Cruz and Raydan [11] successfully applied the famous spectral gradient method to solve unconstrained equations by using some merit function. Then, Zhang and Zhou [12] presented a spectral gradient projection method (SGP) for solving unconstrained monotone equations, which does not utilize any merit function. Later, the SGP was extended by Yu et al. [9] to solve monotone constrained equations. However, the study of conjugate gradient methods for large-scale (un)constrained equations is relatively rare. Cheng [13] proposed a PRP type method (PRPT) for systems of monotone equations, which is a combination of the well-known PRP method and the hyperplane projection method, and the numerical results in [13] show that the PRPT method performs better than the SGP method in [12].Different from the methods in [7, 8, 10], the methods in [9, 11–13] do not need to solve a linearized equation at each iteration; however, the latter do not investigate the convergent rate, and even we do not know whether they possess the linear convergence rate. In this paper, motivated by the projection methods in [7, 8, 10] and the gradient methods in [9, 12, 13], we propose a matrix-free method for solving nonlinear constrained equations, which can be viewed as a combination of the well-known PRP conjugate gradient method and the famous hyperplane projection method, and it possesses linear convergence rate under standard conditions. The remainder of this paper is organized as follows. Section 2 describes the new method and presents its global convergence analysis. The linear convergence rate of the new method is established in Section 3. Numerical results are reported in Section 4. Finally, some final remarks are included in Section 5.


## 2. Algorithm and Convergence Analysis

LetC
* denote the solution set of CES (F
,
C). Throughout this paper, we assume that C
* is nonempty and F
(
·
) is monotone; that is,
(2)
〈
F
(
x
)
-
F
(
y
)
,
x
-
y
〉
≥
0
,
∀
x
,
y
∈
R
n
,
which implies that the solution set C
* is closed. Then let P
C
(
x
) denote the orthogonal projection of a point x
∈
R
n onto the convex set C, which has the following nonexpansive property:
(3)
∥
P
C
(
x
)
-
P
C
(
y
)
∥
≤
∥
x
-
y
∥
,
∀
x
,
y
∈
R
n
.
Now, we describe the matrix-free method for nonlinear constrained equations.Algorithm 1.
Consider the following.
Step  0. Given an arbitrary initial point x
0
∈
C, the parameters 0
<
ρ
<
1, 0
<
σ
<
r
≤
1, 0
<
γ
<
2, and 0
<
β
min
⁡
<
β
max
⁡. Given the initial steplength β
0
=
1 and set k
:
=
0.
Step  1. If F
(
x
k
)
=
0, then stop; otherwise go to Step  2.
Step  2. Compute d
k by
(4)
d
k
=
{
-
F
(
x
k
)
,
if
k
=
0
,
-
F
(
x
k
)
+
β
k
d
k
-
1
-
θ
k
(
F
(
x
k
)
-
F
(
x
k
-
1
)
)
,
if
k
≥
1
,
where
(5)
β
k
PRP
=
〈
F
(
x
k
)
,
F
(
x
k
)
-
F
(
x
k
-
1
)
〉
∥
F
(
x
k
-
1
)
∥
2
,
111
1
11
1
1
θ
k
=
F
(
x
k
)
⊤
d
k
-
1
∥
F
(
x
k
-
1
)
∥
2
,
1111
1
111
1
1111111
∀
k
≥
1
.
If ∥
d
k
∥
>
∥
F
(
x
k
)
∥
/
r, set d
k
=
-
F
(
x
k
).
Step  3. Find the trial point y
k
=
x
k
+
α
k
d
k, where α
k
=
β
k
ρ
m
k with m
k being the smallest nonnegative integer m such that
(6)
-
〈
F
(
y
k
)
,
d
k
〉
≥
σ
∥
d
k
∥
2
.
Step  4. Compute
(7)
x
k
+
1
=
P
C
[
x
k
-
γ
ξ
k
F
(
y
k
)
]
,
where
(8)
ξ
k
=
〈
F
(
y
k
)
,
x
k
-
y
k
〉
∥
F
(
y
k
)
∥
2
.
Choose an initial steplength β
k
+
1 such that β
k
+
1
∈
[
β
min
⁡
,
β
max
⁡
]. Set k
:
=
k
+
1 and go to Step  1.Remark 2.
Obviouslyd
k, defined by (4), is motivated by [14], and it is not difficult to deduce that d
k satisfies
(9)
F
(
x
k
)
⊤
d
k
=
-
∥
F
(
x
k
)
∥
2
.
Therefore, by Cauchy-Schwartz inequality, we have ∥
d
k
∥
≥
∥
F
(
x
k
)
∥. This together with Step  2 of Algorithm 1 implies
(10)
∥
F
(
x
k
)
∥
≤
∥
d
k
∥
≤
∥
F
(
x
k
)
∥
r
.Remark 3.
In (7), we attach a relax factor γ
∈
(
0,2
) (better when close to 2) to F
(
y
k
) based on numerical experiences.Remark 4.
Line search (6) is different from that of [12, 13], which is well-defined by the following Lemma.Lemma 5.
For allk
≥
0, there exists a nonnegative number m
k satisfying (6).Proof.
In fact, ifd
k
=
0, then from (10), we have ∥
F
(
x
k
)
∥
=
0, which means that Algorithm 1 terminates with x
k being a solution of CES (F
,
C). Now, we consider d
k
≠
0 for all k. For the sake of contradiction, we suppose that there exists k
0
≥
0 such that (6) is not satisfied for any nonnegative integer m; that is,
(11)
-
〈
F
(
x
k
0
+
β
k
0
ρ
m
d
k
0
)
,
d
k
0
〉
<
σ
∥
d
k
0
∥
2
,
∀
m
≥
1
.
Letting m
→
∞ and using the continuity of F
(
·
) yield
(12)
-
〈
F
(
x
k
0
)
,
d
k
0
〉
≤
σ
∥
d
k
0
∥
2
.
On the other hand, by (10), we obtain
(13)
-
〈
F
(
x
k
0
)
,
d
k
0
〉
=
∥
F
(
x
k
0
)
∥
2
≥
r
∥
d
k
0
∥
2
,
which together with (12) means that σ
≥
r; however, this contradicts the fact that r
>
σ
>
0. Therefore the assertion holds. This completes the proof.Lemma 6.
Suppose thatF
(
·
) is monotone and let {
x
k
} and {
y
k
} be the sequences generated by Algorithm 1; then {
x
k
} and {
y
k
} are both bounded; furthermore, it holds that
(14)
lim
⁡
k
→
∞
α
k
∥
d
k
∥
2
=
0
.Proof.
From (6), we have
(15)
〈
F
(
y
k
)
,
x
k
-
y
k
〉
≥
σ
α
k
∥
d
k
∥
2
>
0
.
For any x
*
∈
C
*, from (3), the nonexpansiveness of the projection operator, it holds that
(16)
∥
x
k
+
1
-
x
*
∥
2
=
∥
P
C
[
x
k
-
γ
ξ
k
F
(
y
k
)
]
-
x
*
∥
2
≤
∥
x
k
-
γ
ξ
k
F
(
y
k
)
-
x
*
∥
2
=
∥
x
k
-
x
*
∥
2
-
2
γ
ξ
k
〈
F
(
y
k
)
,
x
k
-
x
*
〉
+
γ
2
ξ
k
2
∥
F
(
y
k
)
∥
2
.
By the monotonicity of mapping F
(
·
), we have
(17)
〈
F
(
y
k
)
,
x
k
-
x
*
〉
=
〈
F
(
y
k
)
,
x
k
-
y
k
〉
+
〈
F
(
y
k
)
,
y
k
-
x
*
〉
≥
〈
F
(
y
k
)
,
x
k
-
y
k
〉
+
〈
F
(
x
*
)
,
y
k
-
x
*
〉
=
〈
F
(
y
k
)
,
x
k
-
y
k
〉
.
Substituting (15) and (17) into (16), we have
(18)
∥
x
k
+
1
-
x
*
∥
2
≤
∥
x
k
-
x
*
∥
2
-
2
γ
ξ
k
〈
F
(
y
k
)
,
x
k
-
y
k
〉
+
γ
2
ξ
k
2
∥
F
(
y
k
)
∥
2
=
∥
x
k
-
x
*
∥
2
-
γ
(
2
-
γ
)
〈
F
(
y
k
)
,
x
k
-
y
k
〉
2
∥
F
(
y
k
)
∥
2
≤
∥
x
k
-
x
*
∥
2
-
γ
(
2
-
γ
)
σ
2
α
k
2
∥
d
k
∥
4
∥
F
(
y
k
)
∥
2
,
which together with γ
∈
(
0,2
) indicates that, for all k,
(19)
∥
x
k
+
1
-
x
*
∥
≤
∥
x
k
-
x
*
∥
,
which shows that the sequence {
x
k
} is bounded. By (10), it holds that {
d
k
} is bounded and so is {
y
k
}. Then, by the continuity of F
(
·
), there exists a constant M
>
0 such that ∥
F
(
y
k
)
∥
≤
M for all k. Therefore it follows from (18) that
(20)
γ
(
2
-
γ
)
σ
2
M
2
∑
k
=
0
∞
‍
α
k
2
∥
d
k
∥
4
≤
∑
k
=
0
∞
‍
(
∥
x
k
-
x
*
∥
2
-
∥
x
k
+
1
-
x
*
∥
2
)
<
∞
,
which implies that the assertion (14) holds. The proof is completed.Now, we prove the global convergence of Algorithm1.Theorem 7.
Suppose that the conditions in Lemma6 hold. Then the sequence {
x
k
} generated by Algorithm 1 globally converges to a solution of CES (F
,
C).Proof.
We consider the following two possible cases.
Case  1. Consider liminf
⁡
k
→
∞
∥
d
k
∥
=
0. Thus, by (10), we have liminf
⁡
k
→
∞
∥
F
(
x
k
)
∥
=
0. This together with the continuity of F
(
·
) implies that the sequence {
x
k
} has some accumulation point x
¯ such that F
(
x
¯
)
=
0. From (19), it holds that {
∥
x
k
-
x
¯
∥
} converges, and since x
¯ is an accumulation point of {
x
k
}, it must hold that {
x
k
} converges to x
¯.
Case  2. Consider liminf
⁡
k
→
∞
∥
d
k
∥
>
0. Then by (14), it follows that lim
⁡
k
→
∞
α
k
=
0. Therefore, from the line search (6), for sufficiently large k, we have
(21)
-
〈
F
(
x
k
+
β
k
ρ
m
k
-
1
d
k
)
,
d
k
〉
<
σ
∥
d
k
∥
2
.
Since {
x
k
}
,
{
d
k
} are both bounded, we can choose a sequence and letting k
→
∞ in (21), we can obtain
(22)
-
〈
F
(
x
¯
)
,
d
¯
〉
≤
σ
∥
d
¯
∥
2
,
where x
¯ and d
¯ are limit points of corresponding subsequences. On the other hand, by (10), we obtain
(23)
-
〈
F
(
x
k
)
,
d
k
〉
=
∥
F
(
x
k
)
∥
2
≥
r
∥
d
k
∥
2
.
Letting k
→
∞ in the above inequality, we obtain
(24)
-
〈
F
(
x
¯
)
,
d
¯
〉
≥
r
∥
d
¯
∥
2
.
Thus, by (22) and (24), we get r
≤
σ, and this contradicts the fact that r
>
σ
>
0. Therefore, liminf
⁡
k
→
∞
∥
d
k
∥
>
0 does not hold. This completes the proof.


## 3. Convergence Rate

By Theorem7, we know that the sequence {
x
k
} generated by Algorithm 1 converges to a solution of CES (F
,
C). In what follows, we always assume that x
k
→
x
* as k
→
∞, where x
*
∈
C
*. To establish the local convergence rate of the sequence generated by Algorithm 1, we need the following assumption.Assumption 8.
Forx
*
∈
C
*, there exist three positive constants δ, c, and L such that
(25)
c
dist
⁡
(
x
,
C
*
)
≤
∥
F
(
x
)
∥
,
∀
x
∈
N
(
x
*
,
δ
)
,
(26)
∥
F
(
x
)
-
F
(
y
)
∥
≤
L
∥
x
-
y
∥
,
∀
x
,
y
∈
N
(
x
*
,
δ
)
,
where dist(
x
,
C
*
) denotes the distance from x to the solution set C
*, and
(27)
N
(
x
*
,
δ
)
=
{
x
∈
R
n
∣
∥
x
-
x
*
∥
≤
δ
}
.
Now, we analyze the convergence rate of the sequence {
x
k
} generated by Algorithm 1 under conditions (25) and (26).Lemma 9.
If the conditions in Assumption8 hold, then the sequence {
α
k
} generated by line search (6) has a positive bound from below.Proof.
We only need to prove that for sufficiently largek, α
k has a positive bound from below. If α
k
≤
β
k, then by the construction of α
k, we have
(28)
-
〈
F
(
x
k
+
β
k
α
k
ρ
-
1
d
k
)
,
d
k
〉
<
σ
∥
d
k
∥
2
.
In addition, by (10), we have
(29)
-
〈
F
(
x
k
)
,
d
k
〉
=
∥
F
(
x
k
)
∥
2
≥
r
∥
d
k
∥
2
.
Then, by the above two inequalities, we can obtain
(30)
〈
F
(
x
k
+
β
k
α
k
ρ
-
1
d
k
)
-
F
(
x
k
)
,
d
k
〉
≥
(
r
-
σ
)
∥
d
k
∥
2
.
On the other hand, from (26), we have
(31)
〈
F
(
x
k
+
β
k
α
k
ρ
-
1
d
k
)
-
F
(
x
k
)
,
d
k
〉
≤
L
β
k
α
k
ρ
∥
d
k
∥
2
.
By (30) and (31), for k sufficiently large we obtain
(32)
α
k
≥
ρ
(
r
-
σ
)
L
β
k
≥
ρ
(
r
-
σ
)
L
β
max
⁡
.
Therefore, there is a positive constant α, such that
(33)
α
k
≥
α
,
for all k. The proof is completed.Theorem 10.
In addition to the assumptions in Theorem7, if conditions (25) and (26) hold, then the sequence {
dist
⁡
(
x
k
,
C
*
)
} Q-linearly converges to 0; hence the whole sequence {
x
k
} converges to x
* R-linearly.Proof.
Letz
k
∈
C
* be the closest solution to x
k. That is, ∥
x
k
-
z
k
∥
=
dist
⁡
(
x
k
,
C
*
). By (18), we have
(34)
∥
x
k
+
1
-
z
k
∥
2
≤
∥
x
k
-
z
k
∥
2
-
γ
(
2
-
γ
)
〈
F
(
y
k
)
,
x
k
-
y
k
〉
2
∥
F
(
y
k
)
∥
2
.
For sufficiently large k, it follows from (10) and (26) that
(35)
∥
F
(
y
k
)
∥
=
∥
F
(
y
k
)
-
F
(
z
k
)
∥
≤
L
∥
y
k
-
z
k
∥
≤
L
(
∥
x
k
-
y
k
∥
+
∥
x
k
-
z
k
∥
)
≤
L
(
β
max
⁡
∥
d
k
∥
+
∥
x
k
-
z
k
∥
)
≤
L
(
β
max
⁡
∥
F
(
x
k
)
∥
r
+
∥
x
k
-
z
k
∥
)
=
L
(
β
max
⁡
∥
F
(
x
k
)
-
F
(
z
k
)
∥
r
+
∥
x
k
-
z
k
∥
)
≤
L
(
β
max
⁡
L
r
+
1
)
∥
x
k
-
z
k
∥
=
L
(
β
max
⁡
L
r
+
1
)
dist
⁡
(
x
k
,
C
*
)
.
Thus, from (6), (10), (25), and (33), for sufficiently large k, it holds that
(36)
〈
F
(
y
k
)
,
x
k
-
y
k
〉
≥
σ
α
k
∥
d
k
∥
2
≥
σ
α
∥
d
k
∥
2
≥
σ
α
∥
F
(
x
k
)
∥
2
≥
σ
α
c
2
dis
t
2
(
x
k
,
C
*
)
.
Substituting the above two inequalities into (34), we have
(37)
dis
t
2
(
x
k
+
1
,
C
*
)
≤
∥
x
k
+
1
-
z
k
∥
2
≤
(
1
-
σ
α
r
2
c
2
γ
(
2
-
γ
)
L
2
(
β
max
⁡
L
+
r
)
2
)
dis
t
2
(
x
k
,
C
*
)
,
which implies that the sequence {
dist
⁡
(
x
k
,
C
*
)
}
Q-linearly converges to 0. Therefore, the whole sequence {
x
k
} converges to x
*
R-linearly. The proof is completed.


## 4. Numerical Results

In this section, we test Algorithm1 and compared it with the projection method in [7] and the spectral gradient projection method in [9]. We give the following two simple problems to test the efficiency of the three methods.Problem 11.
The mappingF
(
·
) is taken as F
(
x
)
=
(
f
1
(
x
)
,
f
2
(
x
)
,
…
,
f
n
(
x
)
)
⊤, where
(38)
f
i
(
x
)
=
e
x
i
-
1
,
for
i
=
1,2
,
…
,
n
and C
=
R
+
n. Obviously, this problem has a unique solution x
*
=
(
0,0
,
…
,
0
)
⊤.Problem 12.
The mappingF
(
·
) is taken as F
(
x
)
=
(
f
1
(
x
)
,
f
2
(
x
)
,
…
,
f
n
(
x
)
)
⊤, where
(39)
f
i
(
x
)
=
x
i
-
sin
|
x
i
-
1
|
,
for
i
=
1,2
,
…
,
n
and C
=
{
x
∈
R
+
n
∣
∑
i
=
1
n
‍
x
i
≤
n
,
x
i
≥
-
1
,
i
=
1,2
,
…
,
n
}. Obviously, Problem 12 is nonsmooth at x
=
(
1,1
,
…
,
1
)
⊤.
The codes are written in Mablab7.0 and run on a personal computer with 2.0 GHZ CPU processor. The parameters used in Algorithm1 are set as ρ
=
0.6, r
=
1
0
-
4, σ
=
5
×
1
0
-
5, and γ
=
1.65. The initial steplength in Step  2 of Algorithm 1 is set to be the spectral coefficient
(40)
β
k
+
1
=
s
k
⊤
s
k
s
k
⊤
z
k
,
where s
k
=
x
k
+
1
-
x
k and z
k
=
F
(
x
k
+
1
)
-
F
(
x
k
)
+
0.01
(
x
k
+
1
-
x
k
). By the monotonicity and the Lipschitz continuity of F
(
·
), it is not difficult to show that
(41)
0.01
s
k
⊤
s
k
≤
s
k
⊤
z
k
≤
(
L
+
0.01
)
s
k
⊤
s
k
,
where L is the Lipschitz constant. If β
k
∉
[
β
min
⁡
,
β
max
⁡
], we replace the spectral coefficient by
(42)
β
k
=
{
1
,
if
∥
F
(
x
k
)
∥
≥
1
,
∥
F
(
x
k
)
∥
-
1
,
if
1
0
-
5
≤
∥
F
(
x
k
)
∥
≤
1
,
1
0
5
,
if
∥
F
(
x
k
)
∥
≤
1
0
-
5
,
where β
min
⁡
=
1
0
-
10 and β
max
⁡
=
1
0
10. This parabolic model is the same as the one described in [15]. We stop the iteration if the iteration number exceeds 1000 or the inequality ∥
F
(
x
k
)
∥
≤
1
0
-
6 is satisfied. The method in [7] (denoted by WPM) is implemented with the following parameters: G
k
≡
0, σ
=
0, λ
=
0.95, β
=
0.5, and μ
k
≡
2.5. The method in [9] (denoted by YSGP) is implemented with the following parameters: β
=
0.5, σ
=
0.01, and r
=
0.001.For Problem11, the initial point is set as x
0
=
(
1,1
,
…
,
1
), and Table 1 gives the numerical results by Algorithm 1, WPM, and YSGP with different dimensions, where Iter. denotes the iteration number, Fn denotes the number of function evaluations, and CPU denotes the CPU time in seconds when the algorithm terminates. Table 2 lists the numerical results of Problem 12 with different initial points. The numerical results given in Tables 1 and 2 show that Algorithm 1 performs a little better than YSGP in [9] and obviously better than WPM in [7], since it requires much lower number of iterations or less CPU time than WPM in [7] and a little lower number of iterations or less CPU time than YSGP in [9]. So the proposed method is promising.Table 1
Numerical results with different dimensions of Problem11.
Dimension
Method
Iter.
Fn
CPU
50
Algorithm1
1
5
0.0000
WPM
614
3703
0.5313
YSGP
6
15
0.0000
500
Algorithm1
1
5
0.0156
WPM
657
3959
1.1406
YSGP
6
15
0.0156
5000
Algorithm1
1
5
0.4219
WPM
692
4157
8.4063
YSGP
7
17
0.4375
50000
Algorithm1
1
5
52.4688
WPM
724
4336
148.9219
YSGP
7
17
53.6719Table 2
Numerical results with different initial points of Problem12 with n
=
64.
Initial point
Method
Iter.
Fn
CPU
(
1,1
,
…
,
1
)
Algorithm1
10
115
0.3438
WPM
330
1981
13.1094
YSGP
16
49
0.6094
(
2,2
,
…
,
2
)
Algorithm1
9
91
0.2031
WPM
331
1985
11.7969
YSGP
16
47
0.6875
(
3,3
,
…
,
3
)
Algorithm1
7
88
0.2031
WPM
331
1987
13.5781
YSGP
16
47
0.6250
(
4,4
,
…
,
4
)
Algorithm1
11
117
0.2500
WPM
331
1987
11.5469
YSGP
16
48
0.6094
(
5,5
,
…
,
5
)
Algorithm1
9
79
0.3281
WPM
331
1987
11.2031
YSGP
16
48
0.6719


## 5. Conclusions

A globally convergent matrix-free method to solve constrained equations has been developed, which is not only derivative-free but also completely matrix-free. Consequently, it can be applied to solve large-scale nonsmooth constrained equations. We established the global convergence without the requirement of differentiability of the equations and presented the linear convergence rate under standard conditions. We also report some numerical results to show the efficiency of the proposed method.Numerical results indicate that the parametersr and γ influence the performance of the method, so the choice of the positive constants r and λ is our future work.



---

*Source: 386030-2014-05-19.xml*