An exact global optimization method for deriving weights from pairwise comparison matrices

J Glob Optim (2007) 38:237–247 DOI 10.1007/s10898-006-9073-5 ORIGINAL PAPER An exact global optimization method for deriving weights from pairwise comparison matrices Emilio Carrizosa ·Frédéric Messine Received: 12 January 2006 / Accepted: 27 July 2006 / Published online: 6 September 2006 © Springer Science+Business Media B.V. 2006 Abstract Some multiple-criteria decision making methods rank actions by associating weights to the different criteria or actions, which are pairwise compared via a positive reciprocal matrix A. There is a vast literature on proposals of different mathematical-programming methods to infer weights from such matrix A. However, it is seldom observed that such optimization problems may be multimodal, thus the standard local-search resolution techniques suggested may be trapped in local optima, yielding a wrong ranking of alternatives. In this note we show that standard tools of global optimization based on interval analysis, lead to globally optimal weights in reasonable time. Keywords Branch and bound ·Interval analysis ·Multiple-criteria decision making ·Pairwise comparison matrices 1 Introduction In order to rank a set of ngiven decisions, a possible strategy consists of associating with them weights x1,...,xnthrough a nonlinear regression model of the form xi xj =aij,i,j=1, 2, ...,n,(1) where the elements aij, representing the relative preference of the ith action over the jth action, are positive constant assumed to be given. E. Carrizosa Facultad de Matemáticas, Universidad de Sevilla, Tarfia s/n 41012, Sevilla, Spain e-mail: [email protected] F. Messine (B ) ENSEEIHT-IRIT, 1 rue C Camichel 31071, Toulouse Cedex, France e-mail: F[email protected] 238 J Glob Optim (2007) 38:237–247 Different numerical procedures, not absent of controversy, have been proposed in the literature to infer xfrom model (1) (see e.g. [1, 2, 8, 12, 14, 16–18] and the references therein). A broad class of such procedures yields as xan optimal solution of a problem of the form min x∈Rn ++ γ|xi/xj−aij|n i,j=1,(2) where R++ denotes the set of strictly positive reals, and γis a monotonic norm in the non-negative orthant Rn×n +, i.e., γ(u)≤γ(v),∀u,v,0≤u≤v(3) such as any (weighted) pnorm, 1 ≤p≤+∞. This approach has been proposed, amon others, in [2, 5, 9, 11]. In particular, in [5] it is claimed that “such methods are more appealing because it is undisputable that the matrix (xi/xj)associated to the solution vector xmust be close to the judgment matrix A.” Three particular instances of (2) are included in [5, 11], namely the unweighted 1,2and ∞, under the names of Least Absolute Error,Least Square and Least Worst Absolute Error (see e.g. [16–18]) for a different strategy to infer weights from Awhich is not based on optimization techniques and may lead, as shown in [2], to non-efficient weights. In spite of the fact that (2) may be multimodal, local-search methods seem to be the only proposal so far, the exception being the recent works [3, 4, 7], where different global-optimization strategies for solving (2) for the Euclidean norm γare analysed. The methods in these references exploit extensively the properties of the Euclidean norm, and thus cannot be easily extended to more general settings, as the one considered in this paper. Forchoicesofγdifferent of the Eculidean norm, the literature is rather scarce. In [5], for the particular instances of (2) addressed, (1and ∞norm together with the 2), the problems are considered to be “difficult to solve”, although no attempt to solving them is given. These three problems are later addressed in [11] via “a simple genetic algorithm”, which may yield solutions far from optimality. From a practical viewpoint, it has been observed that it is not so easy to provide precise values for the scalar aij, which are thus replaced by intervals Aij =[aL ij ,aU ij ] (see e.g. [6, 10, 19, 20] and the references therein). In this case, the regression model (1) is replaced by the more general model xi xj ∈Aij,i,j=1, 2, ...,n(4) and (2) is replaced by the problem of finding strictly positive weights xjminimizing a norm γof the distances between the ratios xi/xjand the corresponding intervals Aij, min x∈Rn ++ f(x):=γεxi/xj,Aijn i,j=1(5) J Glob Optim (2007) 38:237–247 239 with ε(s,[aL,aU])defined as the distance between sand the closest point in [aL,aU], ε(s,[aL,aU])=⎧ ⎨ ⎩ 0, if s∈[aL,aU], aL−s,ifs≤aL, s−aU,ifs≥aU (6) =max 0,  s−aU+aL 2 −aU−aL 2.(7) Obviously (5) coincides with (2) in what we call the degenerate case in which all intervals Aij are degenerate, Aij ={aij}. The aim of this paper is to show that well-known strategies of deterministic global optimization, namely, branch and bound algorithms with bounds based on interval analysis, provide in reasonable time globally optimal solutions to problem (5), and thus, as particular case, to model (2), for different choices of γ.Hence,ifa decision-maker infers weights from a pairwise comparison matrix following model (5) for a given norm γ, we show that, although, according to [5], the problems are “difficult to solve”, they are solvable exactly by standard Global Optimization methods. At the same time, uncertainties in the comparisons, modelled as intervals Aij for the weights ratios, are included naturally. This is an interesting advantage of model (5) against, for instance, the methods recommended in [5, 11, 16]. In what follows γis an arbitrary norm monotonic in Rn×n +, as defined in (3). We stress that γis assumed to be given, thus we do not make any suggestion about the choice of γ. Instead, we refer to [5], where it is said that “no one method is expected to outperform other methods at all situations.” For i,j=1, ...,n,Aij =[aL ij ,aU ij ]is an interval satisfying 0 <aL ij ≤aU ij . Again, this interval matrix is assumed to be given. With these elements at hand, we now show how to infer the weights from such (interval) pairwise comparison matrix. 2 General results According to its definition (5), f(λx)=f(x)for any x∈Rn ++ and any choice of λ>0. Hence, a normalization constraint can be added to (5) without loss of generality. In particular, we can fix x1to unity, and replace (5)by min f(x):=γεxi/xj,Aijn i,j=1, s.t. x1=1, x∈Rn ++. (8) Observe that this way, we lower in one the dimensionality of the problem, which may be of great importance when deterministic global optimization techniques are to be used. Since norms are continuous, we have by (7) that the objective function (8) is continuous in the feasible region. However, such feasible region is not closed, which makes at first glance unclear whether the optimal value is attained. However, one can easily construct a compact interval containing all optimal solutions of (8), which are shown to exist. Indeed, one has 240 J Glob Optim (2007) 38:237–247 Proposition 1 Let z be strictly greater than the optimal value of (8). Then, any optimal solution x to (8)satisfies xj∈Xj:=[xL j,xU j],j=1, 2, ...,n, where Xj=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ {1},ifj=1, ⎡ ⎣max i<j xL i aU ij +z γ(eij) ,min i<jxU iaU ji +z γ(eji)⎤ ⎦,ifj>1(9) and eij is the n ×n matrix with all elements null except the element (i,j),with value 1. In particular, (8), and thus (5), always has optimal solutions. Proof We will show that for any x, feasible for (8), but not satisfying (9), f(x)is at least as great as z, known to be strictly greater than the optimal value. Hence, x cannot be optimal. Having shown this, we have in particular that (8) is equivalent to the problem where the constraint x∈Rn ++ is replaced by (9), and thus we are optimizing a continuous function over a compact domain, guaranteeing then the existence of optimal solutions. Since, by assumption, γis monotonic in Rn×n +, one has f(x)≥γ(ε(xi/xj,Aij)eij),∀i,j,∀x∈Rn ++. (10) Take xnot satisfying (9), and let jbe the smallest index for which xj∈ Xj. Then, there must exist i<jsuch that either xj<xL i aU ij +z γ(eij) (11) or xj>xU iaU ji +z γ(eji). (12) If (11) holds, then, ε(xi/xj,Aij)=xi xj −aU ij >z γ(eij). By (10), f(x)≥γ(ε(xi/xj,Aij)eij) ≥γz γ(eij)eij=z thus xcannot be optimal. Analogously, if (12) holds then ε(xj/xi,Aji)=xj xi −aU ji >z γ(eji). J Glob Optim (2007) 38:237–247 241 Hence, f(x)≥γ(ε(xj/xi,Aji)eji) ≥γz γ(eji)eji=z, showing that such xcannot be optimal.  Following Proposition 1, we can reformulate (8) as a nonlinear nonconvex boxconstrained problem of the form min f(x):=γεxi/xj,Aijn i,j=1, s.t. xj∈Xj,j=1, 2, ...,n.(13) Let us now discuss how to come up with a globally optimal solution to (13), first in the particular case in which the norm γis the ∞norm, and then for arbitrary monotonic norms γ. 2.1 The ∞norm For γequal the ∞norm, we have by (7)that f(x)=max 1≤i,j≤nε(xi /xj,Aij) =max 1≤i,j≤nmax 0,  xi xj − aU ij +aL ij 2 − aU ij −aL ij 2 =max 0, max 1≤i,j≤n xi xj − aU ij +aL ij 2 − aU ij −aL ij 2. It is then seen that fis quasiconvex and strictly quasiconvex on jXj, thus local optima are global optima. Moreover, the search of (globally) optimal solutions can be reduced to solving a series of linear programs. Indeed, a simple calculation shows that, for β>0andx∈jXj, the statement f(x)<βis equivalent to the statement 0<min ij min xj(aU ij +β) −xi,xi−xj(aL ij −β). Hence, for δ>0 given, a δ-approximate optimal solution xwill be obtained after solving (via a binary search), O(log(1/δ ) ) linear problems of the form max t s.t. t≤xj(aU ij +β) −xi,∀i,j, t≤xi−xj(aL ij −β),∀i,j, xj∈Xj,∀j. Observe that this problem had been (erroneously) claimed to be “difficult” in [5, 11]. 2.2 Arbitrary monotonic norms It is not hard to obtain (convergent) bounds in fwithin intervals using basic tools of interval analysis [13]. Indeed, for intervals S=[sL,sU],T=[tL,tU], it is easily 242 J Glob Optim (2007) 38:237–247 shown, by enumeration of all possible cases, that the range ε(s,t):={ε(s,T):s∈S} is given by ε(s,t)=[α(s,t),β(s,t)] with α(s,t)=max tL−sU,sL−tU,0 , β(s,t)=max sU−tU,tL−sL,0 . Hence, by the monotonicity of γ, we have that F((X1,...,Xn)) :=γ(α(Xi/Xj,Aij))n i,j=1,γ(β(Xi/Xj,Aij))n i,j=1 (14) is an inclusion function for f. Moreover, since all norms are equivalent, it turns out that Fexhibits first-order convergence in the domain of (13). Hence, a branch and bound algorithm using γ(α(Xi /Xj,Aij))n i,j=1as lower bound of fon intervals is convergent. 3 Numerical experiments A series of tests has been performed to check the feasibility of interval branch- and-bound methods to solve (8) and to explore how perturbations in model (1) affect accuracies and computing times. For a fixed dimension n, we generate a weight vector x=(x1,...,xn)randomly, with components independently and identically distributed as a uniform variable in the set {1, 2, ...,9}; for each pair i,j,i<j, a perturbation factor ξij is then randomly generated following a uniform distribution in an interval [1−c,1+c]with c∈R+, and the interval pairwise comparison matrix in model (4)is defined with Aij =⎧ ⎨ ⎩ [(1−r)xi xjξij,(1+r)xi xjξij],ifi<j, {1},ifi=j, 1/Aji,ifi>j, (15) where ris a fixed positive real number which denotes the wished radius for the intervals Aij,i<j:i.e.2ris the width of the intervals Aij for i<j. For each size nand constant cdelimiting the perturbations ξij in the interval [1−c,1+c], 50 interval pairwise matrices with different widths (defined by their radius r) were generated according to (15). The analysis has been performed using as γthe Euclidean norm; exactly the same methodology can be used for other norms or other objective functions. Mean, best and worst-case results are shown in Tables 1–3. The used algorithm derives from a standard interval branch-and-bound algorithm of Moore–Skelboe type, where the midpoint test is added (this generates the current solution named ˜ fin Algorithm 1), [15]. These kinds of exact global optimization methods work by bisecting the initial domain into smaller and smaller subboxes. Thus, a list of boxes (named Lin Algorithm 1) is generated; boxes are discarded as soon as it is shown that they cannot contain a global optimum (computing a lower bound of a function over a box using interval analysis [13] it can be proved that the best solution in the box cannot be superior to the current one ˜ f). The stopping criterion for a standard interval branch-and-bound algorithm is generally the accuracy reached between the current minimum ˜ fand the lowest lower bound remaining in L(˜ f−min(z,Z)∈Lz<e, J Glob Optim (2007) 38:237–247 243 where eis the wished precision). The bounds are computed by using the inclusion function defined by Eq. (14) and using the definition of interval arithmetic, [13]. For all our numerical tests (2, 250), it was not possible to fix the same accuracy for all of them, because the values of the objective functions vary from 0 up to 100 depending on the size and of the perturbation of the problem (parameters n,r,c). Thus, the algorithm was modified in order to compute the efficient accuracy which can be reached by the optimal solutions (within a fixed number of iterations NBmax). The parameter co represents the step by which the precision eis decreased during the iterations of the algorithm (it could be fixed to 0.0001). The maximal number of iterations NBmax is fixed to 200, 000 in Algorithm 1in order to solve quickly all the numerical tests presented in Tables 1–3. Therefore one has to discuss about the precision reached during this imposed time. This method is detailed as follows: Algorithm 1 •Set X:=the initial domain in which the global minimum is searched, X⊆Rn. •Set ˜ f:=+∞ (the current minimum). •Set L:=(+∞,X)(the lower bound, the box). •Set e:=10 (the wished accuracy). •Set co :=1. •Set NBmax :=200000. •Set k:=1. •While k≤NBmax and e≥10−4Do (1) Extract from Lthe box which has the lowest lower bound. (2) Bisect the considered box by the middle of its edge which has the maximal length, yielding V1,V2. (3) For m:=1to2do *Compute vm:=γα(Vm)i (Vm)j ,Aijijthe lower bound of f on Vm,by using (14). *if vm<˜ f(only one global solution is sought) then .Insert (vm,Vm)in L. .Set ˜ f:=min(˜ f,f(mid(Vm))), where mid(Vm)is the midpoint of Vm. .If ˜ fhas changed then remove from Lall (z,Z) where z≥˜ f,and set ˜ y:=mid(Vm). (4) if ˜ f−min(z,Z)∈Lz<ethen *if co ≥ethen co :=co/10. *e:=e−co. *Compute the CPU-time in seconds until this instruction. (5) k:=k+1. End Do •Results: ˜ f,˜ y,e,k,time. The results given at the end of the algorithm, are the global minimum value ˜ f, only one solution ˜ yand the reached accuracy e, the CPU-time time, the number of iterations 244 J Glob Optim (2007) 38:237–247 Table 1 Mean time and mean precision for 50 problems Pb size n=5n=6n=7n=8n=9n=10 r,ctime (s) Prec time (s) Prec time (s) Prec time (s) Prec time (s) Prec time (s) Prec 0, 0.1 551.83 4.09 ×10−3598.08 1.82 ×10−2521.24 7.19 ×10−2409.94 1.48 ×10−1304.73 2.28 ×10−1211.56 4.55 ×10−1 0, 0.2 573.74 1.67 ×10−2471.04 1.07 ×10−1415.63 3.32 ×10−1324.38 8.10 ×10−1313.54 1.18 ×100327.89 1.67 ×100 0, 0.3 585.75 2.23 ×10−2486.93 2.28 ×10−1451.43 6.17 ×10−1473.96 1.41 ×100439.37 1.89 ×100188.35 2.91 ×100 0.05,0.1 21.62 1.08 ×10−463.33 2.31 ×10−4151.98 5.39 ×10−4203.64 1.56 ×10−3304.38 2.83 ×10−3228.84 3.46 ×10−3 0.05, 0.2 463.63 2.34 ×10−3571.14 1.22 ×10−2406.58 5.27 ×10−2414.92 1.07 ×10−1292.56 2.17 ×10−1294.40 2.88 ×10−1 0.05, 0.3 547.99 8.17 ×10−3451.54 4.78 ×10−2475.70 2.17 ×10−1422.24 3.59 ×10−1243.99 8.16 ×10−1295.19 1.11 ×100 0.1, 0.1 0.03 1.00 ×10−41.37 1.00 ×10−40.44 1.00 ×10−44.60 1.43 ×10−421.42 1.84 ×10−411.32 1.05 ×10−4 0.1, 0.2 0.85 1.00 ×10−451.34 2.15 ×10−4124.24 8.86 ×10−451.34 2.15 ×10−4153.29 2.36 ×10−3158.38 6.18 ×10−3 0.1, 0.3 315.51 1.63 ×10−3485.98 1.14 ×10−2286.86 2.63 ×10−2426.37 5.56 ×10−2298.92 1.33 ×10−1256.64 1.92 ×10−1 Table 2 Best time and best precision for 50 problems Pb size n=5n=6n=7n=8n=9n=10 r,ctime (s) Prec time (s) Prec time (s) Prec time (s) Prec time (s) Prec time (s) Prec 0, 0.1 71.63 1.20 ×10−438.00 4.00 ×10−30.40 3.00 ×10−33.36 3.00 ×10−22.30 5.00 ×10−20.19 1.00 ×10−1 0, 0.2 66.74 5.00 ×10−46.81 5.00 ×10−31.58 4.00 ×10−20.21 9.00 ×10−20.45 2.00 ×10−11.19 3.00 ×10−1 0, 0.3 80.83 7.00 ×10−49.97 8.00 ×10−34.52 7.00 ×10−22.93 3.00 ×10−12.94 4.00 ×10−10.80 5.00 ×10−1 0.05, 0.1 0.01 1.00 ×10−40.02 1.00 ×10−40.05 1.00 ×10−40.10 1.00 ×10−40.27 1.00 ×10−40.37 1.00 ×10−4 0.05, 0.2 0.02 1.00 ×10−450.51 1.00 ×10−44.62 1.00 ×10−313.30 7.00 ×10−30.50 2.00 ×10−22.41 5.00 ×10−2 0.05, 0.3 168.80 2.00 ×10−31.86 1.00 ×10−460.97 2.00 ×10−36.79 2.00 ×10−27.61 2.00 ×10−24.00 1.00 ×10−1 0.1, 0.1 0.01 1.00 ×10−40.02 1.00 ×10−40.03 1.00 ×10−40.05 1.00 ×10−40.07 1.00 ×10−40.09 1.00 ×10−4 0.1, 0.2 0.02 1.00 ×10−40.02 1.00 ×10−40.04 1.00 ×10−40.07 1.00 ×10−40.51 1.00 ×10−40.64 1.00 ×10−4 0.1, 0.3 0.02 1.00 ×10−40.09 1.00 ×10−40.04 1.00 ×10−45.19 2.00 ×10−30.37 4.00 ×10−30.74 2.00 ×10−2 J Glob Optim (2007) 38:237–247 245 Table 3 Worst time and worst precision for 50 problems Pb size n=5n=6n=7n=8n=9n=10 r,ctime (s) Prec time (s) Prec time (s) Prec time (s) Prec time (s) Prec time (s) Prec 0, 0.1 1504.04 2.00 ×10−21102.12 1.00 ×10−11167.66 4.00 ×10−11477.52 5.00 ×10−11233.16 6.00 ×10−11015.84 2.00 ×100 0, 0.2 1110.23 2.00 ×10−11135.60 4.00 ×10−11115.34 2.00 ×1001468.28 3.00 ×1001102.48 4.00 ×1001231.78 3.00 ×100 0, 0.3 1081.01 9.00 ×10−21030.16 2.00 ×1001162.78 2.00 ×1001204.36 5.00 ×1001219.34 5.00 ×1001063.33 8.00 ×100 0.05, 0.1 1034.59 5.00 ×10−4991.14 3.00 ×10−31073.25 6.00 ×10−31113.44 2.00 ×10−21007.06 2.00 ×10−2926.07 3.00 ×10−2 0.05, 0.2 1048.06 3.00 ×10−21095.47 9.00 ×10−21114.07 3.00 ×10−11284.31 7.00 ×10−11112.54 1.00 ×100944.01 1.00 ×100 0.05, 0.3 1077.03 6.00 ×10−21140.27 4.00 ×10−11135.59 2.00 ×1001208.64 2.00 ×1001099.78 3.00 ×1001126.40 3.00 ×100 0.1, 0.1 0.18 1.00 ×10−44.93 1.00 ×10−48.01 1.00 ×10−4182.71 2.00 ×10−3617.16 3.00 ×10−3349.49 3.00 ×10−4 0.1, 0.2 39.13 1.00 ×10−4613.89 3.00 ×10−3855.05 2.00 ×10−21160.54 2.00 ×10−2875.15 6.00 ×10−2839.19 6.00 ×10−2 0.1, 0.3 1485.63 3.00 ×10−21131.15 2.00 ×10−1905.77 2.00 ×10−11289.27 3.00 ×10−11070.57 7.00 ×10−11043.63 6.00 ×10−1