Statistics of partial minima

E Ben-Naim¹, M B Hastings¹ and D Izraelevitz²

¹ Theoretical Division and Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
² Decision Applications Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA

E-mail: ebn@lanl.gov, hastings@lanl.gov and izraelevitz@lanl.gov

Received 19 September 2007, in final form 10 October 2007
Published 6 November 2007

Abstract. We study pseudo-optimal solutions to multi-objective optimization problems by introducing partial minima defined as follows. Point x k-dominates x ' when at least k of the coordinates of x are smaller than the corresponding coordinates of x '. A point not k-dominated by any other point in the set is a k-minimum or a partial minimum, generalizing the global minimum. We study statistical properties of partial minima for a set of N points independently distributed inside the d-dimensional unit hypercube using exact probabilistic methods and heuristic scaling techniques. The average number of partial minima, A, decays algebraically with the total number of points, A ~ N–(d–k)/k, when 1 ≤ k < d. Interestingly, there are k – 1 distinct scaling laws characterizing the largest coordinates: the distribution P(yj) of the jth largest coordinate, yj, decays algebraically, , with for 1 ≤ j ≤ k – 1. The average number of partial minima grows logarithmically, , when k  =  d. The full distribution of the number of minima is obtained in closed form in two dimensions. PACS numbers: 02.50.Cw, 05.40.–a, 89.20.Ff, 89.75.Da

1. Introduction

A host of decisions in computer science, economics, politics and everyday life involve multiple criteria or multiple objectives [1–4]. A pedestrian choosing a walking path considers the distance, the number of turns and the number of traffic lights. In business, takeover bids are decided on a multitude of complex conditions in addition to the total monetary offer. In elections, voters examine how candidates stand on multiple issues.

In multi-objective optimization, a solution that is optimal with respect to all criteria is rarely possible and instead, one faces a set of choices that are suboptimal on at least one criterion. Decisions require algorithms to weed out clearly inferior choices, sort through all the remaining imperfect choices and evaluate the relative trade-offs between costs.

Since in multi-objective optimization, a global optimum is unlikely, we are interested in identifying points that are close to optimum. In this paper, we propose a pseudo-optimality criterion, and derive the likelihood of finding pseudo-optimal solutions as a function of the number of choices.

By definition, a global optimum is superior in all cost dimensions. Intuitively, one may define a pseudo-optimum as superior to all alternatives along a large number of cost dimensions. For example, in a three-cost scenario, there may not be any choice that is optimal with respect to all three costs, but we may be able to find choices that are better than any alternative along two costs, see figure 1. In a voting scenario, no candidate may have the most attractive position on all issues to a given voter. In this case, a voter might naturally restrict her attention to a candidate, or candidates, who have the most attractive position on as many issues as possible.

Let us represent our choices as N points in d dimensions with coordinates x ≡ (x₁, x₂,  ..., x_d). Each coordinate x_i > 0 represents a distinct cost. By convention, small-x values are superior and are considered dominant. Partial minima, a formalization of the pseudo-optima concept discussed above, are defined as follows. A point x is said to be k-dominated by x ' when at least k of the coordinates of x are larger than the corresponding coordinates of x '. A point is said to be a partial minimum, or formally a k-minimum, when it is not k-dominated by any other point in the set as illustrated in figure 1. We stress that a partial minimum is not required to dominate all other points on the same d – k coordinates and may dominate different points along different coordinates. The parameter 1 ≤ k ≤ d quantifies the quality of the partial minimum: a smaller k value represents a more stringent condition. The two extremes are the global minimum, k  =  1, where every coordinate is a minimum of the point set, and the efficient set, k  =  d, that includes all points that are not obviously dominated by other points as shown later in figure 3. Partial minima are conditional multivariate extrema and their properties are amenable to analysis using a statistical physics perspective [5–9].

