Finding GM-estimators with global optimization techniques

Journal of Global Optimization 21: 223–237, 2001. © 2001 Kluwer Academic Publishers. Printed in the Netherlands. 223 Finding GM-estimators with global optimization techniques RAFAEL BLANQUERO, EMILIO CARRIZOSA and EDUARDO CONDE Departamento de Estadística e Investigación Operativa, Universidad de Sevilla, Sevilla, Spain (e-mail: [email protected]; [email protected]; [email protected]) Abstract. In this note we address the problem of finding the GM-estimator for the location parameter of a univariate random variable. When this problem is non-convex but d.c. one can use a standard covering method, which, in the one-dimensional case has a simple form. In this paper we exploit the structure of the problem in order to obtain d.c. decompositions with certain optimality properties in the application of the algorithm. Numerical results show that this general-purpose algorithm outperforms previous ad-hoc methods for this problem. Key words: GM-estimators, Robust estimation, D.C. optimization, Covering methods 1. The model Given a sample of nobservations y1,y 2,... ,y nthe determination of a parameter that, in some sense, represents the data is a classical problem in Statistics. In the last decades the traditional least-squares method has been more and more frequently replaced by other approaches with better properties of robustness [19]. In particular, an M-estimator [14, 15], is an optimal solution of an optimization program of the form inf θ∈R 1⩽j⩽n ρ(rj(θ)), (1) where (r1(θ), . . . , rn(θ)) is the vector of residuals, rj(θ) =yj−θj=1,2,... ,n and ρ:R−→ Ris some continuous, even and nondecreasing function in R+. The class of M-estimators has been further enlarged to the class of so-called GM-estimators (generalized M-estimators), in which the influence of each residual is made dependent on the observation yj. In other words, θ∗is said to be a GM- estimator if it solves an optimization problem of the form inf θ∈Rσ(θ), (2) with σ(θ) = 1⩽j⩽n ρj(θ), 224 R.BLANQUERO, E.CARRIZOSA & E.CONDE Table 1. Examples of functions ρ Name ρ(t) 1 Andrews 1−cos tif |t|⩽π 2else 2 Biweight    1−1−t23if |t|⩽1 1 2else 3 Cauchy 1 2log 1+t2 4Fair |t|−log (1+|t|) 5 Huber t2 2if |t|⩽1 |t|−1 2else 6 Logistic log (cosh t) 7p(p ⩾1)|t|p 8Talwar t2 2if |t|⩽1 1 2else 9Welsch 1 21−e−t22 each ρjis typically of the form ρj(t) =ωjρyj−t Cj for predetermined constants ωjand Cj,andρa function as above. A list of typical examples of funtions ρfound in the literature [8, 15, 19], is presented in Table 1. Finding a GM-estimator leads to solve Problem (2) in the case in which the scale parameter is assumed to be known. A much more realistic situation should considerate the estimation of the scale parameter simultaneously with the location coefficient [16]. Usually the scale estimation is computed by means of Eq. (3)  1⩽j⩽n χyj−t s=0(3) where χ:R−→ R. Nevertheless, in some models the estimation process alternates between the computation of the location parameter and the computation of the scale coefficient. For example, in MM-estimation [24] a three-stage procedure FINDING GM-ESTIMATORS 225 is used. In the first stage a tentative approximation to the location parameter is computed assuming a fixed scale coefficient. In the second stage a M-estimate s∗of the errors scale is obtained using, in equation (3), residuals based on the initial estimate. Finally, in the third stage an M-estimate of the location parameter is computed by solving Problem (2) with Cj=s∗for all j=1...n. In this paper, we restrict ourselves to this third stage of the process which is, in general, a difficult task, mainly due to the fact that the functions ρjused in Problem (2) may be non-convex, as is the case of the functions 1, 2, 3, 8, 9 of Table 1. The usual strategy (e.g., [1, 8, 10, 15, 20, 21, 24]) for coping with (2) consists of deriving the normal equation  1⩽j⩽n ρ j(θ) =0(4) a zero of which is sought by means of an iterative procedure (e.g., Newton method). Since (4) is not a sufficient condition for optimality in the non-convex case, such iterative procedures, strongly dependent of the initial guess, may be trapped in the neighborhood of a zero of (4) which does not represent a global minimum of (2). When the problem is non-convex one has to choose between seeking a local minimum (close to the initial value, thus pressumably retaining some of its properties) or a global minimum. For the latter case, one can use well-known numerical procedures such as the covering method, first proposed by Piyavskii, [17], see also [3, 4, 5]. An outline of the algorithm is given below. ALGORITHM 1. Initialization: Find a compact interval R∗known to contain an optimal solution of (2). Set BEST.THETA equal to an arbitrary θ1∈R∗and set BEST.VA LU E := σ(BEST.THETA). Construct E(1)and set k=1. Step k:Set LBOUND =min θ∈R∗E(k)(θ) (5) and let θk+1be its optimal solution. If BEST.VA L U E −LBOUND ⩽εthen GoTo Output. Else, do: 1. Construct E(k+1)from E(k). 2. Update, if required, BEST.VA LU E and BEST.THETA 3. GoTo step k+1. Output: STOP: BEST.THETA is an ε-optimal solution. The three critical issues of the algorithm are 226 R.BLANQUERO, E.CARRIZOSA & E.CONDE 1. Building R∗. 2. Defining the lower envelopes Ekof σin order to guarantee convergence. 3. Solving the auxiliary problems (5). Since ρis assumed to be even and non-decreasing in R+, we can take as R∗the interval given by the extreme observations, R∗=min j=1,... ,n yj,max j=1,... ,n yj.(6) In what follows we customize Algorithm 1 for the one-dimensional problem (2) under the extra assumption that each ρjis the difference of two convex functions, i.e., it is d.c., [11, 12, 22]. Hence, by the algebra of d.c. functions, σis d.c. Let σ=σ1−σ2bead.c. decomposition of σ,andletσ 1(θi;θ−θi)be the directional derivative of σ1at θi in the direction θ−θi, i.e., σ 1(θi;θ−θi)=lim t↓0 σ1(θi+t(θ −θi)) −σ1(θi) t Let hibe given as hi(θ) =σ1(θi)+σ 1(θi;θ−θi)−σ2(θ), (7) and define E(k) as the pointwise maximum of the functions hi, E(k) =max 1⩽i⩽khi Observe that this slightly differs from the standard choice  hi,  hi(θ) =σ1(θi)+ηi(θ −θi)−σ2(θ), (8) with ηi∈∂σ1(θi), the subdifferential of σ1at θi. Since σ(θi;θ−θi)=max{η(θ −θi):η∈∂σ1(θi)}, [18], we see that hi⩾ hi.Moreover, such a choice will yield certain optimality properties, see Section 2. We will show in this section that, with this choice, Problem (5) is solved after inspecting at most two easily located points, and the updating of E(k) is straightforward. Moreover, the finiteness of the algorithm directly follows from [4, 17], as stated in Proposition 1 PROPOSITION 1. For all k=1,..., E(k) ⩽E(k+1)⩽σ. FINDING GM-ESTIMATORS 227 Moreover for i=1,... ,k, E(k)(θi)=σ(θ i). If BEST.VA LU E and LBOUND are as above, then: (a) lim k↑∞ BEST.VA LU E (k) =min θ∈R∗σ(θ) (b) lim k↑∞(BEST.VA LU E (k) −LBOUND(k)) =0 We discuss now how to solve (5), showing that finding an optimal solution of (5) and updating E(k) can be done in a straightforward manner. PROPOSITION 2. Let k⩾i. The sets Ik i+:= {θ∈R∗|hi(θ) ⩾hs(θ), θ ⩾θi,∀s=1,... ,k} Ik i−:= {θ∈R∗|hi(θ) ⩾hs(θ), θ ⩽θi,∀s=1,... ,k} are intervals with θi∈Ik i+∩Ik i−. Proof. First, θi∈Ik i+∩Ik i−since hi(θi)=σ(θ i)⩾E(k)(θi)⩾hs(θi)for all s. Now, observe that σ 1(θs;θ−θs)is convex in θ, with σ 1(θs;θ−θs)=σ 1+(θs)(θ −θs), if θ ⩾θs σ 1−(θs)(θ −θs), if θ ⩽θs , σ 1+and σ 1−denoting the right and left derivatives of σ1. We have {θ⩾θi|hi(θ) ⩾E(k)(θ)}= ={θ⩾θi|σ1(θi)+σ 1(θi;θ−θi)⩾max 1⩽s⩽k σ1(θs)+σ 1(θs;θ−θs)} ={θ⩾θi|σ1(θi)+σ 1+(θi)(θ −θi)⩾max 1⩽s⩽k σ1(θs)+σ 1(θs;θ−θs)} ={θ⩾θi|0⩾max 1⩽s⩽kσ1(θs)+σ 1(θs;θ−θs)−σ1(θi)−σ 1+(θi)(θ −θi)} which is an interval since it is the lower level set of a convex function. Hence Ik i+=R∗∩{θ⩾θi|hi(θ) ⩾E(k)(θ)}is an interval. In the same manner Ik i−is another interval, with θias endpoint.  From this we have PROPOSITION 3. E(k) is a piecewise concave function, having the intervals Ik 1+,Ik 1−, ...Ik k+,Ik k−as domains of concavity. 228 R.BLANQUERO, E.CARRIZOSA & E.CONDE REMARK 1. Proposition 2 enables us to easily update E(k+1)from E(k),since they only differ in one interval, namely Ik+1 k+1=Ik+1 k+1,−∪Ik+1 k+1,+: E(k+1)(θ) =   E(k)(θ), if θ /∈Ik+1 k+1 σ1(θk+1)+σ 1+(θk+1)(θ −θk+1)−σ2(θ), if θ ∈Ik+1 k+1,+ σ1(θk+1)+σ 1−(θk+1)(θ −θk+1)−σ2(θ), if θ ∈Ik+1 k+1,− Hence, in order to solve (5) we have min θ∈R∗E(k)(θ) =min 1⩽s⩽k{min θ∈Ik s− hs(θ), min θ∈Ik s+ hs(θ)}. This joint with Proposition 1 imply the following PROPOSITION 4. At Step kof Algorithm 1 the optimal solution θk+1of Problem (5) either coincides with θk(in which case it is optimal for Problem (2)), or it is an endpoint of an interval Ik i:= Ik i−∪Ik i+. REMARK 2. By Proposition 4, in order to minimize the new envelope function E(k+1)we only need to inspect two new points, namely the endpoints of the interval Ik+1 k+1. 2. Constructing the d.c. decomposition Solving Problem (2) by using Algorithm 1 requires the knowledge of a d.c. representation for the objective function. Although obtaining explicitly a decomposition for an arbitrary d.c. function may be far from trivial, this task turns out to be very simple in our case. Indeed, using the standard techniques of [12, 22], d.c. decompositions for all the non-convex functions ρof Table 1 can be directly derived. For instance, for Function 8 in Table 1, we have ρ(t) =min{t2 2,1 2} =t2 2+min{0,1 2−t2 2} =t2 2−max{0,t2 2−1 2} On the other hand, for Function 3 in Table 1, we have that ρ(t) =1−t2 (1+t2)2, FINDING GM-ESTIMATORS 229 thus ρ(t) ⩾min 0⩽s 1−s (1+s)2=−1 8, yielding the d.c. decomposition of ρ ρ(t) =ρ(t) +1 2 t2 8−1 2 t2 8(9) A similar analysis can be carried out for the remaining non-convex functions ρ of Table 1. Once a d.c. decomposition of ρ, ρ(t) =α(t) −β(t) is given, one immediately obtains a d.c. decomposition for ρj, ρj(θ) =αyj−θ Cj−βyj−θ Cj, and then a d.c. decomposition for σis given by σ(θ) = 1⩽j⩽n αyj−θ Cj− 1⩽j⩽n βyj−θ Cj For a given function ρ, more than one d.c. decomposition may be found. In spite of the fact that the choice of the decomposition influences the speed of convergence of the algorithm, see [3, 4] and Section 3, this is usually ignored in the literature. Here, good d.c. decompositions are sought, and we will prove that, under certain conditions, an optimal (in a sense described later) decomposition for this particularization of Algorithm 1 can be obtained. This result will be used for finding an optimal d.c. decomposition for every term in (2) and, from here, a heuristic good representation for σ. Let us consider the set Dof real d.c. finite functions defined over an open interval of Rthat contains the closed bounded interval Iwhere the optimization will be carried out. We recall the reader that any function f∈Dhas side derivatives and these are of bounded variation on I,[7]. The following result is given in [9]: LEMMA 1. Let f∈D,ˆ θ∈I, and let Nˆ θ(t) denote the negative variation of f −, the left derivative of f, in the interval with endpoints ˆ θand t. Then the functions ψ[ˆ θ](θ) := f(θ)+θ ˆ θ Nˆ θ(t)dt γ[ˆ θ](θ) := θ ˆ θ Nˆ θ(t)dt are convex. In particular, ψ[ˆ θ]−γ[ˆ θ]is a d.c. decomposition of f. 230 R.BLANQUERO, E.CARRIZOSA & E.CONDE Our aim is to use the decomposition given in Lemma 1 for the functions ρjin Problem (2). The suitability of such a choice in Algorithm 1 is motivated in the following results. First we show that the choice of ˆ θis not relevant, since we always obtain the same lower envelope. LEMMA 2. Let f∈Dand θ0,ˆ θ∈I.Letψ[ˆ θ]and γ[ˆ θ]be defined as in Lemma 1. Then, the function θ∈I−→ ψ[ˆ θ](θ0)+ψ[ˆ θ] −(θ0)(θ −θ0)−γ[ˆ θ](θ) does not depend on ˆ θ. Proof. First of all, note that ([18], Theorem 24.2): ψ[ˆ θ](θ0,θ −θ0)=f −(θ0)(θ−θ0)+N− ˆ θ(θ0)(θ−θ0)if θ<θ0 f +(θ0)(θ−θ0)+N+ ˆ θ(θ0)(θ−θ0)if θ⩾θ0 where f +and f −are respectively the right and left derivatives of fand N− ˆ θ(θ0)=lim θ↑θ0 Nˆ θ(θ) N+ ˆ θ(θ0)=lim θ↓θ0 Nˆ θ(θ) Given θ<θ 0, one has that: ψ[ˆ θ](θ0)+ψ −[ˆ θ](θ0)(θ −θ0)−γ[ˆ θ](θ) =f(θ 0)+θ0 ˆ θ Nˆ θ(ξ)dξ +f −(θ0)(θ−θ0)+N− ˆ θ(θ0)(θ−θ0) −θ ˆ θ Nˆ θ(ξ)dξ =f(θ 0)+f −(θ0)(θ−θ0)−θ θ0 Nˆ θ(ξ)dξ +N− ˆ θ(θ0)(θ−θ0) =f(θ 0)+f −(θ0)(θ−θ0)−θ θ0 {Nˆ θ(ξ) −N− ˆ θ(θ0)}dξ =f(θ 0)+f −(θ0)(θ−θ0)−θ θ0 {Nθ0(ξ) −N− θ0(θ0)}dξ =f(θ 0)+f −(θ0)(θ−θ0)+N− θ0(θ0)(θ −θ0)−θ θ0 Nθ0(ξ)dξ =ψ[θ0](θ0)+ψ −[θ0](θ0)(θ −θ0)−γ[θ0](θ) since Nˆ θ(ξ) −N− ˆ θ(θ0)=Nˆ θ(ξ) +Nˆ θ(θ0)−Nˆ θ(θ0)−N− ˆ θ(θ0) =Nθ0(ξ) +Nˆ θ(θ0)−N− ˆ θ(θ0) =Nθ0(ξ) −N− θ0(θ0) FINDING GM-ESTIMATORS 231 In the same manner, one can cope with θ⩾θ0,so that the function θ∈I−→ ψ[ˆ θ](θ0)+ψ[ˆ θ] −(θ0)(θ −θ0)−γ[ˆ θ](θ) does not depend on ˆ θ. The following result establishes the optimal behaviour of the decomposition of Lemma 1 for applying Algorithm 1, in the following sense: given θ0∈I,the bounding function that it provides is a pointwise majorant for the bounding function obtained from any other d.c. decomposition. In this way, the minimum value of the envelope E(k) provided for this d.c. decomposition will be greater or equal than the minimum of the envelope built from any other representation, so the algorithm should converge more quickly. PROPOSITION 5. Let f∈Dand θ0∈I.Letψ−γbe the d.c. decomposition of Lemma 1 and τ−λany other d.c. decomposition of f. Then, one has for all θ∈I ψ(θ0)+ψ −(θ0)(θ −θ0)−γ(θ)⩾τ(θ0)+τ −(θ0)(θ −θ0)−λ(θ). Proof. Let θ<θ 0. Then, since f=τ−λ, with τand λconvex, f −is the difference of two non-decreasing functions, namely f −=τ −−λ −. Making use of the properties of the bounded variation functions, [7], one has that: Nθ0(ξ) ⩽Tθ0(ξ) =λ −(ξ) −λ −(θ0) where Nθ0(ξ) denotes the negative variation of f −and Tθ0(ξ) is the total variation of λ −in the interval [θ0,ξ].Sinceλis convex, integrating the previous expression between θ0and θ, it follows that ([18], Corollary 24.2.1): θ θ0 Nθ0(ξ)dξ ⩽λ(θ) −λ(θ0)−λ −(θ0)(θ −θ0) from where, since λ=τ−f, we obtain: f(θ 0)+f −(θ0)(θ−θ0)−θ θ0 Nθ0(ξ)dξ ⩾τ(θ0)+τ −(θ0)(θ −θ0)−λ(θ) A similar analysis can be carried out for the case θ⩾θ0, yielding the result.  The following result establishes the optimality of a classical d.c. decomposition for certain types of functions. Indeed, the d.c. decomposition for a C2-function proposed in [2], is a particular case of the representation of Lemma 1 and, using Proposition 5, it is optimal for applying Algorithm 1.