# A Newton-Type Algorithm for Solving Problems of Search Theory

Authors: Liping Zhang
Journal: Advances in Operations Research
 (2013)

---

## Abstract

In the survey of the continuous nonlinear resource allocation problem, Patriksson pointed out that Newton-type algorithms have not been proposed for solving the problem of search theory in the theoretical perspective. In this paper, we propose a Newton-type algorithm to solve the problem. We prove that the proposed algorithm has global and superlinear convergence. Some numerical results indicate that the proposed algorithm is promising.


---

## Body

## 1. Introduction

We consider the problem(1)maxx∈X∑i=1n‍ai(1-exp(-bixi)),
where X={x∈ℝ+n∣eTx=c}, a,b∈ℝ++n, c>0, and e∈ℝn is the vector of ones. The problem described by (1) is called the theory of search by Koopman [1] and Patriksson [2]. It has the following interpretation: an object is inside box i with probability ai, and -bi is proportional to the difficulty of searching inside the box. If the searcher spends xi time units looking inside box i, then he/she will find the object with probability 1-exp(-bixi). The problem described by (1) represents the optimum search strategy if the available time is limited to c time units. Such problems in the form of (1) arise, for example, in searching for a lost object, in distribution of destructive effort such as a weapons allocation problem [3], in drilling for oil, and so forth [2]. Patriksson [2] surveyed the history and applications as well as algorithms of Problem (1); see [2, Sections 2.1.4, 2.1.5, and 3.1.2]. Patriksson pointed out that Newton-type algorithms have not been theoretically analyzed for the problem described by (1) in the references listed in [2].Recently, related problems and methods were considered in many articles, for example, [4–6]. For example, a projected pegging algorithm was proposed in [5] for solving convex quadratic minimization. However, the question proposed by Patriksson [2] was not answered in the literature. In this paper, we design a Newton-type algorithm to solve the problem described by (1). We show that the proposed algorithm has global and superlinear convergence.According to the Fischer-Burmeister function [7], the problem described by (1) can be transformed to a semismooth equation. Based on the framework of the algorithms in [8, 9], a smoothing Newton-type algorithm is proposed to solve the semismooth equation. It is shown that the proposed algorithm can generate a bounded iteration sequence. Moreover, the iteration sequence superlinearly converges to an accumulation point which is a solution to the problem described by (1). Numerical results indicate that the proposed algorithm has good performance even for n=10000.The rest of this paper is organized as follows. The Newton-type algorithm is proposed in Section2. The global and superlinear convergence is established in Section 3. Section 4 reports some numerical results. Finally, Section 5 gives some concluding remarks.The following notation will be used throughout this paper. All vectors are column ones, the subscriptT denotes transpose, ℝn (resp., ℝ) denotes the space of n-dimensional real column vectors (resp., real numbers), and ℝ+n and ℝ++ denote the nonnegative and positive orthants of ℝn and ℝ, respectively. Let Φ′ denote the derivative of the function Φ. We define N:={1,2,…,n}. For any vector x∈ℝn, we denote by diag{xi:i∈N} the diagonal matrix whose ith diagonal element is xi and vec{xi:i∈N} the vector x. The symbol ∥·∥ stands for the 2-norm. We denote by S the solution set of Problem (1). For any α,  β>0,  α=O(β) (resp., α=o(β)) means α/β is uniformly bounded (resp., tends to zero) as β→0.


## 2. Algorithm Description