In this study, we obtain exact statistical properties of partial minima including the multivariate density and its asymptotic behavior as well as scaling properties such as the typical size and average number of minima. We present two major results. First, as a function of the set size N, the average number of minima decays algebraically when 1 ≤ k < d, and grows logarithmically when k  =  d. Second, there are k – 1 different scaling laws for the largest coordinates, each following a power-law distribution with k – 1 distinct exponents. The rest of the d  +  1 – k coordinates are characterized by distributions with sharp tails. We also discuss the relevance of these results to the multi-objective shortest path on graphs, a central problem in multi-objective optimization.

We consider the situation where there are no correlations between the coordinates. That is, each coordinate is independently drawn from some distribution. As discussed below, this situation is equivalent to a uniform distribution in the unit hypercube. Thus, we conveniently assume that x_i is uniformly distributed in [0:1] for all 1 ≤ i ≤ d.

2. Heuristic arguments

Elementary scaling laws for the typical size of a partial minimum and the average number of minima are derived heuristically. We assume that (i) the partial minimum is dominant on a fixed set of k coordinates and (ii) all its coordinates are equal, x_i  =  x, for all i. By the partial minimum definition, the corresponding k-dimensional hypercube contains only the partial minimum itself. The volume of this hypercube is x^k and the expected number of points inside this hypercube must be of order 1, Nx^k ~ 1. Consequently, the typical size x decays algebraically with N,

This characteristic scale decreases as the minimum condition becomes more stringent, that is, as k decreases.

The expected number of partial minima

follows from the expected number of points inside the d-dimensional hypercube with linear dimension x, Nx^d. Partial minima are asymptotically rare and the scale 1 decays indefinitely. Furthermore, with a small probability, there is only one minimum when N is large. The scaling estimate 2 coincides with the exact value A  =  N^–(d–1) for k  =  1, since any point is a global minimum with probability N^–d. For k  =  d, the minimum in any one coordinate is a partial minimum and thus there is at least one partial minimum. Indeed, the decay exponent in 2 vanishes. This special case is discussed separately.

3. The density of minima

The density P_d,k(x) of k minima located at x is obtained analytically through a formal generalization of the heuristic argument above. For example, in two dimensions, the density is

The factor N is the number of ways to choose the minimum, and the second factor guarantees that the rest of the points do not dominate the minimum at (x₁, x₂). These points must not fall inside an L-shaped region of area x₁  +  x₂ – x₁ x₂ or equivalently 1 – (1 – x₁)(1 – x₂) when k  =  1, or a rectangle of area x₁ x₂ when k  =  2, as illustrated in figure 2.

In general, the density of minima

reflects that the N – 1 points are excluded from a d-dimensional region of volume G_d,k(x). The excluded volume obeys the recursion

In our notation, the dimensional index of a function dictates the dimension of its vectorial argument so the vectors on the right-hand side of 4 have d – 1 components. We obtain the recursion relation 4 by separating the excluded region into two regions: one in which the dth coordinate is dominant and one in which it is not. Using the boundary conditions G_d,0  =  1 and G_d,k  =  0 when k > d, we recover G_1,1  =  x₁ and G_2,1  =  x₁  +  x₂ – x₁ x₂. Furthermore,

In general, G_d,d  =  ∏^d_i = 1 x_i and G_d,1  =  1 – ∏^d_i = 1(1 – x_i).

4. Scaling

In the limit N → ∞, the product term x₁ x₂ in P_2,1  =  N[1 – (x₁  +  x₂ – x₁ x₂)]^N–1 is negligible compared with the linear term x₁  +  x₂ and thus,

Generally, only the kth degree terms are asymptotically relevant and the leading behavior is

The auxiliary function F_d,k(x) contains terms, each a distinct product of degree k. For example,

The auxiliary function equals the sum, F_d,1  =  ∑^d_i = 1 x_i, and the product, F_d,d  =  ∏^d_i = 1 x_i, in the two extremes. The function F_d,k(x) is defined recursively

for 1 ≤ k ≤ d with the boundary condition F_0,k  =  δ_k,0. This recursion follows from 4 by dropping the higher-degree term x_dG_d–1,k(x).

The asymptotic behavior 5 can be recast in the scaling form

as N → ∞. The scaling variable is z  =  x N^1/k, in accord with 1, and the scaling function is

The average number of k minima equals the integral of the density, A_d,k  =  ∫dx P_d,k(x), where ∫dx ≡ ∏^d_i = 1∫¹₀dx_i.^Note1 When k < d, the asymptotic behavior of the average follows from the scaling form 7, , and is in agreement with 2. The proportionality constant a_d,k equals the integral of the scaling function, a_d,k  =  ∫dzΦ_d,k(z), although now, the integration range is unrestricted, ∫dz ≡ ∏^d_i = 1∫^∞₀ dz_i. The prefactor is trivial for perfect minima, a_d,1  =  1, and otherwise, it can be obtained analytically only in a few exceptional cases.

5. Extreme statistics

Since global minima are constrained along all cost coordinates, extremely large costs are exponentially rare, whenever such a global minimum exists. Because we have relaxed the minimality condition, this may not necessarily be the case for partial minima. In our voting example, a candidate who is attractive to a voter on a multitude of issues may be extremely unattractive on a particular one. How likely is such a scenario?

We begin our study of extremal statistics [10–12] by first considering the distribution of the largest coordinate in a partial minimum. Without loss of generality, we order the coordinates x₁ < x₂ < · · · < x_d–1 < x_d. Our focus is on the tail of the distribution of the variable x_d, corresponding to the regime x_d x_d–1. We also restrict our attention to the limit N → ∞. The distribution Q₁(x_d) of the largest coordinate x_d equals the integral of the multivariate distribution with respect to the rest of the coordinates,

The second line is obtained by substituting the leading asymptotic behavior 5 and the third line reflects that only the first term in 6 is relevant when x_d x_i for all i < d. Our last step is to multiply and divide the third line by x_d and then invoke the scaling law 2 for the average number of k – 1 minima in d – 1 dimensions. In essence, we utilize the fact that when one of the coordinates is very large, the partial minima criterion involves one less constraint in one less dimension^Note2. The power-law decay of the distribution 9 shows that there is a substantial likelihood that x_d is relatively large.

The distribution Q₂(x_d–1) of the second largest coordinate x_d–1 is obtained using the bivariate distribution ,

The distribution Q₂(x_d–1) equals the integral of the bivariate distribution with respect to the largest coordinate, . This integral is dominated by the divergence at the lower limit of integration, and consequently,

The power-law tail is now steeper.

A similar calculation applies to the distributions of the k – 1 largest elements. In general, the distribution Q_j(y_j) of the jth largest element, y_j, with the definition y_j ≡ x_d + 1–j, decays as a power law,

for 1 ≤ j ≤ k – 1. The decay exponent increases monotonically with the index j,

We can verify the decay law 2 using , where the lower limit of integration is set by the typical size scale 1. Interestingly, there are k – 1 distinct scaling behaviors for the k – 1 largest elements. Each of these extremal coordinates is distributed according to a power-law distribution that is characterized by a distinct exponent.

This multiscaling behavior affects the behavior of the moments y^m_j defined as follows, y^m_j  =  I_m/I₀, where . The integral I_m is dominated by the divergence at the lower cutoff when the order is small, m ≤ α_j, but otherwise, the integral I_m is finite. Consequently, the moments have the following scaling dependences on N:

Low-order moments exhibit ordinary scaling behavior as they are characterized by the typical size scale 1 that underlies the multivariate distribution function 8. As usual, there is a logarithmic correction at the crossover. High-order moments plateau at a fixed value that is independent of the index m, an indication that there is a significant probability that the extreme elements are of order 1. Interestingly, the average size of the different coordinates may follow different scaling laws. For example, there are two scaling laws, y₁ ~ N^–1/6 and y₂ ~ N^–1/3 when d  =  4 and k  =  3. Of course, the sum ∑^d_i = 1 x_i has the same extremal statistics as does x_d.