In this section, we formulate the problem described by (1) as a semismooth equation and develop a smoothing Newton-type algorithm to solve the semismooth equation.We first briefly recall the concepts of NCP, semismooth and smoothing functions [10–12].Definition 1.
A functionϕ:ℝ2→ℝ is called an NCP function if
(2)ϕ(u,v)=0⟺u≥0,v≥0,uv=0.Definition 2.
A locally Lipschitz functionF:ℝn→ℝm is called semismooth at x∈ℝn if F is directionally differentiable at x and for all Q∈∂F(x+d) and d→0,
(3)F(x+d)-F(x)-Qd=o(∥d∥),
where ∂F is the generalized Jacobian of F in the sense of Clarke [13].
F is called strongly semismooth at x∈ℝn if F is semismooth at x and for all Q∈∂F(x+d) and d→0,
(4)F(x+d)-F(x)-Qd=O(∥d∥2).Definition 3.
Letμ≠0 be a parameter. Function Fμ(x) is called a smoothing function of a semismooth function F(x) if it is continuously differentiable everywhere and there is a constant c~>0 independent of μ such that
(5)∥Fμ(x)-F(x)∥≤c~μ,∀x.The Fischer-Burmeister function [7] is one of the well-known NCP functions:
(6)ϕ(u,v)=u+v-u2+v2.
Clearly, the Fischer-Burmeister function defined by (6) is not smooth, but it is strongly semismooth [14]. Let φ:ℝ3→ℝ be the perturbed Fischer-Burmeister function defined by
(7)φ(μ,u,v)=u+v-u2+v2+μ2.
It is obvious that for any μ>0,  φ is differentiable everywhere and for each μ≥0, we have
(8)|φ(μ,u,v)-ϕ(u,v)|≤μ,∀(u,v)∈ℝ2.
In particular, φ(0,u,v)=ϕ(u,v) for all (u,v)∈ℝ2. Namely, φ defined by (7) is a smoothing function of ϕ defined by (6).According to Kuhn-Tucker theorem, the problem described by (1) can be transformed to
(9)eTx=c,xi≥0,s-aibiexp(-bixi)≥0,∀i∈N,xi(s-aibiexp(-bixi))=0,∀i∈N,
where s∈ℝ.Define(10)Ψ(s,x):=(⋮ϕ(xi,s-aibiexp(-bixi))⋮)i∈N,H(s,x):=(eTx-cΨ(s,x)).
According to the Fischer-Burmeister function defined by (6), we formulate (9) as the following semismooth equation:
(11)H(s,x)=0.
Based on the perturbed Fischer-Burmeister function defined by (7), we obtain the following smooth equation:
(12)G(y):=G(μ,s,x):=(μeTx-cΦ(μ,s,x))=0,
where
(13)Φ(μ,s,x):=(⋮φ(μ,xi,s-aibiexp(-bixi))⋮)i∈N.
Clearly, if y*=(0,s*,x*) is a solution to (12) then x* is an optimal solution to the problem described by (1).We give some properties of the functionG in the following lemma, which will be used in the sequel.Lemma 4.
LetG be defined by (12). Then G is semismooth on ℝn+2 and continuously differentiable at any y=(μ,s,x)∈ℝ++×ℝn+1 with its Jacobian
(14)G′(y)=(10000eTΦμ′Φs′Φx′),
where
(15)Φμ′:=vec{-μwi:i∈N},Φs′:=vec{1-(s-aibiexp(-bixi))wi:i∈N},Φx′:=diag{(1-xiwi)+aibi2exp(-bixi)(Φs′)i:i∈N},
with wi:=(xi)2+(s-aibiexp(-bixi))2+μ2, i∈N. Moreover, the matrix G′(y) is nonsingular on ℝ++×ℝn+1.Proof.
G is semismooth on ℝn+2 due to the strong semismoothness of ϕ(u,v). G is continuously differentiable on ℝ++×ℝn+1. For any μ>0, (14) can be obtained by a straightforward calculation from (12). Clearly, we have for any μ>0,
(16)0<1-xiwi<2,0<(Φs′)i<2,aibi2exp(-bixi)>0,∀i∈N,
which implies that (Φx′)ii>0 for all i∈N. Let
(17)(10000eTΦμ′Φs′Φx′)(uvz)=0.
Then, we have u=0 and
(18)∑i=1n‍zi=0,(Φs′)iv+(Φx′)iizi=0,∀i∈N.
The second equality in (18) implies
(19)zi=-(Φs′)i(Φx′)iiv,i∈N.
Since (Φs′)i>0 and (Φx′)ii>0 for i∈N, the first equality in (18) yields v=0 and hence z=0. Therefore, the matrix G′(y) defined by (14) is nonsingular for μ>0.We now propose a smoothing Newton-type algorithm for solving the smooth equation in (12). It is a modified version of the smoothing Newton method proposed in [8]. The main difference is that we add a different perturbed item in Newton equation, which allows the algorithm to generate a bounded iteration sequence. Let y=(μ,s,x)∈ℝn+2 and γ∈(0,1). Define a function ρ:ℝn+2→ℝ+ by
(20)ρ(y):=γ∥G(y)∥min{1,∥G(y)∥}.Algorithm  5Step  0. Choose δ,σ∈(0,1) and μ0>0. Let u-:=(μ0,0,0)∈ℝ++×ℝ×ℝn. Let s0∈ℝ and x0∈ℝn be arbitrary points. Let y0:=(μ0,s0,x0). Choose γ∈(0,1) such that γ∥G(y0)∥<1 and γμ0<1. Set k:=0.Step  1. If G(yk)=0, stop. Otherwise, let ρk:=ρ(yk).Step  2. Compute Δyk:=(Δμk,Δsk,Δxk)∈ℝn+2 by
(21)G(yk)+G′(yk)Δyk=ρku-.Step  3. Let mk be the smallest nonnegative integer such that
(22)∥G(yk+δmkΔyk)∥≤[1-σ(1-γμ0)δmk]∥G(yk)∥.
Let λk:=δmk.Step  4. Set yk+1:=yk+λkΔyk and k:=k+1. Go to Step  1.The following theorem proves that Algorithm  5 is well defined.Theorem 5.
Algorithm  5 is well defined. If it finitely terminates atkth iteration then xk is an optimal solution to the problem described by (1). Otherwise, it generates an infinite sequence {yk=(μk,sk,xk)} with μk>0 and μk≥ρkμ0.Proof.
Ifμk>0 then Lemma 4 shows that the matrix G′(yk) is nonsingular. Hence, Step  2 is well defined at the kth iteration. For any 0<α≤1, define
(23)R(α):=G(yk+αΔyk)-G(yk)-αG′(yk)Δyk.
It follows from (21) that
(24)Δμk=-μk+ρkμ0.
Hence, for any 0<α≤1, we have
(25)μk+αΔμk=(1-α)μk+αρkμ0>0.
From Lemma 4, G is continuously differentiable around yk. Thus, (23) implies that
(26)∥R(α)∥=o(α).
On the other hand, (20) yields
(27)ρk≤γ∥G(yk)∥,ρk≤γ∥G(yk)∥2.
Therefore, for any sufficiently small α∈(0,1],
(28)∥G(yk+αΔyk)∥≤∥R(α)∥+(1-α)∥G(yk)∥+αρkμ0≤(1-α)∥G(yk)∥+αγμ0∥G(yk)∥+o(α)=[1-(1-γμ0)α]∥G(yk)∥+o(α),
where the first inequality follows from (21) and (23), and the second one follows from (26) and (27). Inequality in (28) implies that there exists a constant 0<α~≤1 such that
(29)∥G(yk+αΔyk)∥≤[1-σ(1-γμ0)α]∥G(yk)∥,∀α∈(0,α~].
This inequality shows that Step  3 is well defined at the kth iteration. In addition, by (24), Steps  3 and 4 in Algorithm  5, we have
(30)μk+1=(1-λk)μk+λkρkμ0>0
holds since 0<λk≤1 and μk>0. Consequently, from μ0>0 and the above statements, we obtain that Algorithm  5 is well defined.
It is obvious that if Algorithm  5 finitely terminates atkth iteration then G(yk)=0, which implies that μk=0 and (sk,xk) satisfies (9). Hence, xk is an optimal solution to the problem described by (1).
Subsequently, we assume that Algorithm  5 does not finitely terminate. Let{yk=(μk,sk,xk)} be the sequence generated by the algorithm. It follows that μk>0. We want to prove that {yk} satisfies μk≥ρkμ0 through the induction method. Clearly, ρ(y0)≤γ∥G(y0)∥<1, which yields y0∈Ω. Assume that μk≥ρkμ0; then (24) yields
(31)μk+1-ρk+1μ0=μk+λkΔμk=(1-λk)μk+λkρkμ0-ρk+1μ0≥(1-λk)ρkμ0+λkρkμ0-ρk+1μ0≥μ0(ρk-ρk+1).
Clearly,
(32)ρk={γ∥G(yk)∥2,if∥G(yk)∥<1,γ∥G(yk)∥,otherwise.
It follows from (20) and (22) that
(33)∥G(yk+1)∥≤∥G(yk)∥,ρk+1≤γ∥G(yk+1)∥,ρk+1≤γ∥G(yk+1)∥2.
Hence, combining (31), (32), and (33), we obtain that μk+1≥ρk+1μ0 which gives the desired result.


## 3. Convergence Analysis

In this section we establish the convergence property for Algorithm  5. We show that the sequence{yk=(μk,sk,xk)} generated by Algorithm  5 is bounded and its any accumulation point yields an optimal solution to the problem described by (1). Furthermore, we show that the sequence {yk} is superlinearly convergent.Theorem 6.
The sequence{yk=(μk,sk,xk)} generated by Algorithm  5 is bounded. Let y*=(μ*,s*,x*) denote an accumulation point of {yk}. Then μ*=0, and x* is the optimal solution of the problem described by (1).Proof.
By Theorem5,
(34)μk+1=(1-λk)μk+λkρkμ0≤(1-λk)μk+λkμk=μk,
which implies that {μk} is monotonically decreasing and hence converges. It follows from (22) that
(35)∥G(yk+1)∥≤∥G(yk)∥,∀k.
Hence, {∥G(yk)∥} also converges. Consequently, there exists a constant M>0 such that ∥G(yk)∥≤M for all k. This implies that {μk} and {xk} are bounded, and that for any k and i∈N,
(36)|(xik)2+(sk-aibiexp(-bixik))2+(μk)2xik+(sk-aibiexp(-bixik))(xik)2+(sk-aibiexp(-bixik))2+(μk)2-(xik)2+(sk-aibiexp(-bixik))2+(μk)2|≤M.
This shows that {sk} is also bounded. Consequently, {yk} is bounded. Let y*=(μ*,s*,x*) be an accumulation point of {yk}. We assume, without loss of generality, that {yk} converges to y*. Then, we have
(37)limk→∞∥G(yk)∥=∥G(y*)∥,limk→∞μk=μ*.
By (20),
(38)limk→∞ρk=ρ*:=γ∥G(y*)∥min{1,∥G(y*)∥}.
Suppose ∥G(y*)∥>0 by contradiction. Then ρ*>0. From Theorem 5 and (22), we have
(39)μ*≥ρ*μ0>0,limk→∞λk=0.
Hence, from Step  3 of Algorithm  5, we obtain
(40)∥G(yk+(λkδ)Δyk)∥>[1-σ(1-γμ0)(λkδ)]∥G(yk)∥.
Taking k→∞ in (40) and then combining with (21), we have
(41)-∥G(y*)∥2+ρ*G(y*)Tu-≥-σ(1-γμ0)∥G(y*)∥2,
which yields
(42)ρ*μ0≥[1-σ(1-γμ0)]∥G(y*)∥.
Since ρ*≤γ∥G(y*)∥, (42) implies
(43)(1-σ)(1-γμ0)≤0,
which contradicts σ<1 and γμ0<1. Therefore, ∥G(y*)∥=0 and then μ*=0. Consequently, x* is the optimal solution of the problem described by (1).We next analyze the rate of convergence for Algorithm  5. By Theorem6, we know that Algorithm  5 generates a bounded iteration sequence {yk} and it has at least one accumulation point. The following lemma will be used in the sequel.Lemma 7.
Suppose thaty*=(μ*,s*,x*) is an accumulation point of the iteration sequence {yk} generated by Algorithm  5. Let Q be a matrix in ∂H(s*,x*). Then the matrix Q is nonsingular.Proof.
Define two index sets:(44)N1:={i∈N:(xi*)2+(s*-aibiexp(-bixi*))2≠0},N2:={j∈N:xj*=0,s*=ajbj}.
By a direct computation, we get
(45)Q=(0eTαAβB),
where A is an |N1|×n matrix with all elements being zero except the (i,i)-th as Ai for i∈N1, B is an |N2|×n matrix with all elements being zero except the (j,j)-th as Bj for j∈N2, and
(46)α=vec{1-s*-aibiexp(-bixi*)(xi*)2+(s*-aibiexp(-bixi*))2:i∈N1},Ai=1-xi*(xi*)2+(s*-aibiexp(-bixi*))2+aibi2exp(-bixi*)αi,i∈N1,(βj,Bj)∈{(u,v+ajbj2u):(u-1)2+(v-1)2≤1},j∈N2.
Obviously,
(47)0<αi<2,0<Ai<2(1+aibi2),i∈N1,βj≥0,Bj>0,j∈N2.
Let Q(w,z)=0. Then we have
(48)∑i∈N1‍zi+∑j∈N2‍zj=0;αiw+Aizi=0,i∈N1;βjw+Bjzj=0,j∈N2,
which implies
(49)zi=-αiAiw,i∈N1;zj=-βjBjw,j∈N2.
Therefore,
(50)(∑i∈N1‍αiAi+∑j∈N2‍βjBj)w=0,
which, together with (47), yields w=0. Thus, (49) implies z=0, and hence the matrix Q is nonsingular.Theorem 8.
Let{yk} be the iteration sequence generated by Algorithm  5. Then {yk} superlinearly converges to y*, that is, ∥yk+1-y*∥=o(∥yk-y*∥).Proof.
By Theorem6, {yk} is bounded and then let y*=(μ*,s*,x*) be its any accumulation point. Hence, G(y*)=0 and all matrices Q∈∂G(y*) are nonsingular from Lemma 7. By Proposition 3.1 in [12],
(51)∥G′(yk)-1∥=O(1)
for all yk sufficiently close to y*. From Lemma 4, we know that G is semismooth at y*. Hence, for all yk sufficiently close to y*,
(52)∥G(yk)-G′(yk)(yk-y*)∥=o(∥yk-y*∥).
Since G is locally Lipschitz continuous near y*. Therefore, for all yk sufficiently close to y*,
(53)∥G(yk)∥2=O(∥yk-y*∥2),
which implies that
(54)ρkμ0≤γμ0∥G(yk)∥2=O(∥G(yk)∥2)=O(∥yk-y*∥2).
This inequality, together with (51) and (52), yields
(55)∥yk+Δyk-y*∥=∥yk+G′(yk)-1(-G(yk)+ρku-)-y*∥≤∥G′(yk)-1∥(∥G(yk)-G′(yk)(yk-y*)∥+ρkμ0)=o(∥yk-y*∥).
Following the proof of Theorem 3.1 in [15], we obtain ∥yk-y*∥=O(∥G(yk)∥) for all yk sufficiently close to y*. This implies, from (55), that for all yk sufficiently close to y*,
(56)∥G(yk+Δyk)∥=O(∥yk+Δyk-y*∥)=o(∥yk-y*∥)=o(∥G(yk)∥).
Since limk→∞∥G(yk)∥=0, it follows from (56) that λk=1 can satisfy (22) when yk is sufficiently close to y*. Therefore, for all yk sufficiently close to y*, we have
(57)yk+1=yk+Δyk,
which, together with (55), proves ∥yk+1-y*∥=o(∥yk-y*∥). Namely, {yk} superlinearly converges to y*.


## 4. Computational Experiments

In this section, we report some numerical results to show the viability of Algorithm  5. First, we compare the numerical performance of Algorithm  5 and the algorithm in [5] on two randomly generated problems. Second, we apply Algorithm  5 to solve two real world examples. Throughout the computational experiments, the parameters used in Algorithm  5 were δ=0.75,  σ=0.25,  μ0=0.001, and γ=min{1/∥G(y0)∥,0.99}. In Step  1, we used ∥G(yk)∥≤10-8 as the stopping rule. The vector of ones is all the starting points.Firstly, problems in the form of (1) with 100, 500, 1000, 5000, and 10000 variables were computed. In the first randomly generated example, ai was randomly generated in the interval [10,20], bi was randomly generated in the interval [1,2] for each i∈N, and c was randomly generated in the interval [50,51]. In the second randomly generated example, ai with ∑i=1n‍ai=1 and bi was randomly generated in the interval [0,1] for each i∈N, and c was also randomly generated in the interval [0,1]. Each problem was run 30 times. The numerical results are summarized in Tables 1 and 2, respectively. Here, Dim denotes the number of variable; AT [5] denotes the average run time in seconds used by the algorithm in [5]. In particular, we list more items for Algorithm  5, Inter denotes the average number of iterations, CPU (sec.) is the average run time, and Gval denotes the average values of ∥G(yk)∥ at the final iteration.Table 1
Problem (1) with a, b,c randomly generated in different intervals.
DIM
Inter
Gval
CPU (sec.)
AT [5]
100
7
4.0397
e
-
011
0.0016
0.0008
500
13.6
5.7544
e
-
010
0.0078
0.0021
1000
16.1
9.3739
e
-
010
0.0167
0.0290
5000
22.6
9.2203
e
-
010
0.0844
0.0980
10000
26
9.6499
e
-
010
0.2182
0.3946Table 2
Problem (1) with with a, b,c randomly generated in [0,1].
DIM
Inter
Gval
CPU (sec.)
AT [5]
100
12.3
2.0675
e
-
009
0.006
0.0017
500
15.7
1.3929
e
-
009
0.015
0.0021
1000
16.8
1.2247
e
-
009
0.037
0.0524
5000
19.2
1.0316
e
-
009
0.119
0.1082
10000
21.6
1.9056
e
-
009
0.246
0.4693The numerical results reported in Tables1 and 2 show that the proposed algorithm solves the test problems much faster than the algorithm in [5] when the size of problem is large.Secondly, we apply Algorithm  5 to solve two real world problems. The first example is described in [16]. This problem is how to allocate a maximum amount of total effort, c, among n independent activities, where ai(1-exp(-bixi)) is the return from the ith activity, that is, effort xi, to yield the maximum total return. Note that here ai is the potential attainable and bi is the rate of attaining the potential from effort xi. When no effort is devoted to the ith activity, the value of xi is zero. This example usually arises in the marketing field; the activities may correspond to different products, or the same product in different marketing areas, in different advertising media, and so forth. In this example we wish to allocate one million dollars among four activities with values of ai and bi as in Table 3.Table 3
Activity
A
B
C
D
a
4
×
1
0
6
3
×
1
0
6
2
×
1
0
6
1
0
6
b
2
×
1
0
-
6
3
×
1
0
-
6
1
0
-
6
1
0
-
6For this problem Algorithm  5 obtained the maximum total return4.8952×106 after 25 iterations and elapsing 0.0625 second CPU time. The total effort one million dollars was allocated, 0.5764 million and 0.4236 million to activities A and B, respectively.The second example is to search an object in5 regions, where ai is the prior probability with ∑i=15‍ai=1 of an object of search being in the ith region and 1-exp(-bixi) is the probability of finding an object known to be in the ith region with xi time units. The data of this example, listed in Table 4, come from Professor R. Jiang of Jiaozhou Bureau of Water Conservancy: the available time is c=30 time unite, and after 16 iterations and elapsing 0.0781 second CPU time Algorithm  5 computed the result x=(0,16.6037,0,13.3963,0). This shows that we will spend 17 time units in Region II and 13 time units in Region IV to find water.Table 4
Region
I
II
III
IV
V
a
0.1013
0.3205
0.1323
0.2730
0.1730
b
0.01
0.02
0.01
0.02
0.01


## 5. Conclusions

In this paper we have proposed a Newton-type algorithm to solve the problem of search theory. We have shown that the proposed algorithm has global and superlinear convergence. Some randomly generated problems and two real world problems have been solved by the algorithm. The numerical results indicate that the proposed algorithm is promising.



---

*Source: 513918-2013-01-31.xml*