The crossover moment or equivalently the exponent α_j diverges as k → j. Therefore, the smallest d  +  1 – k coordinates exhibit the ordinary scaling behavior

for k ≤ j ≤ d, and all moments of the respective distribution functions must be finite. In these cases, the distribution functions Q_j have tails that are as sharp as or sharper than an exponential. In the aforementioned case d  =  4 and k  =  3, the third and the fourth largest coordinates exhibit the ordinary scaling, y₃ ~ y₄ ~ N^–1/3.

6. Efficient sets

The set of points that are not dominated on all coordinates by any other point are partial minima when k  =  d (figure 3). We refer to this set as the `efficient set'. The efficient set, also termed the efficient frontier or Pareto equilibria, plays a central role in multi-objective optimization and has been studied extensively in economics, computer science, operations research and game theory [13, 14]. Since there is no objective trade-off between costs, every point in the efficient set is potentially a solution to the multi-objective optimization problem. The study of the properties of efficient sets was the original motivation for our research.

In the special case k  =  d, the expected size of the efficient set, E_d(N) ≡ A_d,d(N), obeys the recursion

The point with the largest x_d coordinate certainly does not dominate any other point. Furthermore, this point is on the efficient set if and only if the rest of its d – 1 coordinates are not dominated by any other point. This event occurs with probability and hence the second term in the recursion. We note that the recursion 16 can also be obtained by performing the integration over x_d in E_d(N)  =  N∫dx[1 – x₁ x₂ · · · x_d]^N–1. This integration is analytically feasible only if k  =  1 or k  =  d.

The recursion relation 16 is subject to the boundary condition E₁(N)  =  1. In two dimensions,

or alternatively, E₂(N)  =  H(N), where is the harmonic number. The average size of the efficient set grows logarithmically, E₂(N)  =  ln N  +  γ  +  · · ·, where γ  =  0.57721 is Euler's constant. In three dimensions, we have , and asymptotically, . The large-N behavior is obtained in general by converting the difference equation 16 into a differential equation dE_d/dN  =  E_d–1/N. The expected size of the efficient set grows logarithmically,

This logarithmic growth reflects that the integral of the scaling function, ∫dzΦ_d,d(z), is divergent at the upper limit. A straightforward generalization of the calculation above shows that the distribution of the extremal coordinates has a logarithmic correction,

for 1 ≤ j ≤ d – 1. We can verify that the average number of points is consistent with the exact behavior as in 18. The crossover moment vanishes and the moments decay logarithmically,

where m > 0 and 1 ≤ j ≤ d – 1.

7. Two dimensions

For the case d  =  2, we obtain closed form expressions for the distribution function of partial minima. This permits us to establish central limit-type behaviors for the distribution of the size of the efficient set.

In two dimensions, the distribution function p_n(N) for the event that the efficient set includes n points, where 1 ≤ n ≤ N, satisfies the recursion^Note3

and is subject to the boundary condition P_n(0)  =  δ_n,0. On the square, there are two coordinates: x₁ and x₂. We can derive 21 by alluding to the same reasoning behind 16, i.e., the point with the largest x₂ coordinate will be on the efficient set if and only if its x₁ coordinate is minimal, an event that occurs with probability N^–1.

Recursion equations for the average E(N)  =  n and the variance V(N)  =  n² – n² with f(n) ≡ ∑^N_n = 1 f(n)P_n are obtained by summing 21. The average satisfies E(N)  =  E(N – 1)  +  N^–1 in accord with 16 and the variance satisfies V(N)  =  V(N – 1)  +  N^–1 – N^–2. Thus, the variance equals the difference between the first and the second harmonic numbers

where H⁽²⁾(N)  =  ∑^N_n = 1 n^–2. The variance and the average have identical leading asymptotic behaviors, .

With the transformation , the auxiliary function satisfies the recursion

with . This recursion defines the Stirling numbers [15], so . Therefore, the full probability distribution is expressed in closed form,

for 0 ≤ n ≤ N.

The general asymptotic behavior, derived in [16],

applies in the limit n → ∞, N → ∞ with the ratio n/ln N finite. For small n ln N, the distribution is Poissonian, P_n(N)  =  N^–1(ln N)^n–1/(n – 1)! and for large n, the distribution approaches a Gaussian centered at the average E(N) ≊ ln N with the variance V(N) ≊ ln N,

We note that the convex hull, a subset of the efficient set (see figure 3), is characterized by similar statistical properties including a limiting Gaussian distribution and logarithmic growths, albeit with different prefactors, of the average and the variance [17–19].

8. Multi-objective shortest path

The multi-objective shortest path on a graph is defined as follows. Consider a graph, possibly with multiple edges connecting pairs of nodes, with d different costs on each edge. Fix the source and the destination nodes, and then consider all paths from source to destination, assigning d total costs to each path computed as the sum of the d individual costs of the path's constituent edges. The multi-objective shortest path problem consists of finding the efficient set of paths. Generally, finding the efficient set is an NP-hard problem, although less demanding approximation schemes exist [20, 21]. Nevertheless, the computation time of the approximation scheme depends crucially on the size of the efficient set.

Suppose the edge costs are independent, random draws from a common distribution. We can consider two limiting topologies. First, for a graph of two nodes connected by N edges, the number of elements in the efficient set grows poly-logarithmically in the number of edges as shown in 18. Second, for a one-dimensional chain of nodes where each pair of neighboring nodes is connected by a pair of edges, the total path costs become correlated [20], even though the individual edge costs are not. We have conducted numerical studies that find that the size of the efficient set is highly sensitive to the distribution of edge costs. Assuming each edge has two costs (w₁, w₂), both chosen from some continuous distribution, the convex hull grows linearly in the length of the chain. Interestingly, we observed various behaviors for the size of the efficient set, ranging from linear in the length of the chain, to power-law behavior, characterized by exponents greater than unity, up to stretched exponential behavior.

Finally, we consider Erdös–Renyi random graphs [22–24]. Using the fact that the shortest path between two randomly chosen nodes grows logarithmically with the total number of nodes in the graph and the fact that paths that are close in length to the shortest path weakly overlap and hence their costs are weakly correlated, the results in this paper can be used to heuristically show [25] that the size of the efficient set of paths grows poly-logarithmically with the number of nodes as in 18. This number is much smaller than the number found for chains where the paths are correlated.

9. Conclusions

We proposed partial minima as a protocol for identifying pseudo-optimal solutions to multi-objective optimization problems. Partial minima are defined by a parameter k: a point in d dimensions that dominates all other points on at least d – k coordinates is a partial minimum. As this optimality criterion becomes more stringent, partial minima improve in quality but are less probable. In the extreme case k  =  d, the number of partial minima grows logarithmically with the total number of points.

Remarkably, there is a series of distinct power-law distributions that characterize the largest coordinates with a consequent multiscaling distribution of the moments, while the rest of the coordinates obey ordinary scaling. Viewed as quasi-optimal solutions to multi-objective optimization problems, partial minima involve a trade-off. When the optimality criterion is relaxed, these quasi-optima become more likely, but are more likely to incur at least one extremely large cost.

Our results hold as long as the set of points are not correlated, that is, as long as they are drawn from independent distributions. These distributions need not be identical. If the ith coordinate is drawn from the distribution f_i(x_i), the transformation and dx_i → f_i(x_i) dx_i maps to a uniform distribution in the unit hypercube. Correlations present an interesting challenge and we anticipate serious modifications to the scaling laws above. For instance, it is simple to show that the size of the efficient set grows as a power of the number of points, ~N^1/2, rather than a logarithm, when the points are uniformly distributed inside the unit circle. Incidentally, this growth is much faster than the N^1/3 for the corresponding number of points in the convex hull [17].

Another interesting issue is the crossover from the algebraic decay 2 to the logarithmic growth 18. The average number of partial minima decreases monotonically with N when k is small, but is a non-monotonic function of N when k is large. For example, when d  =  4 and k  =  3, the average A_d,k peaks at N  =  16. It will be interesting to elucidate how the height and the location of this peak scales with N.

Acknowledgments

We thank Gunes Ercal-Ozkaya for suggesting random graphs and Paul Krapivsky for useful discussions on convex hulls. We acknowledge financial support from DOE grant DE-AC52-06NA25396.

[1] 

Steuer R E 1986 Multiple Criteria Optimization: Theory, Computations, and Application (New York: Wiley)  

[2] 

Sawaragi Y, Nakayama H and Tanino T 1985 Theory of Multiobjective Optimization (Orlando: Academic)  

[3] 

Fudenberg D and Tirole J 1983 Game Theory (Cambridge, MA: MIT Press)  

[4] 

Osborne M J and Rubenstein A 1994 A Course in Game Theory (Cambridge, MA: MIT Press)  

[5] 

Mézard M, Parisi G and Virasoro M 1987 Spin Glass Theory and Beyond (Singapore: World Scientific)  

[6] 

Derrida B and Lebowitz J L 1998 Phys. Rev. Lett. 80 209 
CrossRef

[7] 

Ben-Naim E, Krapivsky P L and Majumdar S N 2001 Phys. Rev. E 64 035101 
CrossRef

[8] 

Korniss G, Toroczkai Z, Novotny M A and Rikvold P A 2000 Phys. Rev. Lett. 84 1351 
CrossRefPubMed

[9] 

Krapivsky P L and Redner S 2002 Phys. Rev. Lett. 89 258703 
CrossRefPubMed

[10] 

Gumbel E J 1958 Statistics of Extremes (New York: Columbia University Press)  

[11] 

Galambos J 1987 The Asymptotic Theory of Extreme Order Statistics (Malabar, FL: Krieger)  

[12] 

Ellis R E 1985 Entropy, Large Deviations and Statistical Mechanics (New York: Springer)  

[13] 

Ottmann T, Soisalon-Soininen E and Wood D 1984 Inform. Sci. 33 157 
CrossRef

[14] 

Fink E and Wood D 1998 J. Geom. 62 99 
CrossRef

[15] 

Graham R L, Knuth D E and Patashnik O 1989 Concrete Mathematics A Foundation for Computer Science (Reading, MA: Addison-Wesley)  

[16] 

Ben-Naim E and Krapivsky P L 2004 J. Phys. A: Math. Gen. 37 5949 
IOPscience

[17] 

Rényi A and Sulanke R 1963 Z. Wahrsch. Verw. Geb. 2 75–84 
CrossRef

[18] 

Efron B 1965 Biometrika 52 331 

[19] 

Vu V 2006 Adv. Math. 207 221 
CrossRef

[20] 

Warburton A 1987 Oper. Res. 35 70 
CrossRef

[21] 

Tsaggouris G and Zaroliagis C 2006 Lect. Notes Comput. Sci. 4288 389 
CrossRef

[22] 

Bollobás B 1985 Random Graphs (London: Academic)  

[23] 

Janson S, Łuczak T and Rucinski A 2000 Random Graphs (New York: Wiley)  

[24] 

Janson S 1999 Comb. Probab. Comput. 8 347 
CrossRef

[25] 

Ben-Naim E, Hastings M B and Izraelevitz D unpublished  

Notes

Note1

 For small N, these integrals can be calculated manually or through recursion formulae. For example, when N  =  2.

Note2

 Formally, when x_d  =  1 then G_d,k(x)  =  G_d–1,k–1(x).

Note3

 Generally, p₁(N)  =  A_d,1  =  N^–(d–1).