Approximation Techniques for the Set of Efficient Points

by Jörg Fliege

Professorial dissertation for obtaining the venia legendi at the Fachbereich Mathematik, Universität Dortmund, Germany

Submitted on January 30, 2001.

Acknowledgments

Chapter 4 of this work is concerned with goal programming problems, a type of problem often appearing in subproblems of multicriteria optimization, especially if the solution strategies described in Chapter 2 and 3 are used. Most of the contents of Chapter 4 have been developed in conjunction with Emilio Carrizosa from the Universidad de Sevilla, Seville, Spain, and have been published as [12] and submitted to Mathematical Programming. Chapter 5 pursues an approach which might be interpreted as a ”dual” strategy to the one followed in Chapter 4. Most of the contents of Chapter 5 have been developed in conjunction with Benar Svaiter from the Instituto Matemática Pura e Aplicada (IMPA), Rio de Janeiro, Brazil, and have been published as [38]. I would hereby like to thank Emilio and Benar for their support, their hospitality, their friendship, and the fun we had while doing research and other things at the various places where we met.

Contents

1 Introduction 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Possible Lines of Attack . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2.1 General Value Functions . . . . . . . . . . . . . . . . . . . . . . 2 1.2.2 Complete Discretization . . . . . . . . . . . . . . . . . . . . . . 3 1.2.3 Total Linearization . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.4 Parameter-Free Techniques . . . . . . . . . . . . . . . . . . . . . 4

1.3 An Approximation Approach . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3.1 Some Thoughts on Parallelization Approaches . . . . . . . . . . 5

2 Multicriteria Problems 7

2.1 The Basic Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.1 Minimal Points . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.2 Several Objective Functions . . . . . . . . . . . . . . . . . . . . 10

2.2 Scalarizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.1 Monotone Functions . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.2 Proper Efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2.3 Approximating the Set of Efficient Points . . . . . . . . . . . . . 17 2.2.4 The Nonconvex Case . . . . . . . . . . . . . . . . . . . . . . . . 23

3 Interior-Point Methods for Convex Multicriteria Problems 27

3.1 The Standard Primal Path-Following Interior-Point Algorithm . . . . . . 28 3.1.1 The First Stage . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.1.2 The Second Stage . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.2 Solving Two Single-Criteria Programs with Linear Objective Functions . 35 3.3 Solving Two Single-Criteria Programs with Quadratic Objective Functions 42 3.4 Solving Two Single-Criteria Programs with General Objective Functions . 46 3.5 Solving Many Convex Programs ”Simultaneously” . . . . . . . . . . . . 46

3.5.1 The General Idea . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.5.2 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.5.3 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

3.6 More than Two Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.6.1 Voronoi-Tesselations . . . . . . . . . . . . . . . . . . . . . . . . 53 3.6.2 Decomposition-based Strategies . . . . . . . . . . . . . . . . . . 54 3.6.3 Adaptive Grids . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4 Generalized Convex Goal Programming 57

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.1.1 Goal Programming . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.1.2 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

4.2 Asymptotically Conical Sets and their Properties . . . . . . . . . . . . . 62 4.3 The Problem Addressed and its Dual . . . . . . . . . . . . . . . . . . . . 63 4.4 Existence of Primal and Dual Solutions . . . . . . . . . . . . . . . . . . 70 4.5 Solving Goal Programs Efficiently . . . . . . . . . . . . . . . . . . . . . 76

4.5.1 Particular Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

5 Nonparametric Descent Methods for Nonconvex Multicriteria Problems 87 5.1 The Problem Considered . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.2 Computing the Search Direction . . . . . . . . . . . . . . . . . . . . . . 89

5.2.1 Other Possibilities for the Search Direction . . . . . . . . . . . . 91 5.3 Computing the Step Length . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.4 The Complete Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.5 Convergence of the Steepest Descent Method . . . . . . . . . . . . . . . 94 5.6 Towards an implementable Method . . . . . . . . . . . . . . . . . . . . . 96 5.7 The Constrained Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.8 A Quadratic Approximation Strategy . . . . . . . . . . . . . . . . . . . . 101

Bibliography 104

Chapter 1

Introduction

1.1 Motivation

In multicriteria optimization, several objective functions have to be minimized simultaneously. Problems of this type are also labeled as belonging to multiobjective optimization, multiple objective optimization, multiple criteria mathematical programming, or, since the function to be minimized has a vector space of dimension strictly greater than 1 as it’s image space, vector optimization. While sometimes the different names mentioned above denote different problem types, we will call all problems with more than one objective function to be belonging to the class of multicriteria optimization problems. Applications for these problems can be found in engineering design [26] (especially truss optimization [17] and antenna design [58]) location science [11], statistics [13], management science [27] (especially portfolio analysis [93], production & project planning [46, 1], and scheduling [5, 73]), environmental analysis [62, 32, 34], cancer treatment planning [53], etc. While methods for solving multicriteria problems have traditionally used ideas from single-criteria optimization (it can be argued that single-criteria optimization is simpler than then general multicriteria case — see below), ideas from the multicriteria case have recently been used successfully for single-criteria problems [29, 30].

Usually, no single point will minimize all of the several objective functions given at once. This is due to the fact that the image space IR^{n} with n > 1 can not be totally ordered. Therefore, the concept of optimality has to be replaced by a weaker concept, usually called efficiency or Pareto-optimality. This concept arises from a prespecified order relation on IR^{n} which is not total, but satisfying some suitable assumptions. Typically, a point is called Pareto-optimal or efficient, if there does not exist a different point with the same or smaller objective function values, such that there is a decrease in at least one objective function value. Note, however, that we will make use of a slightly more general definition, to be introduced in Chapter 2.

In classical single-criterion (scalar) optimization, the situation is comparatively simple: each minimum of the objective function is indistinguishable from each other one in terms of the image space, because the function values of all minima are the same. Since a good model should include all criteria to be minimized which are quantifiable, different solutions of a given single-criterion optimization problem can only be distinguished from each other in terms of non-quantitative criteria. This is in stark contrast to the situation in multicriteria optimization, where different minimal points can be distinguished from each other usually quite well, since they have different objective function values for the objective functions at hand. Moreover, different minimal points will be incomparable to each other, i. e. no single minimal point can be used as an approximation of all the other minimal points. Both effects occur no matter what kind of ordering relation is used in the image space, because these effects are induced by the abovementioned nonexistence of a total order in IR^{n} (n > 1). The next section lists some possible approaches leading out of this dilemma, while Section 1.3 describes which strategy is pursued in the main part of this work.

1.2 Possible Lines of Attack

In this section, some possible solution strategies for continuous multicriteria problems are discussed. Discrete (or combinatoric) multicriteria problems fall outside the scope of this work. For a taxonomy similar to the one presented here, the reader is referred to Benson and Sayin [7].

1.2.1 General Value Functions

One of the main solution strategies for multicriteria optimization problems nowadays is the scalarization approach, first described by Geoffrion [43]. Here, one or several parameterized single-objective (i. e. classical) optimization problems are solved instead of one multicriteria problem. Parameters occur because there are several functions mapping IR^{n}, the image space of the given problem, into IR, the image space of a classical optimization problem. The disadvantage to this approach is that the choice of the parameters is not known in advance, leaving the modeler and the decision-maker with the burden of choosing them. Moreover, this method computes only so-called proper Pareto-optimal points, although it has to be said that under suitable assumptions the set of proper efficient points can not be distinguished from the set of efficient points in floating-point arithmetic. Nevertheless, using the right family of single-criterion optimization problems, one is able (at least in principle) to compute all efficient points.

It is a modified version of this strategy that is used as a tool in Chapter 3 of this work, while the structure of an important subclass of single-criterion optimization problems occurring in value function approaches is discussed in Chapter 4.

1.2.2 Complete Discretization

One might consider a complete discretization of the space in use. Typically, if the value space is the IR^{n}, one uses a Cartesian grid in the IR^{n} and considers in the next step the grid points only to define a discrete multicriteria problem. Then, standard methods from discrete multicriteria analysis might be used. While this approach has some merits for small scale problems — primarily due to the fact that solution methods for discrete multicriteria problems are now rather well developed [82, 80, 81] — it will often be the case that the size of the resulting discretized problem prohibits any further solution approach. As such, a complete discretization approach remains questionable. However, a form of adaptive discretization might able to alleviate the inherent inefficiency of such a method.

1.2.3 Total Linearization

Linearizing all functions involved results in a linear multicriteria optimization problem, for which Simplex-based solver strategies exist (see, e. g., [98, 27, 87, 88]). In singlecriteria optimization, a complete linearization is seldom used, since the resulting linear problem is rather large. From a theoretical point of view, Simplex-based methods of todays standards still lack the efficiency of other methods, especially interior-point approaches.

Moreover, the number of efficient facets and extreme points of the feasible polyhedron usually grows exponentially with the dimension of the problem. But the Simplex-based solvers mentioned above try to enumerate all efficient extreme points resp. facets. It is therefore highly questionable if such a strategy is feasible for problems of larger size. This theoretical argument is augmented by computational tests [7]: with 4 objective functions and 40 unknowns, a polyhedron described by 50 inequalities can have more than 200 000 efficient vertices.

On the other hand, it might not be necessary to linearize the objective functions everywhere. Local linearization is probably the most important technique for single-criteria optimization, and a generalization of this approach is therefore pursued in Chapter 5.

1.2.4 Parameter-Free Techniques

Recall that in value function strategies a function mapping the image space of the given multicriteria problem into the real numbers is used to transform the multicriteria problem into a classical optimization problem. Of course, such a function contains several parameters, making it a priori not clear which ones have to be chosen. Parameter-free techniques do not use such a parameterization.

Some so-called parameter-free multicriteria optimization techniques use an ordering of the different criteria, i. e. an ordering of importance of the components of the objective function vector. In this case, the ordering has to be specified, which is in itself a set of parameters. In other parameter-free techniques, the optimization process is augmented by an interactive procedure (see, e. g. [87, 65]), adding an additional burden to the task of the decision-maker. This particular approach is infeasible for real-time optimization and problematic in large-scale optimization. Moreover, usually only a small amount of interaction ever takes place in practice, although most methods require a rather large amount theoretically [87].

Nevertheless, when the number of criteria is large and the approach described in Section 1.3 might become too time-consuming, parameter-free techniques can become useful. One of these techniques, not needing any kind of interactivity, is described in Chapter 5.

1.3 An Approximation Approach

The strategy proposed in Chapter 2 and 3 for solving multicriteria optimization problems can be described as follows: the discussion in Section 1.1 has shown that knowledge about the whole set of efficient points is necessary to finish the task of solving a multicriteria optimization problem. The actual computation of this set is out of question for realworld problems. What is therefore needed is an approximation to this set, preferably with a well defined approximation quality. The simplest approximation of a given subset of the IR^{n} is a discrete approximation, consisting of a finite set of points. Theoretical results described in Chapter 2 show that under suitable assumptions — mainly convexity — there exists a local Lipshitz-continuous function mapping arguments from a bounded set onto the efficient points. Therefore, it makes sense to evaluate this function at some suitable arguments to build an approximation for it’s graph and thereby for the solution set of the underlying convex multicriteria program. However, each evaluation of this function is rather costly, since such an evaluation involves solving a scalar optimization problem. For reasons of computational efficiency, it makes therefore sense to reuse as much information of previous evaluations of the function under consideration (i. e. of previous scalar optimization runs) as possible. A corresponding warm-start technique of theoretical efficiency for convex multicriteria programs is developed in Chapter 3.

After a discrete approximation to the set of efficient points has been computed, there is still a need to choose one of these points to finish the optimization process. In case the number of discrete points is small (and the approximation therefore comparatively bad), this can be done by hand. Otherwise, methods from discrete multicriteria optimization, already well developed [82, 80, 81], can be used. Alternatively, one might construct a continuous approximation of the set of efficient points by way of the constructed discrete approximation and execute an additional optimization step with some additional objective function. All these possibilities effectively boil down to solving an optimization problem in the set of efficient points, i. e. with the set of efficient points as the set of feasible points. In any case, the discussion of this final step falls outside the scope of this work.

The necessity for a strategy like the one outlined above was also observed by Benson and Sayin [7] and by Das and Dennis [18, 17]. However, in none of these works the methods under consideration are subject to a worst-case analysis with respect to their computational complexity, and while it is claimed that the method developed in [18] generates a discrete approximation to the set of efficient points, a counterexample by the authors themselves [18, Figure 3] shows that this is not the case.

Note that the theoretical efficiency mentioned above has to be interpreted more carefully than in standard scalar optimization. If a multicriteria problem with n objective functions is given and we are in search of a δ-covering of the set of efficient points (i. e. a set of points such that each efficient point lies in a ball of radius δ around one of the given points), the output size of any algorithm for this problem is at least O(1/δ^{n−1}), since the set of efficient points is usually an n − 1-dimensional manifold. The situation is therefore similar to n-dimensional integration or global optimization of (unstructured) functions of n unknowns: there exist exponential lower bounds for all numerical algorithms. But, of course, this does not mean that algorithms of arbitrarily bad run-time behavior should be used. Quite the contrary: the use of efficient algorithms is more important than ever. Moreover, the number of criteria n is usually quite small. In contrast to this, it is the number of unknowns (decision variables) that is increasing strongly when a model is refined to a more realistic one. The class of methods proposed in Chapter 3, based on the results described in Chapter 2, has exponential run-time behavior in the number of criteria (simply because the output size grows exponential), but polynomial time behavior in the number of unknowns.

1.3.1 Some Thoughts on Parallelization Approaches

Many of the strategies outlined above can, in general, be implemented on parallel computers, and the same can be said about the approximation technique outlined above. Especially value function strategies, in which each processor of a parallel machine is able to solve a scalar optimization problem with a different parameter set, are particularly suited for this approach, at least in principle. However, a standard parallelization strategy will simply distribute the task of solving these problems among several processors, thereby guaranteeing a speedup. But the discussion in the following chapters will show that the scalar problems which have to be considered are constructed by perturbing a given problem, and information gained while solving one of them can and should be used for solving the others. Without proper communication between the different tasks, there will be a large amount of superfluous work done on a parallel machine, both in terms of unnecessary communication and in terms of unnecessary optimization steps. Therefore, the approach chosen in this work can be used to make parallelization efforts even more flexible and efficient.

Multicriteria Problems

2.1 The Basic Problem

Let IR denote the set of real numbers and IR_{+} = [0, +∞[ the set of nonnegative real numbers. The nonnegative orthant in the IR^{n} will be denoted by IR_{+}^{n} := (IR_{+})^{n}. In what follows, for given nonempty subsets S_{1}, S_{2} of IR^{n}, we denote by S_{1} + S_{2} the algebraic sum of S_{1} and S_{2},

S_{1} + S_{2} := {s ∈ IR^{n} | s = s_{1} + s_{2} for some s_{1} ∈ S_{1}, s_{2} ∈ S_{2}}

When S_{1} is a singleton, S_{1} = {s_{1}}, we write s_{1}+S_{2} to represent {s_{1}}+S_{2}. For arbitrary sets S ⊆ IR^{n} we will denote by cl(S), int(S), bd(S), relint(S), and aff(S) the closure, the interior, the boundary, the relative interior, and the affine hull of S, respectively. Moreover, for a convex set S and x ∈ S, the normal cone to S at x will be denoted by

N_{S}(x) := {ω ∈ IR^{n} | ∀ y ∈ S : ⟨ ω, y − x ⟩ ≤ 0}.

2.1.1 Minimal Points

In order to speak about minima of an arbitrary set M ⊆ IR^{n}, one has to define an order relation on IR^{n}. Moreover, given two vectors in IR^{n}, it is important for reasons of computational efficiency to be able to check efficiently if these two vectors are ordered according to the given relation. The standard approach consists of the following definition.

Definition 2.1.1 Let K ⊂ IR^{n} be a set. The binary relation ≤_{K} on IR^{n} is defined by

y ≤_{K} x if and only if x − y ∈ K.

A vector x is said to be dominated by y, if y ≤_{K} x holds.

The definition of domination above was made with respect to minimization problems in mind.

The class of binary relations defined includes total orders like the standard lexicographic order. However, a total order relation on IR^{n} (n > 1) has severe disadvantages when optimization problems have to be solved [28, Chapter 7,8] (see also [60, Abschnitt 5.2.2] for the case of the lexicographic order), and it can be argued that the use of a total order amounts to imposing unrealistic assumptions in a multicriteria framework [60, Abschnitt 2.4].

The next theorem is well known, cmp., e. g., [45].

Theorem 2.1.1 Let K ⊂ IR^{n} be a set and let ≤_{K} be the binary relation defined by K. Then, the following statements hold:

1. If 0 ∈ K then ≤_{K} is reflexive.

2. If K + K ⊆ K then ≤_{K} is transitive.

3. If K is a cone containing no lines then ≤_{K} is anti-symmetric.

4. If K ∪ (−K) = IR^{n} then ≤_{K} is total.

5. The set K is closed if and only if the relation ≤_{K} is ”continuous at 0” in the following sense. For all x ∈ IR^{n} and all sequences (x^{(i)})_{i∈IN} in IR^{n} with lim_{i→+∞} x^{(i)} = x and 0 ≤_{K} x^{(i)} for all i ∈ IN it follows that 0 ≤_{K} x holds.

Proof: All claims follow immediately from their assertions. □

Note that K + K ⊆ K holds if K is a convex cone. Moreover, if K is a cone containing no lines and the order ≤_{K} is total then K is not closed, a rather problematic situation with respect to numerical algorithms.

According to the theorem above, it makes sense to choose a closed convex cone K with 0 ∈ K which contains no lines to define the partial order ≤_{K}. Moreover, in our context the space IR^{n} will be the image space of functions to be minimized. As a consequence, it is important for numerical reasons to have scale-invariance of the induced order. This means that

if x ≤_{K} y and λ > 0 then λx ≤_{K} λy,

a property which holds if and only if the set K is a cone.

As a consequence of the discussion above, one might choose as K ⊂ IR^{n} a nonempty closed convex pointed cone. Then, ≤_{K} is reflexive, transitive, anti-symmetric, scaleinvariant and continuous in the sense introduced above. If a relation ≤ on IR^{n} is given and there exists a cone K such that ≤ = ≤_{K} holds, then K is called the ordering cone of ≤. The most important cone used is, of course, the non-negative orthant IR_{+}^{n}. However, we might be interested in using a slightly different cone for numerical reasons [29], [16, Section 15.5], and in semidefinite optimization problems, the cone of positive semidefinite matrices occurs naturally as an ordering cone [91, 56]. Moreover, complications arise as soon as the underlying vector space has infinite dimension. While the discussion above as well the last theorem hold in infinite-dimensional spaces, too, care has to be taken whenever infinite-dimensional spaces are discretized and replaced by finite-dimensional ones.

Example 2.1.1 Consider the space C^{2}([ 0, 1], IR) of twice continuously differentiable functions f : [ 0, 1] ⟶ IR over the unit interval. Then the set of all nonnegative functions f ∈ C^{2}([ 0, 1], IR) is a cone in this space. Let [ 0, 1] be discretized by an arbitrary grid of n ∈ IN subintervals. If a polynomial of degree 0 is used on each subinterval to approximate a given function in one of the usual ways (i. e. by interpolating at the left boundary of each interval, or at the right boundary, or at the midpoint), the cone under consideration is approximated by IR_{+}^{n}. The same holds with respect to IR_{+}^{n+1} if polynomials of degree 1 are used to construct a piecewise linear approximation to a given function, interpolating at the boundary of each interval. Suppose now that polynomials of degree 2 (or higher) are used to define a piecewise polynomial approximation to a given nonnegative smooth function f. Then, the corresponding parameters can be identified with a vector p ∈ IR^{3n}, but usually we can neither guarantee the nonnegativity of all p_{i} nor that an arbitrary vector p ∈ IR_{+}^{3n} defines a nonnegative piecewise polynomial function. Adding smoothness conditions at the endpoints of each interval reduces the dimension of the parameter space, but the infinite-dimensional cone is still not mapped onto the positive orthant in the corresponding finite-dimensional space.

If K is a cone, the dual cone

K^{∗} := {z ∈ IR^{n} | ∀ x ∈ K : ⟨ z, x ⟩ ≥ 0}

to K is denoted by K^{∗}. The quasi-interior of K^{∗} is the set

K^{+} := {z ∈ IR^{n} | ∀ x ∈ K \ {0} : ⟨ z, x ⟩ > 0}.

Next, the notion of minimality or optimality with respect to ≤_{K} has to be defined.

Definition 2.1.2 (Efficiency) Let K ⊂ IR^{n} be a set and ≤_{K} be the binary relation defined by K. Moreover, let M ⊆ IR^{n}. The point x ∈ M is efficient with respect to ≤_{K} if there does not exist a point y ∈ M with

y ≤_{K} x and y ≠ x.

The set of efficient points of M with respect to ≤_{K} will be denoted by E(M, K).

The set E(M, int(K) ∪ {0}) is called the set of weakly efficient points with respect to ≤_{K}.

It is this notion of optimality that will be used further on. (Weakly) efficient points are also called (weakly) nondominated points.

Clearly, x ∈ M is efficient if and only if there does not exist a point y ∈ M with x − y ∈ K \ {0}, i. e. M ∩ (x − K \ {0}) = ∅. Note that E(M, K) = E(M + K, K) always holds and that E(M, K) ⊆ bd(M) holds for K ≠ {0}.

Since E(M, K_{1}) ⊆ E(M, K_{2}) for K_{2} ⊆ K_{1} it might make sense to consider cones enclosing the standard cone used in most applications, IR_{+}^{n}. In this way, a sensitivity analysis of the set of efficient points is possible. Moreover, as Example 2.1.1 has shown, not all discretizations of cones in infinite-dimensional spaces lead to the positive orthant of the IR^{n}.

It is of prime importance to note that in applications the construction of one minimal (efficient) element of the set M is usually not sufficient (the only exception to this rule is the case of a total order ≤_{K}, a case already alluded to and dismissed on page 8). The set E(M, K) is for all practical purposes too large to be represented by only one of it’s elements. Moreover, in stark contrast to the standard scalar optimization case, different efficient points are well distinguished from each other by their differing coordinates. On the other hand, an explicit construction of the whole set of efficient points is usually out of question. As a consequence, one has to approximate E(M, K) in a certain way. This is the subject of Section 2.2.3.

2.1.2 Several Objective Functions

Let f : IR^{m} ⟶ IR^{n} be a function which will serve as the objective function to be minimized in the sense of the last section. Defining M := F (G), one is not only in search for the set of efficient points E(M, K), but also for it’s inverse image f^{−1}(E(M, K)). This inverse image is also called the set of Pareto-optimal points. In the literature, the latter set is sometimes identified with the set of efficient points. Likewise, f^{−1}(E(M, int(K) ∪ {0})) is the set of weakly Pareto-optimal points. In theory, the situation is not much more complicated than above, because if an ordering cone K ⊆ IR^{n} is given, one might always search for E(graph(f), IR^{m} × K). Note, however, that any antisymmetry one had in the relation ≤_{K} gets lost. A different possibility is discussed in Section 2.2.2.1.

2.2 Scalarizations

2.2.1 Monotone Functions

The concept of monotonicity of a function mapping real numbers on real numbers can easily be extended to functions defined on a vector space, provided that two ordering relations are given.

Definition 2.2.1 (Monotonicity) Let K ⊆ IR^{n}, L ⊆ IR^{ℓ} be two cones, S ⊆ IR^{n} a set, and f : S ⟶ IR^{ℓ} a function.

The function f is called (K, L)-monotonically increasing on S, if

x ≤_{K} y ⟹ f(x) ≤_{L} f(y)

holds for all x, y ∈ S.

The function f is called strictly (K, L)-monotonically increasing on S, if

x ≤_{K} y, x ≠ y ⟺ f(x) ≤_{L} f(y), f(x) ≠ f(y)

holds for all x, y ∈ S.

As an abbreviation, if ℓ = 1 and L = IR_{+} = [0, +∞[, a (strictly) (K, L)-mono-tonically increasing function f is called (strictly) K -monotonically increasing on S or simply (strictly) K -monotone on S. If, additionally, K is fixed, one might use an even less wordier phrasing and speak of monotone and strictly monotone functions.

Functions monotone with respect to an arbitrary binary relation ⪯ are also called consistent with respect to ⪯, see [92, Chapter 1], or order-preserving, see [28, Chapter 7] for an overview. These functions play an important role in multicriteria optimization, as the next theorems show.

Theorem 2.2.1 Let K ⊆ IR^{n}, L ⊆ IR^{ℓ} be two cones and M ⊆ IR^{n} be a set. Let f : M ⟶ IR^{ℓ} be a strictly (K, L)-monotone function on M. Then x ∈ E(M, K) if and only if f(x) ∈ E(f(M), L).

Proof: See [91, Kapitel 2]. □

Consequently, mapping a set M whose efficient points are searched for by a strictly monotone function does not destroy the efficiency property of the points in E(M, K). Indeed, the problem at hand might be simplified significantly by choosing a proper mapping f. This is shown by the next two theorems.

Theorem 2.2.2 (Face Destruction) Let K ⊆ IR^{n} be a cone with IR_{+}^{n} ⊆ K and let M ⊆ IR^{n} be a nonempty closed convex set. Let g_{i} : IR ⟶ IR (i = 1, …, n) be strictly convex and strictly monotone (in the usual sense) functions and define g : IR^{n} ⟶ IR^{n} by g := (g_{1}, …, g_{n})^{⊤}. Then, each face of g(M) has at most one point with E(g(M), K) in common.

Proof: Due to IR_{+}^{n} ⊆ K, we have E(M, K) ⊆ E(M, IR_{+}^{n}) as well as E(g(M), K) ⊆ E(g(M), IR_{+}^{n}). The function g is bijective and strictly (IR_{+}^{n}, IR_{+}^{n})-monotone, and therefore g(E(M, IR_{+}^{n})) = E(g(M), IR_{+}^{n}) follows.

Let F be a face of g(M) which meets E(g(M), K). Suppose that there exist two points x, y ∈ M with g(x), g(y) ∈ F ∩E(g(M), K) ⊆ bd(g(M)) but g(x) ≠ g(y). Due to the strict monotonicity of all g_{i}, we have x ≠ y. Moreover, [g(x), g(y)] ⊆ F holds. Now let λ ∈]0, 1[. Then λg(x) + (1 − λ)g(y) ∈ F ⊆ E(g(M), IR_{+}^{n}), and Theorem 2.2.1 shows that z := g^{−1}(λg(x) + (1 − λ)g(y)) ∈ E(M, IR_{+}^{n}). But g_{i}^{−1} is strictly concave, and therefore z_{i} > λx_{i} + (1 − λ)y_{i} for i = 1, …, n. This is a contradiction. □

Neither monotonicity nor strict convexity of the function g_{i} is needed on the whole IR. Instead, it suffices to have this property on the projection of the ith coordinate of elements of M.

According to the preceeding two theorems, it is possible to destroy any polyhedrality that E(M, K) might have by applying a strictly convex transformation. These are rather good news, since it will be shown later that a polyhedral structure is a major problem when approximating the set of efficient points by a discrete subset.

On the other hand, it is clear that under the conditions of the theorem above, the set g^{−1}(E(g(M), K)) can be a strict subset of the set E(M, K). This is so because in the case K ⊃ IR_{+}^{n} the function g is not necessarily strictly (K, K)-monotone, but only strictly (IR_{+}^{n}, IR_{+}^{n})-monotone. Therefore, Theorem 2.2.1 can not be used to conclude equality between the two sets in question. Of course, one might first compute an approximation of the set E(g(M), IR_{+}^{n}), transform this by g^{−1} to an approximation of the set g^{−1}(E(g(M), IR_{+}^{n})), and then check if the points computed lie in E(M, K). This, however, might mean additional work, and it is therefore important to gain knowledge about those points in E(M, K) missing from the set g^{−1}(E(g(M), K)). To this end, suppose that K ⊃ IR_{+}^{n} but K ≠ IR_{+}^{n} and that there exist an x ∈ E(M, K) with g(x) ∉ E(g(M), K). Because of g(x) ∈ E(M, IR_{+}^{n}) there exists a ξ ∈ K \ IR_{+}^{n} with g(x)−ξ ∈ g(M), i. e. g(x)−ξ = g(y) for some y ∈ M, y ≠ x. But x ∈ E(M, K), and therefore y − x ∉ −K. Suppose now in addition that g is continuously differentiable. Then, for small ξ one may write g(x) − ξ = g(y) ≈ g(x) + g′(x)(y − x). Therefore, ξ ≈ g′(x)(x − y), which means that a point not in K is mapped by the linear transformation g′(x) into K \ IR_{+}^{n}. Note that g′(x) is a diagonal matrix, so it represents just a scaling of variables. Accordingly, g(x) is not efficient in g(M) with respect to K, but approximately with respect to the cone

K̂(x, g) := K \ (g′(x) · (K \ IR_{+}^{n})) .

The following heuristic might therefore be of interest. Defining the cone

K̂(g) := ⋂_{x∈M} K̂(x, g),

one sees that E(g(M), K̂(g)) ⊆ g(E(M, K)), and the error induced by replacing the set E(g(M), K) by E(g(M), K̂(g)) can be controlled by choosing suitable functions g. More precisely, the closer g is to the identity mapping (in the sense of an arbitrary norm on the set of convex functions on M), the smaller will be the residual set E(g(M), K) \ E(g(M), K̂(g)).

An interpretation which might be more useful for practitioners can be obtained along the following way. In the setting above, ξ ∈ K \ IR_{+}^{n} holds, which means that there exists at least one index i ∈ {1, …, n} such that ξ_{i} < 0, but possibly ξ_{j} ≥ 0 for all j ≠ i. On the other hand, z := x−y ∉ K, so for ξ small we have that z_{i} < 0 holds. Then, g′_{i}(x) has to be small compared to g′_{j}(x), j ≠ i. If each function g_{i} (i = 1, …, n) is interpreted as a value function mapping a variable with dimension cost to another variable, also with the dimension cost, one sees that the marginal costs at the point x with respect to criterion i are small compared with the marginal costs of the criteria j ≠ i.

The results presented up to now have shown that it is possible to replace a multicriteria optimization problem with an equivalent problem of possibly simpler structure. However, all these equivalent problems are multicriteria problems, too. The next theorem shows a well-known connection between multicriteria and single-criteria problems which will be further exploited in what follows.

Theorem 2.2.3 Let K ⊆ IR^{n} be a cone with 0 ∈ K and {0} ≠ K ≠ IR^{n}. Let v : M ⟶ IR be a K -monotone increasing function and let x ∈ M be a minimum of v in M. If x is unique or if v is strictly K -monotone in M, then x ∈ E(M, K).

Proof: If v is strictly monotone, the result follows with Theorem 2.2.1. For the rest of the proof, see, e. g., [45, Satz 2.20]. □

As long as only one efficient point is searched for, one may therefore choose an arbitrary (strictly) monotone function and minimize this function using standard optimization techniques. This strategy is called scalarization, and the function v is called value function. (More generally, any kind of function (monotone or not) mapping points in M onto real numbers is sometimes called a value function.) If more than one efficient point is searched for, one has to use several monotone functions. Moreover, for reasons of efficiency these functions should be as simple as possible. Clearly, linear forms in K^{∗} are monotone, but for reasons which will become clear in the next section, one might be interested in using slightly more complicated functions. For example, suppose that Theorem 2.2.2 is employed with the most simple strictly convex functions monotone on IR_{+}^{n}, i. e. g_{i}(x_{i}) := x_{i}^{2} (i = 1, …, n), and that after this transformation a linear form ω ∈ K^{∗} is used to scalarize the multicriteria problem with feasible set g(M). The result can be written as a scalar problem of the form min_{x∈M} ∑_{i=1}^{n} ω_{i}x_{i}^{2}. The net result is the use of a quadratic value function instead of a linear one. It is therefore important to discuss the characteristics of value functions of this type.

Theorem 2.2.4 Let K be a closed convex cone in IR^{n}, let M ⊆ IR^{n} be a nonempty convex set and L be the subspace parallel to aff(M). Let Q ∈ IR^{n×n} be a symmetric positive semidefinite matrix and define the function v by v(x) := ⟨ Qx, x ⟩ for all x ∈ M. Then v is monotonically increasing with respect to ≤_{K} on M if and only if

QM ⊆ K^{∗} + L^{⊥}

holds. Furthermore, if K is closed and int(K^{∗}+L^{⊥}) ≠ ∅, then v is strictly monotonically increasing with respect to ≤_{K} on relint(M) if and only if

Qrelint(M) ⊆ relint(K^{∗}) + L^{⊥}

holds.

Proof: Since ∇v(x) = 2Qx for all x, the convex function v is monotonically increasing on M if and only if ⟨ Qx, z ⟩ ≥ 0 for all z ∈ K ∩ (−x + M) and x ∈ M. Now let x ∈ relint(M). Then −x + M is a neighborhood of 0 in L with respect to the induced topology. As a consequence, ⟨ Qx, z ⟩ ≥ 0 holds for all z ∈ K ∩ (−x + M) if and only if ⟨ Qx, z ⟩ ≥ 0 holds for all z ∈ K ∩ L. This is equivalent to Qx ∈ (K ∩ L)^{∗} = K^{∗} + L^{⊥}, see [77, Corollary 16.4.2]. Since K^{∗} + L^{⊥} is closed, the relation Qrelint(M) ⊆ K^{∗} + L^{⊥} is equivalent to QM ⊆ K^{∗} + L^{⊥}.

Using the same way of reasoning, the function v is strictly monotonically increasing on the set relint(M) if and only if ⟨ Qx, z ⟩ > 0 for all z ∈ (K \ {0}) ∩ (−x + relint(M)) and x ∈ relint(M). This, in turn, is equivalent to Qrelint(M) ⊆ (K ∩ L)^{+}. Since K is closed and int(K ∩ L)^{∗} = int(K^{∗} + L^{⊥}) ≠ ∅, we have by using [45, Satz 2.7 c)] and [77, Corollary 16.4.2] that (K ∩ L)^{+} = int(K^{∗} + L^{⊥}). Invoking [77, Corollary 6.6.2], we find that (K ∩ L)^{+} = relint(K^{∗}) + L^{⊥}. □

Depending on the situation at hand, the identity

K^{∗} + L^{⊥} = (K ∩ L)^{∗}

might become useful.

The assumption int(M) ≠ ∅ is rather weak, since E(M, K) = E(M + K, K), and the assumption int(K) ≠ ∅ usually holds in a multicriteria framework. Therefore, the relation Q(M + K) ⊆ K^{∗} might be used to check monotonicity of v on M.

Corollary 2.2.5 Let K be a closed convex cone in IR^{n} and let Q ∈ IR^{n×n} be a symmetric positive semidefinite matrix. Define the function v by v(x) := ⟨ Qx, x ⟩ for all x ∈ IR^{n}. If Qx ∈ K^{+} for an x ∈ IR^{n}, then there exists a neighborhood U of x such that v is monotonically increasing with respect to ≤_{K} on U.

Slightly more complicated candidates for monotone functions are norms on IR^{n}. However, with these monotonicity on IR^{n} can be achieved only if {0} ≠ K ⊆ K^{∗}, see [55, p. 110] and [45, p. 81]. Fortunately, monotonicity is needed only on M, not on IR^{n}. The scalar optimization problems induced by choosing norms or slightly more general value functions will be further investigated in Section 2.2.3.3 and Chapter 4.

2.2.2 Proper Efficiency

It can be said that there are very few optimization methods for continuous multicriteria optimization problems which do not use in some way a scalarization of the problem at hand. With scalarization it is usually meant a rephrasing of the multicriteria problem into a single-criterion one. The following definition and the accompanying theorem can be used to achieve such a rephrasing.

Definition 2.2.2 (Proper Efficiency) Let K ⊆ IR^{n} be a cone and M ⊆ IR^{n}. A point x ∈ M is called properly efficient with respect to K, if there exists a point ω ∈ K^{+} such that x is a solution to the problem

minimize ⟨ ω, y ⟩
subject to y ∈ M.

(2.1)

The set of properly efficient points of M with respect to K will be denoted by E_{p}(M, K).

Because linear functionals ω ∈ K^{∗} are monotone, the next result follows immediately.

Theorem 2.2.6 For K, M ⊆ IR^{n} the following statements hold:

1. Let K be a convex cone with 0 ∈ K and {0} ≠ K ≠ IR^{n}. Then

E_{p}(M, K) ⊆ E(M, K).

2. Let K be a closed convex cone with 0 ∈ K such that K contains no lines. Let M be closed and convex. Then

E(M, K) ⊆ cl(E_{p}(M, K)).

3. Let K = IR_{+}^{n} and M polyhedral. Then

E(M, K) = E_{p}(M, K).

Proof: Different proofs can be found in various textbooks and original articles. See, e. g. [45, Satz 2.22], [83, p. 74], or [21]. The first proof of Part 2 is probably due to Arrow, Barankin, and Blackwell [3]. □

From the above result it is clear that under the assumptions given in the theorem we might focus our attention on the set of properly efficient points. In the case of item 2 of the theorem, one has E_{p}(M, K) ⊆ E(M, K) ⊆ cl(E_{p}(M, K)). But in a numerical solution method one will usually not be able to distinguish between a subset of IR^{n} and it’s closure. As a consequence, E_{p}(M, K) is a good approximation to E(M, K). Moreover, there exist other characterizations of the set of proper efficient points [43], by which it is clear that the set of proper efficient points consists of those efficient points for which an unbounded trade-off in the function values is not possible. Note that the main assumption needed is the convexity of M. Other scalarizations exist which do not need this convexity assumption [56]. However, the corresponding objective function is no longer a simple linear form.

2.2.2.1 Several Objective Functions

The discussion above has shown that we might concentrate on the set E_{p}(M, K), and the next theorem gives just one more motivation.

Theorem 2.2.7 Let G ⊆ IR^{m} be convex and let f_{i} : G ⟶ IR be strictly convex (i = 1, …, n). Define f := (f_{1}, …, f_{n})^{⊤} and M := f(G). Then x ∈ E(M, IR_{+}^{n}) if and only if there exists an ω ∈ IR_{+}^{n} \ {0} such that x is the solution to the problem

minimize ⟨ ω, y ⟩
subject to y ∈ M.

Proof: [45, Satz 2.26] □

The scalarization discussed in the last section is a linear one, i. e. the multiple objective problem was replaced by a problem with one linear objective function, while the set of feasible points stayed the same. In the case of Section 2.1.2 and Theorem 2.2.7, the proper problem consists not only in finding an efficient point z ∈ M = F (G), but also of finding a point x ∈ f^{−1}(z) in the corresponding inverse image. This however, can be achieved by introducing artificial variables t_{i} = f_{i}(x) (i = 1, …, n). Let a scalarization vector ω ∈ K^{+} be given. Instead of the scalarized problem (2.1) of Definition 2.2.2, we can now consider the problem

minimize ⟨ ω, t ⟩
subject to t = f(x),
x ∈ G.

(2.2)

Obviously, each solution of (2.2) induces a solution of (2.1) in a trivial way and vice versa. Moreover, if ω ∈ IR_{+}^{n}, we can replace the equality constraints t = f(x) by much simpler inequality constraints of the type (x, t) ∈ epi(f).

2.2.3 Approximating the Set of Efficient Points

By writing

E_{p}(M, K) = ⋃_{ω∈K^{+}} arg min{⟨ ω, y ⟩ | y ∈ M}  (2.3)

and noting Theorem 2.2.6 and the accompanying discussion, one sees that approximating the set of efficient points is possible by discretizing the set K^{+} and solving a finite number of scalarized problems. Now let a base of K^{+} be given, i. e. a convex nonempty set B ⊂ K^{+} such that for each ω ∈ K^{+} there exists a unique λ ≥ 0 and a unique β ∈ B with ω = λβ. Then (2.3) simplifies to

E_{p}(M, K) = ⋃_{ω∈B} arg min{⟨ ω, y ⟩ | y ∈ M}.  (2.4)

This is an especially important representation if B is bounded, because then only the bounded set B needs to be discretized.

Unfortunately, choosing a specific uniform discrete approximation for B will not result in a uniform approximation of E_{p}(M, K), as computational tests [18, 17] and the following simple example show.

Example 2.2.1 Let n = 2 and M := conv({(1, 0)^{⊤}, (0, 1)^{⊤}}) as well as K := IR_{+}^{2}. Obviously, E(M, K) = M. The standard base for K^{+} = int(K) is B = relint(M). For ω = (1/2, 1/2)^{⊤}, any point in M is a solution to the problem min_{y∈M}⟨ ω, y ⟩. For all other ω ∈ B, the solution to the corresponding scalarized problems is either (1, 0)^{⊤} or (0, 1)^{⊤}. As a consequence, most discretizations of B will result in the poor approximation {(1, 0)^{⊤}, (0, 1)^{⊤}} of M, and increasing the discretization size will usually not be helpful.

Note that this example can easily be generalized to higher dimensions.

Taking a closer look at the example, one of the main problems occurring when using linear scalarizations in a naive way can be summarized as follows. Let C ⊂ B be a dense subset of a base B of K^{+}. For each ω ∈ C, let a solution x(ω) of the problem min_{y∈M}⟨ ω, y ⟩ be given. Then it does not follow that {x(ω) | ω ∈ C} is dense in E_{p}(M, K). Adding the assumption of uniqueness of x(ω) for all ω ∈ C does not help, either. Likewise, if C_{1}, C_{2} ⊆ B with C_{1} ⊂ C_{2} such that C_{2} is a ”better” approximation to B, then it does not follow that ⋃_{ω∈C_{2}} arg min{⟨ ω, y ⟩ | y ∈ M} is a ”better” approximation to E_{p}(M, K) than ⋃_{ω∈C_{1}} arg min{⟨ ω, y ⟩ | y ∈ M}.

Other methods [17] designed to approximate the set of efficient points avoid the nonuniformity drawback depicted above, but fail to approximate the set of efficient points in an even more important respect [17, Figure 3].

2.2.3.1 Linear Scalarizations

Theorem 2.2.2 points the way to sufficient conditions which ensure that a proper discretization of a base B of K^{+} leads to a good approximation of E(M, K). Let M be closed, nonempty, and convex and let all facets of M have at most one point with E_{p}(M, K) in common. Then, for each ω ∈ B the function v with v(x) := ⟨ ω, x ⟩ has a unique minimum in E(M, K). Moreover, the set-valued optimal argument function

φ : ω ⟼ arg min{⟨ ω, y ⟩ | y ∈ M}

can be identified in the standard way with a function φ : B ⟶ M. Note that φ(B) = E_{p}(M, K).

Lemma 2.2.8 (Continuity of the Optimal Argument Function) Let the set M be convex, closed, and nonempty and let B be a base of K^{+}. Suppose that all facets of M have at most one point with E_{p}(M, K) in common. Then, the following holds:

1. The function φ : B ⟶ M is continuous.

2. There exists a partition (B_{i})_{i∈I} of B of nonempty convex sets B_{i} ⊆ B, compact relative to the induced topology on B, such that the sets relint(B_{i}) are mutually disjoint and such that for all i ∈ I the function φ is constant on B_{i}. Moreover, the family (φ(B_{i}))_{i∈I} consists of mutually disjoint sets.

Proof:

1. Suppose that there exists a sequence (ω_{k})_{k∈IN} such that ω_{k} ∈ B (k ∈ IN) and lim_{k→∞} ω_{k} = ω ∈ B, but lim_{k→∞} φ(ω_{k}) ≠ φ(ω). Then, there exists a converging subsequence of (φ(ω_{k}))_{k∈IN}, for otherwise there would exist a facet F of M such that F ∩ E_{p}(M, K) has dimension 1 or higher. Without loss of generality, we may assume lim_{k→∞} φ(ω_{k}) =: x ∈ M. Then x as well as φ(ω) are solutions to the problem min_{y∈M}⟨ ω, y ⟩, a contradiction.

2. Use I := E_{p}(M, K) and define B_{i} := φ^{−1}(i). The compactness follows from standard arguments. The convexity of B_{i} follows by noting that ω ∈ B_{i} if and only if −ω ∈ N_{M}(i). □

Note that the uniqueness assumption of points in a base is not used. Only compactness and convexity of the set B is needed to arrive at the assertions of the lemma.

Remark 2.2.9 In the context of item No. 2 of the last lemma, let J ⊆ I be the set of all indices j for which B_{j} contains exactly one point and define C := ⋃_{j∈J} B_{j}. If B is discretized according to some discretization scheme producing a finite number of points, it is obviously unnecessary to consider more than one point from each set B_{i} with i ∈ I \J. Instead, one has to pay attention to the set C. Fortunately, φ is one-to-one on C. The problem of identifying the boundary between C and B \ C remains.

Corollary 2.2.10 Let M be a convex, compact and nonempty set such that all facets of M have at most one point with E_{p}(M, K) in common. Let B be a base of K^{+} and C be a dense subset of B. Then, φ(C) ⊆ E_{p}(M, K) is dense in E_{p}(M, K).

Proof: Let x ∈ E_{p}(M, K). There exists an ω ∈ B such that x is the unique solution of the problem min_{y∈M}⟨ ω, y ⟩. Let (ω_{k})_{k∈IN} be a sequence in C with ω as an accumulation point. This sequence induces a unique sequence (φ(ω_{k}))_{k∈IN} of solutions to the corresponding problems with φ(ω_{k}) ∈ E_{p}(M, K) (k ∈ IN). According to the last lemma, x is the limit of this sequence. □

Note that the function φ is not Lipshitz-continuous, as the example K := IR_{+}^{2}, M := \{(x, −\sqrt{x})^{⊤} | x ∈ IR_{+}\} + IR_{+}^{2} shows.

Nevertheless, at least for K ⊇ IR_{+}^{n}, we can always use Theorem 2.2.2 to construct a one-to-one mapping g : E(M, K) ⟶ E(g(M), K) such that the density argument holds for E(g(M), K). Since g is one-to-one and continuous, a good approximation of E(M, K) can be constructed by discretizing a base of K^{+} and solving parameterized optimization problems with linear objective function for the feasible set g(M). In these optimization problems, the parameters will only enter in the objective function, and not in the constraints. This will be exploited in the next chapter.

2.2.3.2 Quadratic Scalarizations

The discussion on p. 14 has shown that nonlinear scalarizations are useful in their own right. Moreover, such scalarizations can often be interpreted as a nonlinear transformation followed by a linear scalarization. In what follows, quadratic value functions are discussed.

Theorem 2.2.11 (Quadratic Scalarizations) Let K be a convex cone and M be a convex set with nonempty interior. For an arbitrary set of matrices B ⊆ IR^{n×n} define

A(M, B) := ⋃_{Q∈B} arg min{⟨ Qy, y ⟩ | y ∈ M}.  (2.5)

1. Let

Q_{M,K} := {Q ∈ IR^{n×n} positive definite| QM ⊆ K^{∗}}.

Then

A(M, Q_{M,K}) ⊆ E(M, K).

2. Let K^{+} ⊆ int(IR_{+}^{n}) and M ⊂ int(IR_{+}^{n}). Define

Q̂_{M,K} := {diag(ω_{1}/z_{1}, …, ω_{n}/z_{n}) | z ∈ M, ω ∈ K^{+}}.

Then

E_{p}(M, K) ⊆ A(M, Q̂_{M,K}).

Proof:

1. Apply Theorem 2.2.4 and 2.2.3.

2. Let x ∈ E_{p}(M, K). Then there exists an ω ∈ K^{+} such that x is solution to min_{y∈M}⟨ ω, y ⟩. Since M is convex, this is equivalent to −ω ∈ N_{M}(x), where N_{M}(x) is the normal cone to M at the point x. Define q_{i} := ω_{i}/(2x_{i}) > 0 and Q := diag(q). Then −∇f(x) = −2Qx ∈ N_{M}(x), where f(y) := ⟨ Qy, y ⟩. As a consequence, x is a minimum of the convex program

minimize ⟨ Qy, y ⟩
subject to y ∈ M.

(2.6)

□

Corollary 2.2.12 Let K be a closed convex cone with 0 ∈ K, K^{+} ⊆ int(IR_{+}^{n}) such that K contains no lines. Let M be closed and convex with nonempty interior and M ⊂ int(IR_{+}^{n}). Then

cl(A(M, Q_{M,K})) ⊆ cl(E(M, K)) = cl(E_{p}(M, K)) ⊆ cl(A(M, Q̂_{M,K})).

Obviously, the situation of highest interest occurs whenever

Q̂_{M,K} ⊆ Q_{M,K}  (2.7)

as well as 0 ∈ K, {0} ≠ K ≠ IR^{n} holds, because then

E_{p}(M, K) ⊆ A(M, Q̂_{M,K}) ⊆ E(M, K).

To give some sufficient condition for this, suppose that M ⊂ int(IR_{+}^{n}) and define M^{−1} := {(1/z_{1}, …, 1/z_{n})^{⊤} | z ∈ M}. Moreover, for two points x, y ∈ IR^{n} define x · y := (x_{1}y_{1}, …, x_{n}y_{n})^{⊤}. This notation is extended to subsets of IR^{n} in the usual way, i. e. S_{1} · S_{2} := {x · y | x ∈ S_{1}, y ∈ S_{2}}. With this, the relation

K^{+} · M^{−1} · M ⊆ K^{∗}  (2.8)

is sufficient for (2.7).

Corollary 2.2.13 Let K = IR_{+}^{n}. Moreover, let M be convex with nonempty interior and M ⊂ int(IR_{+}^{n}). Then

E_{p}(M, K) ⊆ A(M, Q̂_{M,K}) ⊆ E(M, K).  (2.9)

Unfortunately, if int(M) ≠ ∅, the relation (1, …, 1)^{⊤} ∈ int(M^{−1} · M) holds, and therefore (2.8) cannot hold for cones K with K^{∗} ⊂ IR_{+}^{n}, K^{∗} ≠ IR_{+}^{n}, K^{+} = int(K^{∗}).

In the context of the last corollary, the set A(M, Q̂_{M,K}) is, as an approximation to E(M, K), at least as good as E_{p}(M, K). From a numerical point of view, stronger results than (2.9) are of no interest, since the three sets in (2.9) can not be distinguished numerically anyway. Due to the scale-invariance of the problems (2.6) one might concentrate on the standard base of int(IR_{+}^{n}) and use

S := {diag(q) |q ∈ int(IR_{+}^{n}), ∑_{i=1}^{n} q_{i} = 1}

instead of Q̂_{M,K}.

In the general case of K ≠ IR_{+}^{n}, K^{+} ⊂ IR_{+}^{n} there remains the difficulty of the rather implicit definition of Q_{M,K}. Moreover, a naive parameterization of Q̂_{M,K} would involve 2n parameters from a rather large set. Smaller parameter sets can only be used if a base of cone(M^{−1} · K^{+}) is known. A typical strategy would therefore be first to calculate points in a superset of A(M, Q̂_{M,K}), and, then, reject those points not belonging to E(M, K). For example, compute a point x ∈ A(M, S) by solving (2.6) with a matrix Q ∈ S. Define ω := 2Qx ∈ −N_{M}(x). If ω ∈ K^{+}, then x ∈ E_{p}(M, K). If ω ∉ K^{+} and N_{M}(x) is one-dimensional, then x ∉ E_{p}(M, K). Note that N_{M}(x) is one-dimensional as soon as M is strictly convex and x ∈ bd(M). The latter is always the case for solutions to (2.6) if 0 ∉ M. If ω ∉ K^{∗}, there exists a z ∈ −K with ⟨ ω, z ⟩ > 0. This means that z is a direction of increase at x for the objective function v of (2.6). But then the function v is not monotone with respect to ≤_{K} in a neighborhood of x, and Theorem 2.2.3 can not be applied to check the efficiency of x.

There is a clear overhead in computation here in the sense that first E(M, IR_{+}^{n}) is approximated and then points of the approximation are rejected in order to calculate an approximation to the set E_{p}(M, K), K^{+} ⊂ IR_{+}^{n}. Nevertheless, the possibility to compute simultaneously an outer and an inner approximation to the set of efficient points may be interesting.

Note that each minimum of the scalar problems (2.6) is unique, due to the strict convexity of the objective functions considered. Moreover, the objective function is coercive. This means that if M is closed, convex, and nonempty, the set-valued optimal argument function

ψ : Q ⟼ arg min{⟨ Qy, y ⟩ | y ∈ M}

for the quadratic scalarizations considered in this section can be identified in the standard way with a function ψ : P (n) ⟶ M, where P (n) is the set of all positive definite n × n-matrices.

Theorem 2.2.14 (Lipshitz-Continuity of the Optimal Argument Function) Let M ⊂ IR^{n} be closed, convex, and nonempty. Then, the function ψ : P (n) ⟶ M is locally Lipshitz-continuous.

Proof: The result holds due to [9, Proposition 4.32]. □

Note that the continuity of ψ follows already from the fact that each objective function of (2.6) is strictly convex for Q ∈ P (n).

Corollary 2.2.15 Let M ⊂ IR^{n} be closed, convex, and nonempty. Let B ⊆ P (n) and R be a dense subset of B. Then A(M, R) is dense in A(M, B).

Proof: The proof follows the same lines as for Corollary 2.2.10 and is therefore omitted here. □

The linear and quadratic scalarizations presented in this and the last subsection have one interesting quality, especially from a numerical point of view: the objective functions are infinitely often differentiable with respect to their arguments from M and their parameters. This makes it especially simple to use them in a numerical algorithm, a subject to which we will return to in Chapter 3. Other scalarization techniques do not necessarily lead to smooth objective functions. The next subsection shows that smooth reformulations of the corresponding scalar optimizations problems are sometimes possible. However, it is not clear how efficient standard solution techniques are when applied to these reformulations. In the next chapter, we will therefore consider only the linear and quadratic scalarizations discussed previously.

2.2.3.3 Norm Scalarizations

All the scalarizations investigated above can be interpreted as special cases of goal programming scalarizations, i. e. a scalarizations where norms (or their non-symmetric counterpart, gauges) are used. This is a trivial observation for problems in the context of Theorem 2.2.11 and its corollaries, but even for (2.3) we simply note that in the case of K^{+}, M ⊆ IR_{+}^{n} the norms in use are of the form

γ_{ω}(x) := ∑_{i=1}^{n} ω_{i}|x_{i}|

with ω ∈ K^{+}. Clearly, γ_{ω}(x) = ⟨ ω, x ⟩ for all x ∈ M. Similar scalarizations different from the ones described above can also be used. We will return to this subject in Chapter 4, cmp. also [56, 12]. Especially useful might be norm scalarizations of the form

min ∥Qx∥_{p}  (2.10)
subject to x ∈ M

with 1 ≤ p ≤ +∞ and Q ∈ P (n). Clearly, (2.10) is a generalization of (2.6). By a proper choice of Q, the objective function of (2.10) can be made monotonically increasing. (For example, if K = IR^{n} one can choose an arbitrary diagonal matrix.) On the other hand, scalarizations like (2.10) with p ≠ 2 might be better suited to generate efficient points of nonconvex sets M. However, the objective function will be nonsmooth, which usually results in serious difficulties when one tries to solve the scalarized problem numerically. Therefore, a smooth reformulation of the objective function might be helpful. To this end, consider the standard idea of replacing (2.10) by

min t
subject to x ∈ M,  (2.11)
(Qx, t) ∈ K_{p},

where

K_{p} := {(y, t) ∈ IR^{n+1} | ∥y∥_{p} ≤ t}.

is the epigraph of the function ∥ · ∥_{p}. Now the parameters Q occur in the constraints of the program, a situation which is better avoided, as the discussion in Chapter 3 will show. For p = 1 or p = +∞, the cone-constraints are linear. We will discuss the case 1 < p < +∞ further in Chapter 4, especially Example 4.5.5.

2.2.4 The Nonconvex Case

It is clear that with the tools introduced in the last section, only locally efficient points can be characterized if the set of feasible points M is not convex. Computing globally efficient points is, of course, a generalization of single-criterion optimization. As such, the problem posed in it’s general form is N P -complete, and therefore there is few hope for efficient solution methods. However, special cases, especially those of low dimension, might be solved efficiently. In what follows, only methods able to find locally efficient points are discussed, albeit briefly. If m discrete points are generated as a candidate set for (global) efficiency, checking for pairs of points dominating each other can be done in time O(m^{2}), and this is automatically done by any method for discrete multicriteria problems.

Theorem 2.2.3 also holds for strictly monotone functions v: local minima of v are locally efficient. The strictly monotone quadratic functions discussed in this chapter can therefore be used to generate locally efficient points.

Note, however, that the set of globally efficient points may not be connected. Even isolated points may occur. This is illustrated by the following simple example.

Example 2.2.2 For n = 2 let K = IR_{+}^{2} and let B(r) := {z ∈ IR^{n} : ∥z∥_{∞} ≤ r} denote the ball of the max-norm with radius r. Define

M := ((5/2, 3/2)^{⊤} + B(1/2)) ∪ ((3/2, 5/2)^{⊤} + B(1/2)) .

For this nonconvex set, the set of (locally as well as globally) efficient points is E(M, K) = {(1, 2)^{⊤}, (2, 1)^{⊤}}.

Nevertheless, the following theorem holds.

Theorem 2.2.16 Let K be a convex cone and M be closed such that E(M, K) ⊆ K. Then

E(M, K) ⊆ cl(A(M, P (n)))

holds.

Proof: Let x ∈ E(M, K). Then M ∩ (x − K) = {x} and 0 ∈ x − K. Let (Q_{k})_{k∈IN} be a sequence of positive definite matrices such that Q_{k} has the eigenvector x to the eigenvalue 1/∥x∥_{2} for all k ∈ IN and all other eigenvalues diverge to +∞ for k ⟶ +∞. Let (x_{k})_{k∈IN} be a sequence with x_{k} ∈ arg min{⟨ Q_{k}y, y ⟩ | y ∈ M} for all k ∈ IN. Then (x_{k})_{k∈IN} converges to x. □

Heuristically, the problems constructed in the proof ”converge” to a problem equivalent to

min ∥tx∥_{2}^{2}
subject to tx ∈ M,

whose only solution is t = 1. The assumption E(M, K) ⊂ K holds, for example, when IR_{+}^{n} ⊆ K and M ⊂ int(IR_{+}^{n}).

Equality between the set of efficient points and the set of solutions of a family of parameterized scalar optimization problems can be shown even in the nonconvex case [56]. However, the scalar optimization problems considered are nonsmooth, and it is not clear how these problems can be solved efficiently or how a smooth reformulation should look like.

Note that the set E(M, K) \ A(M, P (n)) might be nonempty, as can be seen by modifying Example 2.2.2 as follows.

Example 2.2.3 For n = 2 let K = IR_{+}^{2} and let B_{2}(r) := {z ∈ IR^{n} : ∥z∥_{2} ≤ r} denote the ball of the Euclidean norm with radius r. Define

M := ((2, 1)^{⊤} + B_{2}(1)) ∪ ((1, 2)^{⊤} + B_{2}(1)) .

Then E(M, K) is the lower left boundary of M joining the points (1, 0)^{⊤}, (1, 1)^{⊤}, and (0, 1)^{⊤}, but there does not exist a positive definite matrix Q such that (1, 1)^{⊤} is solution to the problem (2.6). Moreover, for a given base B of P (2) one needs a dense subset R ⊂ B in order to guarantee that A(M, R) is dense in E(M, K).

On the other hand, in Example 2.2.2 the set A(M, P (2)) consists of the set of weakly efficient points, which can be written as [a, b] ∪ [b, c] ∪ [c, d] ∪ [d, e] with a = (1, 3)^{⊤}, b = (1, 2)^{⊤}, c = (2, 2)^{⊤}, d = (2, 1)^{⊤}, and d = (3, 1)^{⊤}, respectively.

26

Interior-Point Methods for Convex Multicriteria Problems

Standard interior-point methods are concerned with the minimization of a function f : G ⟶ IR over a set G ⊆ IR^{m}. Based on the discussion in the last chapter, in the context of multicriteria programming this function f will be parameterized according to n parameters q ∈ IR^{n}. If the multicriteria scalarization from Theorem 2.2.11 is used, we can assume that q ∈ S = {q ∈ int(IR_{+}^{n}) | ∑_{i} q_{i} = 1} and f(x) = ⟨ diag(q)x, x ⟩. In case a linear scalarization f(x) = ⟨ q, x ⟩ is used, q ∈ K^{+} holds. If several objective functions f_{i} : G ⟶ IR (i = 1, …, n) are given, we can use either f(x) = ∑_{i=1}^{n} q_{i}(f_{i}(x))^{2} resp. f(x) = ∑_{i=1}^{n} q_{i}f_{i}(x) or the reformulation

minimize ⟨ diag(q)t, t ⟩
subject to f_{i}(x) ≤ t_{i}  (i = 1, …, n)
x ∈ G.

(3.1)

resp.

minimize ⟨ q, t ⟩
subject to f_{i}(x) ≤ t_{i}  (i = 1, …, n)
x ∈ G.

(3.2)

Accordingly, we can assume that the objective function f is linear or quadratic and the set of feasible points is fixed, i. e. independent of the parameterization. This is an important simplification in order to employ to interior-point methods, since quadratic objective functions are 0-compatible with self-concordant barrier functions [67] and constant sets of feasible points reduce the amount of work needed to find strictly feasible starting points.

The rest of this chapter is as follows: in Section 3.1, the standard primal path-following interior-point algorithm as presented in [67] is rehearsed in short. It turns out that the first stage of the algorithm (the construction of a point close to the ”center” of the set of feasible points) can be used without any modification in a multicriteria framework. In Section 3.2, an algorithm is presented in which two scalar linear programs are solved simultaneously by way of an adapted interior-point technique. Moreover, an efficiency estimate will elucidate the dependence of the number of steps of the algorithm on the perturbation size. The corresponding algorithmic technique will be generalized to quadratic scalar programs in Section 3.3 and to smooth convex scalar programs in Section 3.4. Since, on the one hand, the discussion in Section 2.2.1 has shown that only linear scalarizations need to be considered when solving multicriteria problems, and, on the other hand, it will turn out that there exist fundamental differences between linear and quadratic scalarizations when solving multicriteria problems, only linear scalarizations will be considered from Section 3.5 onwards. In this Section 3.5, bicriteria problems are considered, and a general algorithm based on a refinement strategy in the parameter space is proposed. It turns out that the theoretical complexity of this algorithm is lower than the naive approach of using a prespecified number of parameter values to achieve a certain density of approximation in the image space IR^{n} of the multicriteria optimization problem considered. Finally, the results from Section 3.5 are extended to the case of more than two criteria in Section 3.6. There, a general algorithm for approximating the set of efficient points of a smooth convex multicriteria optimization problem is discussed.

3.1 The Standard Primal Path-Following Interior-Point Algorithm

This section recalls in short how the standard primal path-following algorithm presented by Nesterov and Nemirovskii [67, Chapter 3] works. Moreover, the notation used in the rest of this chapter is introduced, together with some slight generalizations and corrections. For most of the proofs of the result presented here, the reader is referred to [67].

3.1.1 The First Stage

Let G ⊆ IR^{m} be a compact convex set with nonempty interior and F : int(G) ⟶ IR a self-concordant barrier for G with self-concordancy parameter ϑ ≥ 1. Moreover, let f : G ⟶ be a β -compatible objective function. In a scalarized multicriteria setting, the discussion of the last chapter has shown that it makes sense to assume that f is linear or at most quadratic. Moreover, a nonquadratic function can be replaced in the usual way by a linear one. We can therefore assume β = 0.

It is known that for compact sets G the Hessian ∇^{2}F (x) at an arbitrary point x ∈ int(G) is positive definite. This allows us to define for x ∈ int(G), the Newton decrement of F at x by

λ(F, x) := (⟨∇F (x)), (∇^{2}F (x))^{−1}∇F (x)⟩)^{1/2} ,

cmp. [67, p. 15–16]. The Newton-decrement can be viewed as a measure for the length of ∇F (x). More precisely, it is the length of the gradient in the norm induced by ∇^{2}F (x). Therefore, the smaller λ(F, x) is, the closer x is to a minimum of F.

Define the auxiliary function ω :]0, 1/3[⟶ IR by

ω(λ) := 1 − (1 − 3λ)^{1/3}.

In the standard path-following scheme, there is needed a 5-tuple

(λ_{1}, λ′_{1}, λ_{2}, λ_{3}, λ′_{3}) ∈ IR^{5}

basically controlling the neighbourhoods to the central path (to be defined below) used, such that the following system of inequalities holds:

0 < (\frac{λ_{i}}{1 − λ_{i}} )^{2} ≤ λ′_{i} < λ_{i} < 2 − \sqrt{3}  (i = 1, 3),  (3.3)
λ′_{1} < λ_{2} < λ_{3},  (3.4)
\frac{ω(λ′_{3})}{1 − ω(λ′_{3})} < 1,  (3.5)
\frac{ω(λ_{2})}{1 − ω(λ_{2})} < \frac{1}{3},  (3.6)
\frac{ω^{2}(λ′_{1})(1 + ω(λ′_{1}))}{1 − ω(λ′_{1})} ≤ \frac{1}{9},  (3.7)

see [67, p. 70]. But (3.5) is equivalent to λ′_{3} ≤ 7/24, while (3.6) is equivalent to λ_{2} < 37/192. Moreover, some simple calculations show that (3.7) is equivalent to

λ′_{1} ≤ \frac{1}{3} − \frac{1}{3} (\frac{4}{3} − \frac{1}{9} ^{3} \sqrt{27 + 3 \sqrt{57}} − \frac{2}{3} \frac{1}{\sqrt[3]{27 + 3 \sqrt{57}}})^{3} =: ϱ.

Due to ϱ > 37/192, (3.4) tells us that we can drop (3.7). As a consequence, the system (3.3)–(3.7) can be simplified to

0 < λ′_{i} < λ_{i} < 2 − \sqrt{3}  (i = 1, 3),  (3.8)
(\frac{λ_{1}}{1 − λ_{1}} )^{2} ≤ λ′_{1},  (3.9)
(\frac{λ_{3}}{1 − λ_{3}} )^{2} ≤ λ′_{3} < \frac{7}{24},  (3.10)
λ′_{1} < λ_{2} < λ_{3},  (3.11)
λ_{2} < \frac{37}{192}.  (3.12)

The first stage of the interior-point method starts at an arbitrary point y^{(0)} ∈ int(G) and calculates an approximation y^{(j)} to the unique minimizer of F such that the Newton decrement at y^{(j)} is bounded by λ_{2}, i. e.

λ(F, y^{(j)}) ≤ λ_{2}.  (3.13)

In these iterations, only standard Newton steps for the function F are performed. The number of iterations j needed to achieve this is bounded from above by

1 + \frac{λ_{1} + \sqrt{ϑ}}{λ_{1} − λ′_{1}} (ln (\frac{ϑ}{a(y^{(0)}, G)}) + ln (\frac{2}{λ_{2} − λ′_{1}} )) ,

see [67, Proposition 3.2.3], and only the parameters λ_{1}, λ′_{1}, and λ_{2} are used within the algorithm. In the run-time estimate above, a(y^{(0)}, G) is the asymmetry coefficient of y^{(0)} with respect to G, defined for arbitrary points x ∈ int(G) by

a(x, G) := sup{α ≥ 0 | x + α(x − G) ⊆ G}.

For lower bounds of this number, the reader is referred to Section 4.5. Obviously, for large ϑ the leading term \sqrt{ϑ} ln ϑ in the estimate above has the coefficient 1/(λ_{1} − λ′_{1}), so it makes sense to maximize λ_{1} − λ′_{1}. Within the bounds given above, the maximum is taken at

λ_{1} = −1/3 ^{3} \sqrt{27 + 3 \sqrt{105}} + 2 \frac{1}{\sqrt[3]{27 + 3 \sqrt{105}}} + 1 = 0.2290 …

and λ′_{1} = λ_{1}^{2}/(1−λ_{1})^{2} = 0.0882 . . .. The leading coefficient then becomes 1/(λ_{1}−λ′_{1}) = 7.103 . . ..

In a multicriteria setting, the computation of an approximation of the analytic center of G can be executed without any modification at all. If we do not have any further information on the parameters of the scalarization, no better starting point for minimizing the scalarized problems is at hand.

3.1.2 The Second Stage

In the second stage, we start at the just obtained x^{(0)} := y^{(j)} and follow the central path to a minimum of f by computing discrete approximations x^{(i)} ∈ G (i ∈ IN) to this central path. More precisely, we consider the function

F_{t}(x) := tf(x) + F (x)

as well as the central path

{x(t) ∈ G | λ(F_{t}, x(t)) = 0; t > 0} = {x(t) ∈ G | x(t) minimum of F_{t}; t > 0}

and compute parameter values t(i) > 0 with lim_{i→∞} t(i) = ∞ as well as

x^{(i+1)} := x^{(i)} − (∇_{x}^{2}F_{t(i+1)}(x^{(i)}))^{−1} ∇_{x}F_{t(i+1)}(x^{(i)}), i = 0, 1, 2, ….  (3.14)

The parameter t plays the role of the barrier parameter and is increased in each step. The generated sequence (x^{(i)})_{i∈IN} is a discrete approximation to the central path of the problem considered.

The sequence of barrier parameters can be defined by

t(i) := t(1) exp ((i − 1)\frac{λ_{3} − λ′_{3}}{λ_{3} + \sqrt{ϑ}}) ,  i = 1, 2, ….  (3.15)

We will take a look at the derivation of this rule in Lemma 3.1.4. Of course, an initialization t(1) of the sequence is needed. A proper way for initializing this starting value is discussed now.

For a given x ∈ int(G), define the norm ∥ · ∥_{x,F} by

∥y∥_{x,F} := (⟨y, ∇^{2}F (x)y⟩)^{1/2} .

(Recall that ∇^{2}F (x) is positive definite.) With this norm, the Dikin ellipsoid of F centered at x ∈ int(G) with radius r is defined by

W (x, r) := {y ∈ IR^{m} | ∥y − x∥_{x,F} ≤ r}.

One of the most useful relations with respect to this ellipsoid is

W (x, 1) ⊆ cl(G)  (3.16)

for all x ∈ int(G), see [67, Proposition 2.3.2]. In fact, it can be shown [75] that selfconcordancy of the function F is equivalent to this inclusion and just one further inequality which we will use to proof Lemma 3.1.2 below. As usual, the dual norm to ∥ · ∥_{x,F} will be denoted by ∥ · ∥_{x,F}^{∘}, i. e.

∥p∥_{x,F}^{∘} := sup{⟨ p, y ⟩ | ∥y∥_{x,F} ≤ 1}.

The following lemmas are well-known. Since the results are crucial in the following derivations, the proofs are given for the sake of completeness.

Lemma 3.1.1 Let there be given x ∈ int(G), μ ≥ 0, and t ≥ 0 with

t ≤ \frac{μ − λ(F, x)}{∥∇f(x)∥_{x,F}^{∘}} .

Then

λ(F_{t}, x) ≤ μ.

Proof: See also [67, Proposition 3.2.4]. We have the following chain of inequalities:

λ(F_{t}, x) = sup_{h≠0} \frac{|⟨ h, ∇F_{t}(x) ⟩|}{(⟨ h, ∇^{2}F_{t}(x)h ⟩)^{1/2}}
≤ sup_{h≠0} \frac{|⟨ h, ∇F (x) ⟩| + t|⟨ h, ∇f(x) ⟩|}{(⟨ h, ∇^{2}F (x)h ⟩)^{1/2}}
≤ λ(F, x) + t sup_{h≠0} \frac{|⟨ h, ∇f(x) ⟩|}{(⟨ h, ∇^{2}F (x)h ⟩)^{1/2}}
= λ(F, x) + t∥∇f(x)∥_{x,F}^{∘} ≤ μ.

□

As a consequence, with μ = λ_{3} > λ_{2} ≥ λ(F, x^{(0)}), a sufficiently small t(1) will always be sufficient to start with the path-following method.

The denominator of the estimate for t in Lemma 3.1.1 can be estimated itself. This is the purpose of the next two lemmas.

Lemma 3.1.2 Let there be given x ∈ int(G) with λ(F, x) ≤ δ < 1/3 and r := ω(δ)/(1 − ω(δ)) < 1/3. Denote the minimum of F by y. Then

W (x, τ) ⊆ W (y, r + τ/(1 − r))

for all τ ≥ 0.

Proof: According to [67, Theorem 2.2.2], we have

∥x − y∥_{y,F} ≤ \frac{ω(λ(F, x))}{1 − ω(λ(F, x))} ≤ \frac{ω(δ)}{1 − ω(δ)}

which means x ∈ W (y, r). Due to Theorem 2.1.1 from [67], it follows that

\frac{1 − 2ω(δ)}{1 − ω(δ)}∥h∥_{y,F} = (1 − r)∥h∥_{y,F} ≤ (1 − ∥x − y∥_{y,F})∥h∥_{y,F} ≤ ∥h∥_{x,F}

for all h ∈ IR^{m}. Now let z ∈ IR^{m}. Then

∥z − y∥_{y,F} ≤ ∥z − x∥_{y,F} + ∥x − y∥_{y,F} ≤ ∥z − x∥_{y,F} + r ≤ \frac{1}{1 − r}∥z − x∥_{x,F} + r,

from which the result follows. □

Lemma 3.1.3 Let there be given x ∈ int(G) with λ(F, x) ≤ δ < 1/3 and r := ω(δ)/(1 − ω(δ)) < 1/3. Denote the minimum of F by y and let p ∈ IR^{m} be a subgradient of a convex function f : G ⟶ IR at x. Then

∥p∥_{x,F}^{∘} ≤ \frac{1}{τ} (  sup_{z∈W (y,r+τ/(1−r))} f(z) − f(x))  (3.17)

for all τ ∈]0, (1 − r)^{2}[.

Proof: Due to the subgradient inequality and Lemma 3.1.2,

∥p∥_{x,F}^{∘} = sup_{h:∥h∥_{x,F} ≤1}⟨ p, h ⟩
= \frac{1}{τ} sup_{h:∥h∥_{x,F} ≤τ}⟨ p, h ⟩
= \frac{1}{τ} sup_{z∈W (x,τ)}⟨ p, z − x ⟩
≤ \frac{1}{τ} sup_{z∈W (x,τ)}(f(z) − f(x))
≤ \frac{1}{τ} (  sup_{z∈W (y,r+τ/(1−r))} f(z) − f(x))

□

One can easily get rid of f(x) in (3.17) by using

f(x) ≥ inf_{z∈W (y,r)} f(z),  (3.18)

see again the proof of Lemma 3.1.2.

Numerous simple bounds for ∥p∥_{x,F}^{∘} can be derived from this lemma. For example, using τ := (1 − r)^{2}/2 and r < 1/3 leads to

∥p∥_{x,F}^{∘} ≤ \frac{2}{(1 − r)^{2}} ( sup_{z∈W (y,(1+r)/2)} f(z) − f(x))
< \frac{9}{2} ( sup_{z∈W (y,2/3)} f(z) − f(x)) ,  (3.19)

while τ := (1 − r)(1/2 − r) leads to

∥p∥_{x,F}^{∘} ≤ \frac{1}{(1 − r)(1/2 − r)} ( sup_{z∈W (y,1/2)} f(z) − f(x))
< 9 ( sup_{z∈W (y,1/2)} f(z) − f(x)) ,

which, together with (3.18), leads to the bound for ∥p∥_{x,F}^{∘} used in [67].

We now turn our attention to the efficiency of Newton’s method applied to the function F_{t}.

Lemma 3.1.4 Let there be given a 5-tuple (λ_{1}, λ′_{1}, λ_{2}, λ_{3}, λ′_{3}) feasible for the system (3.8)–(3.12), let x ∈ int(G) be given such that

λ(F_{t}, x) ≤ λ_{3}  (3.20)

holds, and define y by

y := x − (∇_{x}^{2}F_{t}(x))^{−1} ∇_{x}F_{t}(x).

Then y ∈ int(G) and

λ(F_{t}, y) ≤ λ′_{3}.  (3.21)

Proof: Due to (3.20), λ_{3} < 2 − \sqrt{3}, and [67, Theorem 2.2.2], we have y ∈ int(G) as well as

λ(F_{t}, y) ≤ (\frac{λ(F_{t}, x)}{1 − λ(F_{t}, x)})^{2} .

Using (3.20) again, it follows that λ(F_{t}, y) ≤ λ_{3}^{2}/(1−λ_{3})^{2}. The result follows with (3.10).

□

As a consequence of this lemma, defining the sequence (x^{(i)})_{i∈IN} by the Newtonprocess (3.14) and the update rule (3.15) will lead to λ(F_{t(i)}, x^{(i−1)}) ≤ λ_{3} as well as λ(F_{t(i)}, x^{(i)}) ≤ λ′_{3} for all i ≥ 1, provided that the starting point generated during the first stage is used. More precisely, we need λ(F_{t(1)}, x^{(0)}) < λ_{3} and a starting value t(1) for the barrier parameter bounded by the right-hand side of the estimate in Lemma 3.1.1.

It can be shown [67] that the number of iterations i needed to achieve an absolute accuracy of

f(x^{(i)}) − min_{x∈G} f(x) ≤ ε,  (3.22)

0 < ε < 1/∥∇f(x^{(0)})∥_{x^{(0)},F}^{∘}, is bounded above by

1 + \frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (ln (2ϑ + \frac{ζ(λ′_{3})}{2}) + ln(1/ε) + ln(∥∇f(x^{(0)})∥_{x^{(0)},F}^{∘})

+ ln (\frac{9}{λ_{3} − λ_{2}} ))

(3.23)

where the auxiliary function ζ is defined with the help of the function ω by

ζ(λ) := (ω(λ))^{2} \frac{1 + ω(λ)}{1 − ω(λ)}.

We will call all x ∈ G satisfying (3.22) an ε-solution to the optimization problem at hand. (Note that (3.8) implies λ_{3} ≠ λ′_{3}.) The term ∥∇f(x^{(0)})∥_{x^{(0)},F}^{∘} in this estimate can be bounded using one of the various estimates given after Lemma 3.1.3. All of these upper bounds have the form sup_{z∈W} f(z) − f(x^{(0)}) with a suitable set W ⊆ G containing the minimum of F. Note also that f(x^{(0)}) can be replaced by using (3.18). As a consequence, a scaling of our absolute accuracy measure ε reflects a scaling of the term ∥∇f(x^{(0)})∥_{x^{(0)},F}^{∘}.

Similar to the first stage, for large ϑ the leading term \sqrt{ϑ} ln ϑ has the coefficient 1/(λ_{3} − λ′_{3}), and it makes sense to choose λ_{3} = λ_{1} and λ′_{3} = λ′_{1}, with λ_{1}, λ′_{1} as above. It makes then sense to choose λ_{2} = (λ_{1} + λ′_{1})/2 = 0.1586 . . ..

3.2 Solving Two Single-Criteria Programs with Linear Objective Functions

Let c, ĉ ∈ IR^{m} with c ≠ 0 ≠ ĉ and let δ := ∥c − ĉ∥_{2} > 0. Again, let G be a compact convex set with nonempty interior. Moreover, suppose that ∥x∥_{2} ≤ R holds for a number R > 0 and all x ∈ G. We are now interested in solving the two linear optimization problems

minimize ⟨ c, x ⟩  (3.24)
subject to x ∈ G

and

minimize ⟨ ĉ, x ⟩  (3.25)
subject to x ∈ G.

Note that we do not yet want to consider the bicriteria program consisting of the two objective functions ⟨ c, · ⟩ and ⟨ ĉ, · ⟩. This is the subject of Section 3.5. Right now, we consider two different single-criteria problems and we want to solve both of them. Again, note that the discussion in Chapter 2 and that at the start of this chapter has shown that we have to solve many single-criteria linear programs, not just two of them. We discuss here the case of exactly two for the sake of simplicity.

Since in general we can not solve an optimization problem exactly, we have to be content with an ε-solution to each of our problems in the sense of (3.22), where ε > 0 is our prespecified accuracy measure. Now let x ∈ G be an ε-solution to (3.24) and suppose that ∥c∥_{2} = ∥ĉ∥_{2} holds. Then a simple calculation shows that x is an (ε + 2δR)-solution to (3.25). As a consequence, for sufficiently small perturbations δ one may accept the approximative solution to the first problem as an approximative solution to the second one. The following discussion will lead to an efficient algorithm for δ of moderate size.

The two objective functions given by f(x) = ⟨ c, x ⟩ as well as f̂(x) = ⟨ ĉ, x ⟩ are linear and therefore β -compatible with β = 0 for any self-concordant barrier. Now let F be such a self-concordant barrier for G with self-concordancy parameter ϑ. Obviously, we can compute an approximation x^{(0)} to the analytic center of G in the sense of (3.13) first, before proceeding to approximate a solution to (3.24) or (3.25). To solve both problems, we might then start from x^{(0)} to execute the second stage for the first problem, and then restart from x^{(0)} to solve the second problem, again executing the second stage. This, however, means that we are ”recycling” or ”reusing” only the starting point to solve two different problems which are usually in close connection to each other, i. e.

∥c − ĉ∥_{2} = δ ≪ 1. Obviously, the O-constants in (3.23) get multiplied by a factor of 2, independently of the size of δ.

It will now be shown that we can recycle additional steps of the main stage, and how the number of reusable steps depends on δ.

Let there be given a 5-tuple (λ_{1}, λ′_{1}, λ_{2}, λ_{3}, λ′_{3}) feasible for the system (3.8)–(3.12). Choose a sixth parameter λ̂_{3} such that λ_{3} < λ̂_{3} holds and the 5-tuple (λ_{1}, λ′_{1}, λ_{2}, λ̂_{3}, λ′_{3}) is feasible for (3.8)–(3.12), too. As an example, choose λ_{1} = 0.193, λ′_{1} = (λ_{1}/(1 − λ_{1}))^{2} = 0.057 …, λ_{2} = 0.150, λ_{3} = 4/25 = 0.16, λ′_{3} = 1/9, and λ̂_{3} = 1/4. Define the parameterized functions F_{t}, F̂_{t} in the usual way, i. e.

F_{t}(x) := tf(x) + F (x),
F̂_{t}(x) := tf̂(x) + F (x).

When solving (3.25), we can simply replace f by f̂, i. e. F_{t} by F̂_{t}. Suppose that we try to solve (3.24) and (3.25) with the parameter tuple (λ_{1}, λ′_{1}, λ_{2}, λ_{3}, λ′_{3}), using x^{(0)} as a starting point. According to Lemma 3.1.4, inequality (3.20) is sufficient for (3.21) if updates according to Newton’s scheme are made, and we might therefore concentrate on (3.20). The idea is now to use the points x^{(i)}, computed as approximations to the central path of problem (3.24), as approximations to the central path of the problem (3.25). More precisely, the point x^{(i−1)} will be used as an approximation to a point on the central path of problem (3.25) as long as

λ(F̂_{t(i)}, x^{(i−1)}) < λ̂_{3}  (3.26)

holds.

To start with this scheme, we need a starting point w := x^{(0)} ∈ G and a barrier parameter value t(1) with λ(F_{t(1)}, w) ≤ λ_{3} as well as λ(F̂_{t(1)}, w) ≤ λ_{3}. According to Lemma 3.1.1, the value

t(1) := \frac{λ_{3} − λ(F, w)}{max\{∥∇f(w)∥_{w,F}^{∘} , ∥∇f̂(w)∥_{w,F}^{∘} \}} = \frac{λ_{3} − λ(F, w)}{max\{∥c∥_{w,F}^{∘} , ∥ĉ∥_{w,F}^{∘} \}}  (3.27)

can be used as a starting value for both optimization problems. Then, (3.26) holds for i = 1. Assuming appropriate scaling of the objective functions, by (3.19) we clearly have

t(1) > \frac{λ_{3} − λ_{2}}{9R max\{∥c∥_{2}, ∥ĉ∥_{2}\}} = O (\frac{1}{R}) .  (3.28)

Suppose now that (3.26) does not hold for an i > 1, but that we still have that x^{(i−2)} is close to the central path of the second problem, i. e. λ(F̂_{t(i−1)}, x^{(i−2)}) < λ̂_{3}. How is it now possible to obtain a point x̂^{(i−2)} for which

λ(F̂_{t(i−1)}, x̂^{(i−2)}) ≤ λ_{3}  (3.29)

holds without too much work? (Especially, without running through the whole Stage 2 of the standard path-following algorithm again?) Inequality (3.29) means that (3.20) holds with respect to x̂^{(i−2)}. With such a point, we can go on with Stage 2 for problem (3.25) in the following way. Define the points y_{i−2}^{(j)} ∈ IR^{m} by making Newton steps for the function F̂_{t(i−1)}, starting at the point x^{(i−2)}: set y_{i−2}^{(0)} := x^{(i−2)} and

y_{i−2}^{(j+1)} := y_{i−2}^{(j)} − (∇^{2}F̂_{t(i−1)}(y_{i−2}^{(j)}))^{−1}∇F̂_{t(i−1)}(y_{i−2}^{(j)}),  (3.30)

j = 0, 1, . . .. According to Lemma 3.1.4 (cmp. also Theorem 2.2.3 and especially the implication (2.2.39) from [67]), one has due to λ̂_{3} < 2 − \sqrt{3} that

λ(F̂_{t(i−1)}, y_{i−2}^{(j+1)}) ≤ (\frac{λ(F̂_{t(i−1)}, y_{i−2}^{(j)})}{1 − λ(F̂_{t(i−2)}, y_{i−2}^{(j)})})^{2} ≤ \frac{λ(F̂_{t(i−1)}, y_{i−2}^{(j)})}{2} < \frac{λ̂_{3}}{2}

(j = 0, 1, …) holds. Obviously, λ(F̂_{t(i−1)}, y_{i−2}^{(j)}) < λ̂_{3}/2^{j} for j = 0, 1, … follows. Therefore, if we make at least

k := ⌈\frac{ln λ̂_{3} − ln λ_{3}}{ln 2} ⌉  (3.31)

Newton steps in the scheme to compute the y_{i−2}^{(j)}, we arrive at (3.29) by defining x̂^{(i−2)} := y_{i−2}^{(k)}. Obviously, this is a worst-case estimate, it may be the case that one needs fewer iterations.

We might envision these additional Newton steps as corrector steps. (In the example with the numerical values of the λ-parameters given as above, we obtain k > 0.64 … as a sufficient condition for k, so we need just one corrector step and immediately use x̂^{(i−2)} := y_{i−2}^{(1)}.) After the corresponding corrector steps, Stage 2 of the interior-point algorithm can be used on problem (3.25) as usual. In this case, we have saved the computations needed for the first i − 1 iterations (i. e. we have saved i Newton steps), and we need an additional amount of work of at most k iterations. The number of additional operations k, however, is independent of the size of the perturbation δ, since the points already generated are still relatively close to the central path of the second problem in the sense that λ(F̂_{t(i−1)}, x^{(i−2)}) < λ̂_{3} holds.

But how long does (3.26) hold, i. e. how long are we allowed to skip iterations to solve (3.25) while iterating to solve (3.24)? According to (3.16), we have

{z ∈ IR^{m} | ⟨ (z − x), ∇^{2}F_{t}(x)(z − x) ⟩ ≤ 1} ⊆ G

for all t > 0 and x ∈ int(G), because F_{t} is self-concordant for all t > 0 (compare Proposition 3.2.2, Definition 3.1.1, and Definition 2.1.1 of [67]). Since G is bounded, ∇^{2}F_{t}(x) is positive definite. Suppose that there is given a number R > 0 such that ∥x∥_{2} ≤ R holds for all x ∈ G. Denote by μ_{1}(x, t), …, μ_{n}(x, t) the eigenvalues of ∇^{2}F_{t}(x). Then, the inequality

μ_{i}(x, t) ≥ \frac{1}{4R^{2}}  (3.32)

(i = 1, …, n) holds. Obviously, this bound is independent of x and t. All this holds also if one replaces F_{t} by F̂_{t}. But ∇^{2}F_{t}(x) = ∇^{2}F̂_{t}(x) = ∇^{2}F (x) and ∇F_{t}(x) = tc+∇F (x) as well as ∇F̂_{t}(x) = tĉ + ∇F (x). As a consequence, we have

λ(F̂_{t}, x))^{2} = ⟨∇F̂_{t}(x), (∇^{2}F̂_{t}(x))^{−1}∇F̂_{t}(x)⟩
= ⟨t(ĉ − c) + tc + ∇F (x), (∇^{2}F (x))^{−1} (t(ĉ − c) + tc + ∇F (x))⟩
= t^{2} ⟨ĉ − c, (∇^{2}F (x))^{−1}(ĉ − c)⟩ + 2t ⟨ĉ − c, (∇^{2}F (x))^{−1}∇F_{t}(x)⟩
+ ⟨∇F_{t}(x), (∇^{2}F_{t}(x))^{−1}∇F_{t}(x)⟩
= t^{2} ⟨ĉ − c, (∇^{2}F (x))^{−1}(ĉ − c)⟩ − 2t ⟨ĉ − c, Δx⟩ + (λ(F_{t}, x))^{2},

where Δx is the Newton step defined by ∇^{2}F (x)Δx = −∇F_{t}(x). Consider now the case x = x^{(i−1)} and t = t(i). Since x^{(i)}, x^{(i−1)} ∈ G (cmp. [67, Proposition 3.2.4 (i)]), the inequality ∥Δx∥_{2} = ∥x^{(i)} − x^{(i−1)}∥_{2} ≤ 2R immediately follows. (This is a global, rather conservative estimate. What is important in what follows is that this estimate is independent of i.) Moreover, the largest eigenvalue of (∇^{2}F_{t}^{(2)}(x))^{−1} is bounded above by 4R^{2} uniformly in x and t (see (3.32)), and λ(F_{t(i)}, x^{(i−1)}) ≤ λ_{3} holds for all i (see Lemma 3.1.4). Therefore, the chain of equations above leads to

λ(F̂_{t(i)}, x^{(i−1)}) ≤ (4δ^{2}R^{2}(t(i))^{2} + 4t(i)δR + λ_{3}^{2})^{1/2}

(recall δ = ∥c − ĉ∥_{2}), which means that λ(F̂_{t(i)}, x^{(i−1)}) < λ̂_{3} holds as long as

4δ^{2}R^{2}(t(i))^{2} + 4t(i)δR + λ_{3}^{2} < λ̂_{3}^{2}

holds. Replacing the inequality sign by an equality sign, we obtain a quadratic equation in t(i) with two real-valued solutions, only one of them greater than 0. This solution can be written as

s := \frac{(λ̂_{3}^{2} − λ_{3}^{2} + 1)^{1/2} − 1}{2δR} = O (\frac{1}{δR}) = O (\frac{t(1)}{δ}) .  (3.33)

The update rule for t(i) is given by (3.15), while t(1) is initialized as in (3.27). What is needed is t(i) < s. Solving for i yields

i < 1 + \frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (ln(max\{∥c∥_{w,F}^{∘} , ∥ĉ∥_{w,F}^{∘} \}) − ln R − ln δ
+ ln((λ̂_{3}^{2} − λ_{3}^{2} + 1)^{1/2} − 1) − ln(λ_{3} − λ_{2}) − ln 2) .

Note that in this estimate, the term ln ϑ does not appear. Define

K := max{∥c∥_{w,F}^{∘} , ∥ĉ∥_{w,F}^{∘} }

and abbreviate the constant terms by

τ = τ(λ_{2}, λ_{3}, λ̂_{3}) := ln((λ̂_{3}^{2} − λ_{3}^{2} + 1)^{1/2} − 1) − ln(λ_{3} − λ_{2}) − ln 2.  (3.34)

It is clear that we does not need to take a look at problem (3.25) for the first

1 + [\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}}(ln K − ln R − ln δ + τ)]  (3.35)

iterations of the primal interior-point method applied to solve (3.24). Only after that we need to branch, i. e. take a look at two separate problems. (If the number above is less than zero, it is understood that branching is immediately necessary.) Note that this is a worstcase estimate, i. e. it may be the case that one can make more iterations before switching into a ”parallel” process by branching. Moreover, if δ is replaced by, e. g., δ/2, the number of iterations we are allowed to execute before branching increases by O(\sqrt{ϑ}). We will return to this observation in Section 3.5. (Using a larger value for K when defining t(1), like the one given in estimate (3.28), will increase the number of iterations in which (3.26) holds.)

Recall that ε > 0 is the prespecified solution accuracy that we want to obtain and suppose that

δ ≥ ε^{2} \frac{(λ_{3} − λ_{2})^{2}e^{τ}}{81(2ϑ + ζ(λ′_{3})/2)^{2}KR} = O (\frac{ε^{2}}{ϑ^{2}KR})  (3.36)

holds. Taking the result obtained above together and considering (3.23), we see that

2 + 2\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (ln (2ϑ + \frac{ζ(λ′_{3})}{2}) + ln(1/ε) + ln K + ln (\frac{9}{λ_{3} − λ_{2}} ))

−\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}}(ln K − ln R − ln δ + τ) − 1 + k

= \frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (2 ln (2ϑ + \frac{ζ(λ′_{3})}{2}) + 2 ln(1/ε) + ln K + ln R + ln δ (3.37)

+2 ln (\frac{9}{λ_{3} − λ_{2}} ) − τ ) + k + 1

Newton-steps are needed instead of

2 + 2\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (ln (2ϑ + \frac{ζ(λ′_{3}}{2} ) + ln(1/ε) + ln K + ln (\frac{9}{λ_{3} − λ_{2}} ))

Newton-steps in Stage 2 to achieve an accuracy of ε for both problems. The lower bound (3.36) just ensures that we have not obtained ε-accuracy before branching and there is still some work to do, and that therefore the estimate above on the number of steps is a nonnegative number. Instead of introducing (3.36), one might replace (3.37) by an estimate of the form max{0, …}, where the dots stand for (3.37).

If δ < 1, there are savings in computation time proportional to

−\sqrt{ϑ}(ln R + ln δ) = \sqrt{ϑ} ln(1/Rδ)

iterations. Obviously, substantial savings will take place if δ < 1/R. Again, note that this estimate is rather conservative, mainly due to the fact that the inequality ∥x^{(i)}−x^{(i−1)}∥_{2} ≤ 2R has been used to derive (3.33). Note also that Stage I is still needed, though.

The following algorithm summarizes the main points of the scheme, and the next theorem restates the results found by the analysis above.

Algorithm 3.2.1

1. Choose (λ_{1}, λ′_{1}, λ_{2}, λ_{3}, λ′_{3}) ∈ IR^{5} feasible for the system (3.8)–(3.12).

Choose λ̂_{3} ∈ IR such that λ_{3} < λ̂_{3} holds and such that (λ_{1}, λ′_{1}, λ_{2}, λ̂_{3}, λ′_{3}) is feasible for (3.8)–(3.12), too.

2. Compute a point x^{(0)} ∈ int(G) such that (3.13) holds, i. e. λ(F, x^{(0)}) < λ_{2}. (This is Stage I and can be done along the lines of [67, Section 3.2.3].)

3. (Initialization of barrier parameter, compare (3.27).)

t(1) := \frac{λ_{3} − λ(F, x^{(0)})}{max\{∥c∥_{x^{(0)},F}^{∘} , ∥ĉ∥_{x^{(0)},F}^{∘} \}}

4. i := 1

5. (Solve the first problem.) WHILE NOT stopping criterion for problem (3.24) fulfilled DO

(a) (Do a Newton-step, compare (3.14).)

x^{(i)} := x^{(i−1)} − (∇_{x}^{2}F_{t(i)}(x^{(i−1)}))^{−1} ∇_{x}F_{t(i)}(x^{(i−1)})

(b) (Update barrier parameter compare (3.15).)

t(i + 1) := t(i) exp (\frac{λ_{3} − λ′_{3}}{λ_{3} + \sqrt{ϑ}})

(c) i := i + 1

6. (Compute index of branching step.)

i := max{i ≥ 0 | λ(F̂_{t(i)}, x^{(i−1)}) < λ̂_{3}}

7. y_{i−1}^{(0)} := x^{(i−1)}, j := 1

8. (Branching) WHILE λ(F̂_{t(i)}, y_{i−1}^{(j−1)}) ≥ λ_{3} DO

(a) (Do a corrector Newton-step, compare (3.30).) y_{i−1}^{(j)} := y_{i−1}^{(j−1)} − (∇_{x}^{2}F̂_{t(i)}(y_{i−1}^{(j−1)}))^{−1} ∇_{x}F̂_{t(i)}(y_{i−1}^{(j−1)})

(b) j := j + 1

9. x̂^{(i−1)} := y_{i−1}^{(j−1)}

10. (Solve the second problem.) WHILE NOT stopping criterion for problem (3.25) fulfilled DO

(a) (Do a Newton-step, compare (3.14).)

x̂^{(i)} := x̂^{(i−1)} − (∇_{x}^{2}F̂_{t(i)}(x̂^{(i−1)}))^{−1} ∇_{x}F̂_{t(i)}(x̂^{(i−1)})

(b) (Update barrier parameter, compare (3.15).)

t(i + 1) := t(i) exp (\frac{λ_{3} − λ′_{3}}{λ_{3} + \sqrt{ϑ}})

(c) i := i + 1 Theorem 3.2.2 Let G ⊆ IR^{m} be convex and compact with nonempty interior and let F be a self-concordant barrier for G with self-concordancy parameter ϑ. Moreover, let ∥x∥_{2} ≤ R hold for all x ∈ G for a number R > 0. Let c, ĉ ∈ IR^{m} with c ≠ 0 ≠ ĉ and let there be given an accuracy measure ε > 0. With the λ-parameters chosen in Step 1 of Algorithm 3.2.1 and x^{(0)} computed as in Step 2, define

τ = τ(λ_{2}, λ_{3}, λ̂_{3}) := ln((λ̂_{3}^{2} − λ_{3}^{2} + 1)^{1/2} − 1) − ln(λ_{3} − λ_{2}) − ln 2

as well as

K := max{∥c∥_{x^{(0)},F}^{∘} , ∥ĉ∥_{x^{(0)},F}^{∘} }.

Suppose that

∥c − ĉ∥_{2} ≥ ε^{2} \frac{(λ_{3} − λ_{2})((λ̂_{3}^{2} − λ_{3}^{2} + 1)^{1/2} − 1)}{162(2ϑ + ζ(λ′_{3})/2)^{2}KR}

as well as 0 < εK < 1 holds. Then, the total number of Newton-steps in Steps 5–10 of Algorithm 3.2.1 needed to achieve ε-accuracy for both optimization problems is bounded above by

\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (2 ln (2ϑ + \frac{ζ(λ′_{3})}{2}) + 2 ln(1/ε) + ln K + ln R + ln ∥c − ĉ∥_{2}

+2 ln (\frac{9}{λ_{3} − λ_{2}} ) − τ ) + ⌈\frac{ln λ̂_{3} − ln λ_{3}}{ln 2} ⌉ + 1.

Note again that the lower bound on the perturbation size δ = ∥c − ĉ∥_{2} in the theorem above was introduced just for convenience. For smaller δ it might be the case that the WHILE-loop in Step 10 is never executed.

We would like to remark here that the exact computation of the O-constants of the runtime behaviour of Algorithm 3.2.1, as done above, is, unfortunately, not ”computational gymnastics” [76], although we wished that this would be the case. It turns out that we need precise knowledge of these constants to equilibrate the amount of work necessary for solving a large number of optimization problems (and not just two of them) between different parts of a corresponding algorithm. This is the subject of Section 3.5.

Let us also stress that the number R enters the run-time estimate of the theorem, while it is not used in algorithm. No knowledge on G is needed other than boundedness and the fact we know a self-concordant barrier F for G.

The approach presented above uses only primal information and makes only few assumptions about the feasible set G. (Most importantly, that G is compact and that we know a self-concordant barrier for it.) If primal-dual algorithms are considered, care has to be taken with respect to the (strict) feasibility of the iterates considered. Indeed, if a primal objective function is perturbed, the perturbation shows up in the dual constraints, and it is not a priori clear how dual feasibility might be guaranteed. Up to now, only primal-dual schemes for problems with special structure, most notably with linear constraints and linear objective functions [97] or convex-quadratic objective functions [35] have been considered.

3.3 Solving Two Single-Criteria Programs with Quadratic Objective Functions

Suppose now that we have given c, ĉ ∈ IR^{m} as well as Q, Q̂ ∈ IR^{m×m} positive definite. Define δ_{c} := ∥c − ĉ∥_{2} and δ_{Q} := ∥Q − Q̂∥. Again, let G be a compact convex set with nonempty interior such that there is given a number R > 0 with ∥x∥_{2} ≤ R for all x ∈ G.

We are interested in solving the two quadratic optimization problems

minimize ⟨ c, x ⟩ + \frac{1}{2}⟨ x, Qx ⟩  (3.38)
subject to x ∈ G

and

minimize ⟨ ĉ, x ⟩ + \frac{1}{2}⟨ x, Q̂x ⟩  (3.39)
subject to x ∈ G

The two objective functions f(x) := ⟨ c, x ⟩ + (1/2)⟨ x, Qx ⟩ and f̂(x) := ⟨ ĉ, x ⟩ + (1/2)⟨ x, Q̂x ⟩ are β -compatible with β = 0 for any self-concordant barrier. Let F be such a self-concordant barrier for G with self-concordancy parameter ϑ. We will make use of the following simple observation.

Lemma 3.3.1 Let A, ΔA ∈ IR^{m×m} and b, x, Δx ∈ IR^{m}. Suppose that Ax = b as well as (A + ΔA)(x + Δx) = b holds. Moreover, let A as well as A + ΔA be regular. It then follows that Δx = −(A + ΔA)^{−1}ΔAx holds.

Proof: evident. □

Using the same strategy and notation as in the last section, we see that we need an estimate for t such that λ(F̂_{t}, x) < λ̂_{3} holds for all x ∈ int(G) with λ(F_{t}, x) < λ_{3}. By definition, we have

(λ(F̂_{t}, x))^{2} = ⟨∇F̂_{t}(x), (∇^{2}F̂_{t}(x))^{−1}∇F̂_{t}(x)⟩
= ⟨t∇f̂(x) + ∇F (x), (∇^{2}F̂_{t}(x))^{−1}(t∇f̂(x) + ∇F (x))⟩ ,

and due to t∇f̂(x) + ∇F (x) = t(∇f̂(x) − ∇f(x)) + ∇F_{t}(x)) we arrive at

(λ(F̂_{t}, x))^{2}
= t^{2} ⟨∇f̂(x) − ∇f(x), (∇^{2}F̂_{t}(x))^{−1}(∇f̂(x) − ∇f(x))⟩
+2t ⟨∇f̂(x) − ∇f(x), (∇^{2}F̂_{t}(x))^{−1}∇F_{t}(x)⟩
+ ⟨∇F_{t}(x), (∇^{2}F̂_{t}(x))^{−1}∇F_{t}(x)⟩ .

Since F̂_{t} is self-concordant, the eigenvalues of ∇^{2}F̂_{t}(x) are bounded below by 1/4R^{2}. The estimate

⟨∇f̂(x) − ∇f(x), (∇^{2}F̂_{t}(x))^{−1}(∇f̂(x) − ∇f(x))⟩
≤ 4R^{2}∥∇f̂(x) − ∇f(x)∥_{2}^{2}  (3.40)

immediately follows. Suppose that the sequence (x^{(i)})_{i∈IN} (x^{(i)} ∈ G) is generated by the standard primal path-following algorithm applied to (3.38) with the barrier F. Write ∇^{2}F̂_{t}(x) = t∇^{2}f̂(x)+∇^{2}F (x) = t(∇^{2}f̂(x)−∇^{2}f(x))+∇^{2}F_{t}(x) and use Lemma 3.3.1 with A = ∇^{2}F_{t}(x), ΔA = t(∇^{2}f̂(x) − ∇^{2}f(x)) as well as b = ∇F_{t}(x). Together with x^{(i)} − x^{(i+1)} = (∇^{2}F_{t}(x))^{−1}∇F_{t}(x), this leads to

(∇^{2}F̂_{t}(x))^{−1}∇F_{t}(x)
= (I − t(∇^{2}F̂_{t}(x))^{−1}(∇^{2}f̂(x) − ∇^{2}f(x)))(x^{(i)} − x^{(i+1)})  (3.41)

and therefore

⟨∇f̂(x) − ∇f(x), (∇^{2}F̂_{t}(x))^{−1}∇F_{t}(x)⟩
≤ ∥∇f̂(x) − ∇f(x)∥_{2}∥x^{(i)} − x^{(i+1)}∥_{2}(1 + 4R^{2}t∥∇^{2}f̂(x) − ∇^{2}f(x)∥)
≤ 2R∥∇f̂(x) − ∇f(x)∥_{2}(1 + 4R^{2}t∥∇^{2}f̂(x) − ∇^{2}f(x)∥).  (3.42)

Moreover, using (3.41) twice again leads to

⟨∇F_{t}(x), (∇^{2}F̂_{t}(x))^{−1}∇F_{t}(x)⟩

= (λ(F_{t}, x))^{2} − t ⟨∇F_{t}(x), (∇^{2}F̂_{t}(x))^{−1}(∇^{2}f̂(x) − ∇^{2}f(x))(∇^{2}F_{t}(x))^{−1}∇F_{t}(x)⟩
< λ_{3}^{2} + t∥x^{(i)} − x^{(i+1)}∥_{2}∥(∇^{2}F̂_{t}(x))^{−1}∇F_{t}(x)∥_{2}∥∇^{2}f̂(x) − ∇^{2}f(x)∥
≤ t∥x^{(i)} − x^{(i+1)}∥_{2}^{2}∥I − t(∇^{2}F̂_{t}(x))^{−1}(∇^{2}f̂(x) − ∇^{2}f(x))∥∥∇^{2}f̂(x) − ∇^{2}f(x)∥
+λ_{3}^{2}

which results in

⟨∇F_{t}(x), (∇^{2}F̂_{t}(x))^{−1}∇F_{t}(x)⟩
< λ_{3}^{2} + 4R^{2}t(1 + 4R^{2}t∥∇^{2}f̂(x) − ∇^{2}f(x)∥)∥∇^{2}f̂(x) − ∇^{2}f(x)∥. (3.43)

By noting ∥∇f̂(x) − ∇f(x)∥_{2} ≤ δ_{c} + δ_{Q}R, ∥∇^{2}f̂(x) − ∇^{2}f(x)∥ = δ_{Q} as well as (3.40), (3.42), and (3.43), we arrive at

(λ(F̂_{t}, x))^{2} < 4R^{2}(δ_{c} + δ_{Q}R)^{2}t^{2} + 4Rt(δ_{c} + δ_{Q}R)(1 + 4R^{2}δ_{Q}t)
+λ_{3}^{2} + 4R^{2}δ_{Q}t(1 + 4R^{2}δ_{Q}t)
= 4R^{2}(δ_{c}^{2} + 6Rδ_{c}δ_{Q} + 9R^{2}δ_{Q}^{2})t^{2} + 4R(2Rδ_{Q} + δ_{c})t + λ_{3}^{2}
= 4R^{2}(δ_{c} + 3Rδ_{Q})^{2}t^{2} + 4R(δ_{c} + 2Rδ_{Q})t + λ_{3}^{2}.

Accordingly, λ(F̂_{t}, x) < λ̂_{3} holds as long as

4R^{2}(δ_{c} + 3Rδ_{Q})^{2}t^{2} + 4R(δ_{c} + 2Rδ_{Q})t + λ_{3}^{2} < λ̂_{3}^{2}.

Solving for t yields

t < \frac{δ_{c} + 2Rδ_{Q}}{2R(δ_{c} + 3Rδ_{Q})^{2}} ((1 + λ̂_{3}^{2} − λ_{3}^{2})^{1/2} − 1) .  (3.44)

It is easy to see that (δ_{c} + 2Rδ_{Q})/(δ_{c} + 3Rδ_{Q})^{2} ≥ 1/(δ_{c} + (9/4)Rδ_{Q}) holds. Therefore,

t < \frac{(1 + λ̂_{3}^{2} − λ_{3}^{2})^{1/2} − 1}{2R(δ_{c} + (9/2)Rδ_{Q})} =: s

is sufficient for λ(F̂_{t}, x) < λ̂_{3}. Clearly, this estimate is a generalization of (3.33). In virtually the same way as in the last section, we conclude that

1 + [\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}}(ln K − ln R − ln(δ_{c} + (9/2)Rδ_{Q}) + τ)]  (3.45)

iterations of the primal interior-point method applied to solve problem (3.38) can be made before considering problem (3.39). Here,

K := max{∥c + Qx^{(0)}∥_{x^{(0)},F}^{∘} , ∥ĉ + Q̂x^{(0)}∥_{x^{(0)},F}^{∘} }.

Likewise,

\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (2 ln (2ϑ + \frac{ζ(λ′_{3})}{2}) + 2 ln(1/ε) + ln K + ln R

+ ln (δ_{c} + \frac{9}{2}Rδ_{Q}) + 2 ln (\frac{9}{λ_{3} − λ_{2}} ) − τ ) + k + 1

Newton-steps are needed to solve both problems to ε-optimality. Note that with δ := max{δ_{c}, δ_{Q}} we have that

ln R + ln(δ_{c} + (9/2)Rδ_{Q}) = O (ln δ + ln R) ,

so, from the point of view of complexity theory, the situation is not worse than in the linear case. On the other hand, it has to be noted that the computation of the quantity λ(F̂_{t}, x) needs as many operations as a Newton step for the function F̂_{t}. This is in stark contrast to the linear case, where ∇_{x}^{2}F̂_{t}(x) = ∇_{x}^{2}F_{t}(x) holds and it is therefore possible to reuse the factorization of the corresponding matrices. As such, as long as we do not want to fall back on heuristics to estimate the number of steps we are allowed to recycle, we have that at the present moment only the estimate (3.44) is available for saving unnecessary work during the course of solving two quadratic programs at once. For some heuristics useful in a primal-dual framework, see [49].

3.4 Solving Two Single-Criteria Programs with General Objective Functions

It is easy to see that the ideas from the last section carry over to general case of two smooth convex objective functions f, f̂ ∈ C^{3}(G). Indeed, if uniform bounds δ_{1}, δ_{2} of the form ∥∇f̂(x) − ∇f(x)∥_{2} ≤ δ_{1}, ∥∇^{2}f̂(x) − ∇^{2}f(x)∥ ≤ δ_{2} for all x ∈ G are known, the estimates (3.40), (3.42), and (3.43) immediately lead to

(λ(F̂_{t}, x))^{2} ≤ 4R^{2}δ_{1}^{2}t^{2} + 4Rδ_{1}t(1 + 4R^{2}δ_{2}t) + λ_{3}^{2} + 4R^{2}δ_{2}t(1 + 4R^{2}δ_{2}t)
= 4R^{2}(δ_{1} + 2Rδ_{2})^{2}t^{2} + 4R(δ_{1} + Rδ_{2})t + λ_{3}^{2}

and an analysis similar to the one in the last two sections will lead to similar results, provided that the parameter needed in defining β -compatibility is used in the interior-point method in the standard way (cmp. (3.15) and [67, Chapter 3]). This case, however, is of minor interest when discussing the use of interior-point methods for multicriteria problems, since the results in Chapter 2 have shown that parameterized linear or parameterized quadratic problems are sufficient for approximating the set of efficient points.

3.5 Solving Many Convex Programs ”Simultaneously”

Suppose now that two cost vectors c_{1}, c_{2} ∈ IR^{m} are given and that we want to solve ”all” the problems with set of feasible points G and with cost vectors c(μ) := μc_{1} + (1 − μ)c_{2} (μ ∈]0, 1[). More precisely, we are in search for an approximation of the set

⋃_{μ∈]0,1[} arg min{⟨ c(μ), x ⟩ | x ∈ G}  (3.46)

by a discrete set of points. This is exactly the situation in standard bicriteria optimization in with as a cone K the dual of the conic hull of the vectors c_{1} and c_{2} is used.

In a standard approach, our task amounts to discretizing the parameter interval [ 0, 1] by a (prespecified) number of parameters μ_{i} such that ∥c(μ_{i}) − c(μ_{i+1})∥_{2} ≤ δ holds with a prespecified accuracy measure δ ∈]0, 1[. Under the assumption ∥c_{1}∥_{2}, ∥c_{2}∥_{2} ≤ 1 (made without loss of generality) we might use any discretization μ_{i} of [0, 1] for which |μ_{i} − μ_{i+1}| ≤ δ/2 holds. With this, we basically have to solve 2/δ different standard optimization problems with the same set of feasible points G. To solve these optimization problems numerically, an accuracy ε > 0, used in a stopping criterion for the numerical algorithm used for the standard optimization problems, should also be given.

In what follows, a numerical scheme is described that uses the ideas from the last section to approximate the set (3.46). First, a short run-down of the general idea is presented, before the actual algorithm is described and its complexity is analyzed.

3.5.1 The General Idea

Using the results from the first section of this chapter, we start by defining δ_{0} := 1/2 and μ_{0} := 1/2. In this way, instead of solving two problems ”simultaneously”, we execute Newton iterations for all problems with cost vectors c(μ) for which |μ − μ_{0}| ≤ δ_{0} holds (i. e. all problems we have to consider) up to the time when we have to ”branch”. (Not counting the ubiquitous Stage 1, this takes

\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (ln K − ln R + ln 2 + τ)

Newton iterations, where

K := max{∥c(μ)∥_{x^{(0)},F} | μ ∈ [0, 1]} ≤ 2R.)

After that, a ”branch” is needed to proceed. Define δ_{1} :== δ_{0}/2 = 1/4, μ_{1}^{(1)} := 1/4, and μ_{1}^{(2)} := 3/4. Now, the two problems corresponding to the parameter values μ_{1}^{(1)}, μ_{1}^{(2)} are considered. In order to reach the neighborhoods of their central paths, at most 2k corrector steps, where k is given as in (3.31), are needed. After these corrector steps, Newton iterations are executed for all problems with parameters μ for which |μ−μ_{1}^{(1)}| ≤ δ_{1} holds until we have to branch and after that, Newton iterations are executed for all problems with parameters μ for which |μ − μ_{1}^{(2)}| ≤ δ_{1} until branching is needed. (Again, this means that all possible problems are considered. Moreover, by the estimate (3.35) on the number of steps saved, it is clear that

\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} ln 2

Newton-steps are made in each branch before the process stops and we have to branch again.)

Then, define δ_{2} := 1/8, μ_{2}^{(1)} = 1/8, μ_{2}^{(2)} = 3/8, μ_{2}^{(3)} = 5/8, μ_{2}^{(4)} = 7/8, make 4k corrector steps, and repeat the scheme.

This scheme of branching and refining will be used until we arrive at a stage N with δ_{N} = 1/2^{N+1} ≤ δ. Then, 1/δ_{N} ≤ 1/δ different problems are considered and we can now solve them to ε-optimality in the standard way. The different parameters μ_{i}, i = 1, …, 2^{N+1}, considered after the final branching step define a δ-covering of the unit interval [0, 1], i. e. μ_{i} ∈ [0, 1] for all i and for each μ ∈ [0, 1] there exists an i with |μ − μ_{i}| < δ. (In fact, if equality is allowed in the last inequality, the points μ_{i} (i = 1, …, 2^{N+1}) define a δ_{N} -covering.)

3.5.2 The Algorithm

The scheme presented in the last section can be described more formally by the following algorithm. Each problem is identified with its parameter μ ∈ [0, 1], and the central path of each of these problems is followed by the standard path-following method. With respect to bookkeeping, for each problem only one vector of variables x_{μ} ∈ G and one barrier parameter t_{μ} > 0 is needed. In what follows, the abbreviations

c(μ) := μc_{1} + (1 − μ)c_{2},
f_{μ}(x) := ⟨ c(μ), x ⟩,
F_{t}^{μ}(x) := tf_{μ}(x) + F (x)

will be used.

Algorithm 3.5.1 (Approximating the Efficient Set of a Bicriteria Program)

1. Choose (λ_{1}, λ′_{1}, λ_{2}, λ_{3}, λ′_{3}) ∈ IR^{5} feasible for the system (3.8)–(3.12).

Choose λ̂_{3} ∈ IR such that λ_{3} < λ̂_{3} holds and such that (λ_{1}, λ′_{1}, λ_{2}, λ̂_{3}, λ′_{3}) is feasible for (3.8)–(3.12), too.

2. Choose a discretization size 0 < δ < 1.

3. Compute a point x ∈ int(G) such that (3.13) holds, i. e. λ(F, x) < λ_{2}. (This is Stage I and can be done along the lines of [67, Section 3.2.3].)

4. (Initialization of the barrier parameter, compare (3.27) and Lemma 3.1.1.)

t := \frac{λ_{3} − λ(F, x)}{max\{∥c(μ)∥_{x,F}^{∘} | μ ∈ [0, 1]\}}

5. δ(0) := 1/2, μ := 1/2, S(0) := {μ}, j := 0, x_{μ} := x, t_{μ} := t

6. WHILE δ(j) > δ DO

(a) FORALL μ ∈ S(j) DO

i. WHILE λ(F_{t_{μ}}^{μ}, x_{μ}) ≥ λ_{3} DO (corrector steps for F_{t_{μ}}^{μ})

x_{μ} := x_{μ} − (∇_{x}^{2}F_{t_{μ}}^{μ}(x_{μ}))^{−1} ∇_{x}F_{t_{μ}}^{μ}(x_{μ})

ii. WHILE λ(F_{t_{μ}}^{μ+δ(j)}, x_{μ}) < λ̂_{3} AND λ(F_{t_{μ}}^{μ−δ(j)}, x_{μ}) < λ̂_{3} AND NOT stopping criterion for problem with parameter μ fulfilled DO (Path following for the problem with parameter μ until branching is necessary.)

A. x_{μ}^{old} := x_{μ}, t_{μ}^{old} := t_{μ} B. (Newton step)

x_{μ} := x_{μ} − (∇_{x}^{2}F_{t_{μ}}^{μ}(x_{μ}))^{−1} ∇_{x}F_{t_{μ}}^{μ}(x_{μ})

C. (Update barrier parameter.)

t_{μ} := t_{μ} exp (\frac{λ_{3} − λ′_{3}}{λ_{3} + \sqrt{ϑ}})

iii. x_{μ} := x_{μ}^{old}, t_{μ} := t_{μ}^{old}

(b) (Now x_{μ} represents the point along the central path for the problem with parameter μ as far away from the analytic center x as possible.) j := j + 1, δ(j) := δ(j − 1)/2

(c) FORALL μ ∈ S(j) DO (branch)

x_{μ−δ(j)} := x_{μ}, x_{μ+δ(j)} := x_{μ}, t_{μ−δ(j)} := t_{μ}, t_{μ+δ(j)} := t_{μ}

(d) (refinement) S(j) := {μ + δ(j), μ − δ(j) | μ ∈ S(j − 1)}

7. (Solve all problems in S(j) to a specified accuracy.)

FORALL μ ∈ S(j) DO

(a) (Corrector steps for the last refinement.) WHILE λ(F_{t_{μ}}^{μ}, x_{μ}) ≥ λ_{3} DO

x_{μ} := x_{μ} − (∇_{x}^{2}F_{t_{μ}}^{μ}(x_{μ}))^{−1} ∇_{x}F_{t_{μ}}^{μ}(x_{μ})

(b) WHILE NOT stopping criterion for problem with parameter μ fulfilled DO

i. (Newton step)

x_{μ} := x_{μ} − (∇_{x}^{2}F_{t_{μ}}^{μ}(x_{μ}))^{−1} ∇_{x}F_{t_{μ}}^{μ}(x_{μ})

ii. (Update barrier parameter.)

t_{μ} := t_{μ} exp (\frac{λ_{3} − λ′_{3}}{λ_{3} + \sqrt{ϑ}})

The discrete sets S(j) ⊂ [0, 1] represent the set of parameters which are considered as problem parameters in the path-following scheme. Note that, due to step 6 (b) and (d), we have S(j − 1) ∩ S(j) = ∅.

3.5.3 Analysis

As usual, define

K := max{∥c(μ)∥_{x^{(0)},F}^{∘} | μ ∈ [0, 1]},

where x^{(0)} ∈ int(G) is the starting point computed in Step 3 of Algorithm 3.5.1. (Proper scaling of c_{1} and c_{2} ensures K ≤ 9R.) At the j th step of the WHILE-loop of Step 6 of the algorithm (j = 0, 1, …), the set S(j) is a δ(j)-covering of the parameter set [0, 1], and we have δ(j) = 1/2^{j+1} as well as |S(j)| = 2^{j} Moreover, for j > 0 we have |S(j)| = 1/δ(j − 1). The number of times the WHILE-loop of Step 6 is executed is given by

N := ⌈\frac{ln(1/δ)}{ln 2}⌉ − 1.

The total number of all elements in all sets S(j) considered, i. e. the number of branches made is ∑_{j=0}^{N−1} |S(j)| = ∑_{j=0}^{N−1} 2^{j} = 2^{N}−1 < 1/δ−1. As a consequence, the total number of corrector steps in Step 6 (a) (i) is at most k(1/δ − 1), where k is the constant given by (3.31). An additional number of at most kS(N) ≤ k/δ corrector steps is executed in Step 7 (a). Therefore, the total number of corrector steps executed by the algorithm is bounded above by 2k/δ − k.

Assume now in what follows that with respect to the amount of work done the worst case occurs, i. e. the stopping criterion of the WHILE-loop of Step 6 (a) ii is fulfilled as soon as possible. Then, the smallest possible number of Newton steps has been executed since the last branch. According to the analysis in Section 3.2 (see especially (3.35)), there will be

I := 2 + [\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (ln K − ln R + ln 2 + τ)]

Newton steps made before the first branch occurs. Now let j ∈ {1, …, N − 1}. Since δ(j) was replaced by δ(j − 1)/2 in previous exectution of the WHILE-loop of Step 6, it is clear by the estimate (3.35) on the number of steps saved that at most

1 + ⌈\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} ln 2⌉

Newton-steps are executed for each problem considered before branching occurs again. (With the asymptotically optimal (for large ϑ) values for λ_{3} and λ′_{3} given on p. 34, the bracketed number is approximatively 4.922(0.229 + \sqrt{ϑ}) ≥ 6.050.) This number has to be multiplied with the number of problems considered during this run through the WHILEloop of Step 6 of Algorithm 3.5.1, which is 2^{j}. Using the same estimate as above, the total number of Newton steps made in the WHILE-loop of Step 6 after the first branch is bounded above by

(\frac{1}{δ} − 1) (⌈ln 2\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}}⌉ + 1) .

Therefore, the total number of Newton steps executed in Step 6 and 7 (a) is at most

I + (\frac{1}{δ} − 1) (⌈ln 2\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}}⌉ + 1) + \frac{2k}{δ} − k

≤ 3 + \frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (ln K − ln R + ln 2 + τ)

+ (\frac{1}{δ} − 1) (ln 2\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} + 2) + \frac{2k}{δ} − k

≤ 1 − k + \frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (\frac{ln 2}{δ} + ln K − ln R + τ ) + \frac{2k + 2}{δ}.  (3.47)

But after the execution of Step 7 (a) we have for each of the S(N) = 2^{N} < 1/δ different single-criterion problems considered found a point close to its central path as if we had made at least

I − 1 + N ⌊ln 2\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}}⌋

≥ I − 1 + (\frac{ln(1/δ)}{ln 2} − 1) (ln 2\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} − 1)
≥ 1 + \frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (ln K − ln R + ln 2 + τ) + (\frac{ln(1/δ)}{ln 2} − 1) (ln 2\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} − 1)
= 2 + \frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (ln K − ln R + τ + ln(1/δ)) − \frac{ln(1/δ)}{ln 2}

Newton steps in Stage 2 of the standard primal path-following algorithm. Now suppose that an accuracy measure ε > 0 is specified for all problems considered. By assuming

δ ≥ ε\frac{(λ_{3} − λ_{2})e^{τ}}{9R(2ϑ + ζ(λ′_{3})/2)}

and using the estimate above as well as (3.23), we see that we need at most

\frac{ln(1/δ)}{ln 2}  (3.48)

+\frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}} (ln (2ϑ + \frac{ζ(λ′_{3})}{2}) + ln(1/ε) + ln (\frac{9}{λ_{3} − λ_{2}} ) + ln R − τ + ln δ)

further Newton steps for each of the problems considered to solve it to ε-optimality, which is exactly what is done in Step 7 (b) of Algorithm 3.5.1. The lower bound on δ was introduced for the same reasons as in Section 3.2: it is convenient to assume that none of the problems considered is already solved to ε-accuracy in Step 6 of our algorithm. If this assumption is not satisfied, the algorithm will be even faster. Multiplying (3.48) by 1/δ and adding the result to the actual number of Newton steps already made (see (3.47)), we see that we need in total at most

\frac{1}{δ} \frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}}(ln ϑ+ln(1/ϵ)+ln R+ln δ+C_{1})+ \frac{2k}{δ} + \frac{ln(1/δ)}{δ ln 2} + \frac{λ_{3} + \sqrt{ϑ}}{λ_{3} − λ′_{3}}(ln K−ln R+C_{2})

(C_{1}, C_{2} suitable constants depending only on λ_{1}, λ′_{1}, λ_{2}, λ_{3}, λ′_{3}, λ̂_{3}) Newton iterations to solve all problems considered to ε-optimality. Note that, due to ϑ ≥ 1, (3.8), and (3.10), (λ_{3} + \sqrt{ϑ})/(λ_{3} − λ′_{3}) > 1/ ln 2 holds. Moreover, if the asymptotically optimal (for large ϑ) values for the λ-parameters are used, the (ln δ)/δ-term from above has approximatively the form

(7.102(\sqrt{ϑ} + 0.229) − 1.442)\frac{ln δ}{δ}

and is therefore negative as long as δ < 1. As a consequence, we can drop the term ln(1/δ)/δ ln 2 in the estimate above completely by re-adjusting the constant C_{1}.

We summarize our findings in the next theorem.

Theorem 3.5.2 (Complexity Estimate) Let G ⊆ IR^{m} be convex and compact with nonempty interior and let F be a self-concordant barrier for G with self-concordancy parameter ϑ. Moreover, let ∥x∥_{2} ≤ R hold for all x ∈ G and a number R > 0. Let c_{1}, c_{2} ∈ IR^{m} with c ≠ 0 ≠ c_{2} and let there be given an accuracy measure ε > 0. With x^{(0)} computed in Step 3 of Algorithm 3.5.1, define

K := max{∥μc_{1} + (1 − μ)c_{2}∥_{x^{(0)},F}^{∘} | μ ∈ [0, 1]},

and with the λ-parameters chosen as in Step 1 suppose that

δ ≥ ε\frac{(λ̂_{3}^{2} − λ_{3}^{2} + 1)^{1/2} − 1}{18R(2ϑ + ζ(λ′_{3})/2)}

as well as 0 < εK < 1 holds. Then, the total number of Newton-steps in Steps 6– 7 of Algorithm 3.2.1 needed to compute a δ-covering (μ_{i})_{i} of [0, 1], μ_{i} ∈ [0, 1], and points x_{μ_{i}}, each of them an ε-solution to the problem

minimize ⟨ μ_{i}c_{1} + (1 − μ_{i})c_{2}, x ⟩
subject to x ∈ G,

is

O(1) (\frac{\sqrt{ϑ}}{δ}(ln ϑ + ln(1/ε) + ln R + ln δ) + \sqrt{ϑ}(ln K − ln R)) ,

where the O-constant depends on λ_{1}, λ′_{1}, λ_{2}, λ_{3}, λ′_{3}, λ̂_{3}.

In contrast to this, the standard approach consists of solving each problem by starting from the same starting point (the approximation to the analytic center) without reusing any information. This strategy needs

O(1)\frac{1}{δ}\sqrt{ϑ}(ln ϑ + ln(1/ϵ) + ln K)

Newton steps. Clearly, for δ small enough, the proposed method has a higher efficiency.

Note that in each step of the algorithm, each of the points x_{μ} is an ε̃-solution to the problem parameterized by μ, where ε̃ > 0 is an accuracy measure implicitly defined by the number of Newton steps already made. More precisely, if i is the number of Newton steps made, then ε̃ = O(ϑ/ exp(i/\sqrt{ϑ})) is the accuracy achieved.

We have chosen here a branching into two different branches. Due to this, the δcovering of the parameter set [0, 1] is refined in each step to a δ/2-covering, multiplying the number of parameters to be considered by a factor of 2. Clearly, other branching schemes can also be used, in which a δ-covering is refined to a δ/ℓ-covering, ℓ > 2.

3.6 More than Two Criteria

Consider the case in which n > 2 criteria are given, i. e. a cone K ⊂ IR^{n} is specified together with a base B ⊂ K^{+} of K^{+}. It is rather easy to adapt the algorithm described in Section 3.5.2 to this general case. However, we need to discretize the n − 1-dimensional set B ⊂ IR^{n}, which is computationally a nontrivial task. Moreover, it has already been mentioned in the introduction that we will not get rid of the factor O(1/δ^{n−1}) in the computational complexity of any algorithm considered. Therefore, instead of pursuing further theoretical run-time estimates, some strategies about proper implementation techniques for the scheme presented are in order. The main question that we have to answer is: how do we discretize the base B?

3.6.1 Voronoi-Tesselations

Let μ denote the coordinates of a parameterization of B. Usually, μ ∈ IR^{n−1}. For example, if K = IR_{+}^{n}, the set cl(B) can be parameterized in barycentric coordinates by the convex hull of the origin and the n − 1 Euclidean unit vectors in IR^{n−1}. (And then, B is parameterized by the relative interior of this hull.) Let ω(μ) ∈ B be the element of B parameterized by μ and define

f_{μ}(x) := ⟨ ω(μ), x ⟩,
F_{t}^{μ}(x) := tf_{μ}(x) + F (x).

We can start with one parameter μ and define δ(0) := sup_{ω∈B} ∥ω−ω(μ)∥_{2}. Again, the set S = S(j) denotes the set of parameters μ which are considered as problem parameters in the path-following scheme. This set induces a Voronoi-tesselation (V_{μ})_{μ∈S} of B. Since the variable x_{μ} ∈ G (μ ∈ S) is supposed to represent the solution to all problems in the Voronoi cell V_{μ}, the stopping criterion of the WHILE-loop in Step 6 (a) ii has to be adapted to

6 (a) ii. WHILE λ(F_{t_{μ}}^{ν}, x_{μ}) < λ̂_{3} for all ν ∈ V_{μ} AND NOT stopping criterion for problem with parameter μ fulfilled DO

Moreover, the set S has to be refined in Step 6 (d) according to some strategy, i. e. points have to be added or at least one point has to be replaced by two or more points. If a standard Cartesian grid is used, we arrive at similar complexity estimates as their counterparts which have been derived in Section 3.5.3.

However, this algorithmic concept is rather problematic, since it is in no way clear how the abovementioned stopping criterion can be implemented in an efficient way, given that we have to compute (or rather update) Voronoi-tesselations in IR^{n−1}. Therefore, different strategies, not necessarily based on Voronoi diagrams might be in order.

3.6.2 Decomposition-based Strategies

A more feasible strategy is to decouple the point ω(μ) from the subset of the parameter set which it represents. To this end, let the set S = S(j) consist of pairs (ω(μ), S_{μ}) with ω(μ) ∈ S_{μ} ⊆ B such that the finite family (S_{μ})_{μ} defines a decomposition of B. Clearly, Voronoi tesselations are a special case, but the additional degrees of freedom imposed by the possibility of choosing ω(μ) ∈ S_{μ} arbitrarily as well as choosing the geometry of the decomposition has some distinct advantages. For example, we might choose an arbitrary refinement of the sets S_{μ} when replacing S(j) by S(j + 1). The most simple example in the standard case K = IR_{+}^{n} would be to start with S(0) = {((1/n)e, B)} and to use simplicial refinements only. Techniques for refining simplicial decompositions are well developed, see, e. g. [52, 51, 94]. Whenever a set S_{μ} which is part of a decomposition of B is polyhedral, Step 6 (a) ii can be replaced by

6 (a) ii. WHILE λ(F_{t_{μ}}^{ν}, x_{μ}) < λ̂_{3} for all vertices ν of S_{μ} AND NOT stopping criterion for problem with parameter μ fulfilled DO ,

a formulation which has obvious advantages for simplices S_{μ}.

3.6.3 Adaptive Grids

The strategies discussed in Section 3.6.1 and 3.6.2 both have their merits when we try to discretize the base B up to a certain discretization size. (Which can be measured, e. g. by max_{S_{μ}} diam(S_{μ}), the maximal diameter of the sets considered.) However, we have to keep in mind that the ultimate goal is not a discretization of B, but of E(M, K). The only interconnection between B and E(M, K) is the function φ : B ⟶ E_{p}(M, K), for which the results from Theorem 2.2.6 and 2.2.8 hold. Therefore, an estimate of the quality of a discrete approximation of E(M, K) should not be based on the quality of a discretization of B (which is a trivial problem). Instead, a bad discretization of E_{p}(M, K) should be refined locally wherever the quality is worst, thereby inducing a further refinement of the given discretization of B.

One possibility to achieve this is to note that usually E(M, K) ⊆ bd(M) holds [45, Satz 2.9], i. e. the set of efficient points has at most the same dimension as B, and that therefore a simplicial decomposition of E(M, K) is in order. For example, let (Ŝ_{i})_{i} be a triangulation of a surface in IR^{n} such that all vertices φ(ω_{j}) of this triangulation lie in E(M, K). Moreover, let there be given a topologically equivalent triangulation of the preimages ω_{j} in B. This triangulation in B induces a one-to-one relation between the set of n − 1-dimensional simplices (Ŝ_{i})_{i} and a corresponding set of n − 1-dimensional simplices (S_{i})_{i} in B. With this, one example of the refinement step is

6. (d) Let Ŝ_{i} be the n − 1-dimensional simplex with largest volume and S_{i} the corresponding simplex in B. Divide S_{i} into subsimplices S_{i}^{(j)} (j = 1, …, k). Denote by ω^{(ℓ)} (ℓ = 1, …, L) the new vertices. Replace S_{i} by the S_{i}^{(j)}, add the ω^{(ℓ)} to the set of ω_{j} and keep track of the refined triangulation.

Denoting by vol_{n−1} the n − 1-dimensional Lebesgue measure of a set, we see that max_{i}vol_{n−1}(Ŝ_{i}) is a global accuracy measure for the approximation to E(M, K) just constructed, and this measure is monotonically decreasing as the algorithm proceeds.

A simplicial approximation of E(M, K) has advantages with respect to a visualization of the approximation. Two possibilities are immediately at hand. First, each discrete point computed lying in E(M, K) can be represented by a point in the plane. These planar points can then be connected according to the edge structure of the decomposition. The second method computes the dual representation of the graph computed by the first method. There, every simplex of the approximation of E(M, K) is represented by a point in the plane, and these points are then connected according to the neighborhood structure of the simplices. Both graph representations are minimal in the sense that each vertex of the computed graph is connected to exactly n other vertices, and both representations map a neighborhood structure in the set E(M, K) onto edges of a graph.

56

Generalized Convex Goal Programming

4.1 Introduction

Recall from the discussion in Section 2.2.3.3 that typical scalarization techniques for multicriteria optimization problems lead to standard single-criterion optimization problems in which a norm has to be minimized over some set. But these problems occur also quite naturally in a large number of different applications (let us mention here just least-squares problems [8] and facility location [20] for the finite-dimensional case and almost arbitrary problems from approximation theory as well as control theory for the infinite-dimensional case). It is therefore important to gain as much insight into their structure as possible, while keeping the theoretical framework as general as possible. This is the purpose of this chapter.

4.1.1 Goal Programming

Goal Programming is a special scalarization technique in which certain norms or gauges γ : IR^{n} ⟶ IR are employed as functions monotone on K. Here, monotonicity is understood in the sense of Definition 2.2.1 with S := K. Consequently, Theorem 2.2.3 can be used to deduce that a minimum of γ over a set M ⊆ K is efficient. Moreover, scalarization techniques leading to a special class of gauges are able to characterize all efficient points as solutions to corresponding goal programs, even if the feasible set is not convex [56]. It is due to this that the study of the corresponding mathematical programs is important for multicriteria optimization. Note, however, that in applications the notion of Goal Programming is often reduced to the computation of just one efficient point, in contrast to the strategy outlined in the last chapter. Goal programs occur also in linearizations of nonconvex multicriteria problems, see Chapter 5 and especially Section 5.6.

The origins of Goal Programming date back to the work of Charnes, Cooper and Ferguson [15], where an ℓ_{1}-estimation regression model was proposed to estimate executive compensation. Since then, and thanks to its versatility and ease of use, it has become the by far most popular technique for tackling (linear) multiple-objective problems, as evidenced by the bulk of literature on theory and applications of the field. See, e. g., [79, 84, 90, 89] and the categorized bibliography of applications therein.

By a Non-Preemptive Goal Programming problem one usually means some particular instance of the following model: a polyhedron P ⊆ IR^{m} is given as the set of decisions; there exist N criteria matrices, C_{1}, …, C_{N}, with C_{j} in IR^{m×n_{j}}; each decision x ∈ P is valued according to criterion C_{j} by the vector C_{j}^{⊤}x, to be compared with a given target set T_{j} ⊆ IR^{n_{j}}. With this, the deviation d_{j}(x) of decision x with respect to the target set T_{j} is defined as

d_{j}(x) = inf_{z_{j}∈T_{j}} γ_{j}(C_{j}^{⊤}x − z_{j})

for some given norm γ_{j}, while the overall deviation at x is measured by

γ(d_{1}(x), …, d_{N}(x)),

where γ is a norm in IR^{N} assumed to be (IR_{+}^{N}, IR_{+})-monotone on IR_{+}^{N}, see [6, 57]. Then, the solution(s) minimizing the overall deviation are sought. In other words, one solves the convex program

inf_{x∈P} γ(d_{1}(x), …, d_{N}(x)).  (4.1)

As pointed out e. g. in [19, 78, 79], Non-Preemptive Goal Programming and related models can be rephrased as minimum-distance problems. This follows from the previous formulation, since (4.1) is equivalent to

minimize γ(γ_{1}(C_{1}^{⊤}x − z_{1}), …, γ_{N}(C_{N}^{⊤}x − z_{N}))
subject to x ∈ P,
z_{j} ∈ T_{j} (j = 1, …, N).

(4.2)

Denoting by γ̃ the norm in IR^{n_{1}} × … × IR^{n_{N}} defined as

γ̃(u_{1}, …, u_{N}) = γ(γ_{1}(u_{1}), …, γ_{N}(u_{N})),

problem (4.2) can be written as the minimum γ̃ -norm problem

minimize γ̃(u_{1}, …, u_{N})
subject to u_{j} = C_{j}^{⊤}x − z_{j}  (j = 1, …, N)
(x, z) ∈ P × ∏_{j=1}^{N} T_{j}

(4.3)

In many applications, each criterion C_{j} is assumed to be a vector c_{j} ∈ IR^{m}, so it values x through the scalar ⟨ c_{j}, x ⟩; each target set T_{j} is then a subset of IR of one the forms

T_{j} = [t_{j}, +∞[,  (4.4)
T_{j} = ] − ∞, t_{j}],  (4.5)
T_{j} = {t_{j}},  (4.6)

or, in Goal Range Programming [47], of the form

T_{j} = [\underline{t}_{j}, \overline{t}_{j}].  (4.7)

This corresponds to a goal constraint of type ⟨ c_{j}, x ⟩ ≥ t_{j}, ⟨ c_{j}, x ⟩ ≤ t_{j}, ⟨ c_{j}, x ⟩ = t_{j}, or \underline{t}_{j} ≤ ⟨ c_{j}, x ⟩ ≤ \overline{t}_{j}, respectively. In other words, one desires to have ⟨ c_{j}, x ⟩ above t_{j}, below t_{j}, exactly at t_{j}, or between \underline{t}_{j} and \overline{t}_{j}, respectively.

Whereas the choice of the aggregating norm γ is crucial, (although, in applications, mostly reduced to the cases ∥ · ∥_{1} or ∥ · ∥_{∞}) the choice of the norms γ_{j} to measure deviations in the case n_{j} = 1 ∀ j is irrelevant, and we can consider each γ_{j} to be equal to the absolute value function. Then, the deviations take on the more familiar form

d_{j}(x) = \{ \begin{matrix}max \{t_{j} − ⟨ c_{j}, x ⟩, 0\} & if T_{j} = [t_{j}, +∞[, \\ max \{⟨ c_{j}, x ⟩ − t_{j}, 0\} & if T_{j} =] − ∞, t_{j}], \\ |⟨ c_{j}, x ⟩ − t_{j}| & if T_{j} = \{t_{j}\}, \\ max \{\underline{t}_{j} − ⟨ c_{j}, x ⟩, 0\} + max \{⟨ c_{j}, x ⟩ − \overline{t}_{j}, 0\} & if T_{j} = [\underline{t}_{j}, \overline{t}_{j}]. \end{matrix}

From these expressions, it should become clear that target sets of type (4.7), (thus also of type (4.6)) are used only for modeling convenience, since they can be derived from sets of types (4.4) and (4.5): splitting criterion j into criteria j_{1}, j_{2}, and defining T_{j}^{1} = [\underline{t}_{j}, +∞[ and T_{j}^{2} = ] − ∞, \overline{t}_{j}], the deviation d_{j}(x) is simply the sum of the deviations with respect to T_{j}^{1} and T_{j}^{2}.

4.1.2 Examples

Applications of Goal Programming abound in the literature; see e. g. the list of 351 applications papers cited in [79]. However, the range of applicability of (4.1) is by no means reduced to what is usually classified as Goal Programming: a vast series of important models in different fields of Optimization can also be seen as particular instances of (4.1), mainly from the perspective of minimum-distance problems. Some of them are briefly discussed below.

4.1.2.1 Overdetermined systems of (in)equalities

If a system of linear equalities and inequalities

\begin{aligned}⟨ a_{1}, x ⟩ & ≥ & b_{1} \\ ⟨ a_{2}, x ⟩ & ≥ & b_{2} \\ ⋮ & ≥ & ⋮ \\ ⟨ a_{p}, x ⟩ & ≥ & b_{p} \\ ⟨ a_{p+1}, x ⟩ & = & b_{p+1} \\ ⋮ & = & ⋮ \\ ⟨ a_{p+q}, x ⟩ & = & b_{p+q} \end{aligned}  (4.8)

is infeasible, one can look for a so-called least infeasible solution, i. e. a point x^{∗} solving

min_{x} γ(max(0, b_{1} − ⟨ a_{1}, x ⟩), …, max(0, b_{p} − ⟨ a_{p}, x ⟩),

|b_{p+1} − ⟨ a_{p+1}, x ⟩|, …, |b_{p+q} − ⟨ a_{p+q}, x ⟩|)

for some norm γ assumed to be (IR_{+}^{p+q}, IR_{+})-monotone on IR_{+}^{p+q}. This is simply a Goal Programming problem in which the vectors a_{i} (i = 1, …, p + q) play the role of the criteria and the components b_{i} (i = 1, …, p + q) of the right hand side vector represent the targets, see Example 4.3.4 in Section 4.3.

When only equalities appear in (4.8), one obtains the problem of solving an overdetermined system of linear equations, classical in Approximation Theory [69, 86], or, equivalently, the Linear Regression problem [85]. Usually, γ is assumed to be an ∥ · ∥_{p}-norm, mainly p = 2 (yielding the well-known Least Squares problem [8]), p = 1, or p = ∞, see [2].

4.1.2.2 Multifacility location

In Continuous Location [64, 72], distances are usually measured by gauges. For simplicity, we will consider throughout this chapter only gauges γ of the form

γ(x) = inf{t ≥ 0 : x ∈ tB}  (4.9)

for some nonempty convex compact B ⊂ IR^{N} (its unit ball) containing the origin in its interior. In applications, this additional assumption is usually fulfilled, see, e. g. [23, 64]. Observe that norms correspond to symmetric gauges. Moreover, since the origin is assumed to be an interior point, the gauge takes always finite values. See e. g. [40] for the case of gauges with values on IR_{+} ∪ {+∞}.

Let F be a nonempty finite set and let ∅ ≠ E ⊆ F × F. Then (F, E) is a directed graph. Following e. g. [31, 61], F represents the set of facilities (some of which may have fixed locations in IR^{N}), whereas E represents the interactions between these facilities.

For each edge e := (f, g) ∈ E, let γ_{e} be a given gauge in IR^{N}, which measures the cost of the interaction between facility f and facility g. Let γ be a gauge in IR^{E} which is (IR_{+}^{E}, IR_{+})-monotone in IR_{+}^{E}.

For a nonempty closed convex set G ⊆ (IR^{N})^{F}, consider the optimization problem

inf_{(x_{f})_{f∈F} ∈G} γ((γ_{(f,g)}(x_{f} − x_{g}))_{(f,g)∈E}).  (4.10)

The most popular instance of (4.10) is the continuous minisum multifacility location problem, see [72, 95, 96] and the references therein. There, the node set F is partitioned into two sets A and V , representing respectively the fixed and the free locations, and a family (a_{f})_{f∈A} ∈ (IR^{N})^{A} of fixed locations is given. The feasible region G is then defined by

G = {x = (x_{f})_{f∈F} ∈ (IR^{N})^{F} | x_{f} = a_{f} for all f ∈ A} ,  (4.11)

while the gauge γ is taken as the ∥ · ∥_{1}-norm, so that one minimizes the sum of all interactions between the facilities,

inf_{{x_{f} =a_{f} \atop ∀f∈A}} ∑_{(f,g)∈E} γ_{(f,g)}(x_{f} − x_{g}).  (4.12)

Let J_{(F,E)} be the incidence matrix of the graph (F, E), i. e. J_{(F,E)} ∈ IR^{E×F} is the matrix in which the row e := (f, g) ∈ E has zeroes in all its positions except in the position f, where the entry is 1, and in position g, where the entry is −1. Moreover, define the matrix C by C := J_{(F,E)} ⊗ I_{N}, the Kronecker product of J_{(F,E)} with the unit matrix I_{N} ∈ IR^{N×N}. Let γ be the gauge in (IR^{N})^{E} defined by

γ : u := (u_{e})_{e∈E} ⟼ γ(u) := ∥(γ_{e}(u_{e}))_{e∈E}∥_{1}
= ∑_{e∈E} γ_{e}(u_{e}).

Then, problem (4.12) can be written as

minimize γ(Cx) = γ(((Cx)_{e})_{e∈E})
subject to x ∈ G ⊆ (IR^{N})^{F} ,

(4.13)

which is a particular instance of (4.1).

A similar representation can be obtained for the continuous minimax multifacility location problem [54], in which expression (4.13) holds for γ defined by

γ : u := (u_{e})_{e∈E} ⟼ γ(u) := ∥(γ_{e}(u_{e})_{e∈E})∥_{∞}
= max_{e∈E} γ_{e}(u_{e})

General monotone gauges γ have been suggested by Durier [22, 23]. In the latter paper, he introduced problems with fixed costs, which can also be accommodated within this framework. Indeed, for

inf_{(x_{f})_{f∈F} ∈G} γ ((ω_{(f,g)} + γ_{(f,g)}(x_{f} − x_{g}))_{(f,g)∈E})

with a given vector (ω_{e})_{e∈E} ∈ IR^{E} with non-negative components, one may write

inf_{(x_{f})_{f∈F} ∈G} γ ((γ̂_{(f,g)}(ω_{(f,g)}, x_{f} − x_{g}))_{(f,g)∈E}) ,

where each γ̂_{e} is a gauge in IR × IR^{N} defined by

γ̂_{e}(ω, z) = |ω| + γ_{e}(z).

With this, again an expression of type (4.13) is obtained.

Our aim is to study a generalized version of Problem (4.1) under some mild assumptions on the feasible set G, namely, G will be assumed to be an asymptotically conical set. To do this, we have structured the rest of this chapter as follows: In Section 4.2 the concept of asymptotically conical sets is introduced, and some elementary properties are discussed. Then, in Section 4.3, the problem under study is formally defined and its dual is derived. In Section 4.4, the existence of primal and dual optimal solutions is studied in detail, giving, in particular, sufficient conditions for the attainment of the optimal value. Then, an Interior-Point method is described in Section 4.5, yielding a unified methodology for solving problems which, until now, were solved by different (some not polynomial) techniques.

4.2 Asymptotically Conical Sets and their Properties

Definition 4.2.1 (Asymptotically Conical Sets) A nonempty set S ⊆ IR^{m} is said to be asymptotically conical if it admits a representation of the form

S = M + E,  (4.14)

for some compact convex set M and some closed convex cone E. In such a case, the pair (M, E) is an asymptotically conical representation (a. c. r.) of S.

The optimization problem addressed in this chapter will have an asymptotically conical set as its set of feasible points, see Section 4.3. Here, we take a quick look at the basic properties of such class of sets.

Denote by G_{∞} the recession cone of G,

G_{∞} = {y ∈ IR^{m} : G + y ⊆ G} ,

see Theorem 8.1 in [77]. We now have the following properties.

Lemma 4.2.1 If (M, E) is an a. c. r. of G then

G_{∞} = E.  (4.15)

Lemma 4.2.2 Let the sets G_{1}, G_{2} ⊆ IR^{m} be asymptotically conical and let (M_{1}, E_{1}), (M_{2}, E_{2}) be the corresponding asymptotical conical representations. Then,

1. (M_{1} × M_{2}, E_{1} × E_{2}) is an a. c. r. of the asymptotically conical set G_{1} × G_{2}.

2. (M_{1} + M_{2}, E_{1} + E_{2}) is an a. c. r. of the asymptotically conical set G_{1} + G_{2}.

Lemma 4.2.3 Let A be an affine transformation of the form A(x) = Ax + b with a matrix A and a vector b. Then, (AM + b, AE) is an a. c. r. of A(G).

Remark 4.2.4 Compact sets, polyhedra, affine spaces, and cones are asymptotically conical. Although each of these classes is closed under intersections, this is not the case of the whole class of asymptotically conical sets. Indeed, take e. g. the following two asymptotically conical sets in IR^{3}:

S_{1} = {(x_{1}, x_{2}, x_{3})^{⊤} | x_{3}^{2} ≥ x_{1}^{2} + x_{2}^{2}, x_{3} ≥ 0}

and S_{2} = {(x_{1}, x_{2}, x_{3})^{⊤} | x_{1} = 1}, whereas, by Property 4.2.1, no a. c. r. for the hyperbola S_{1} ∩ S_{2} exists. Moreover, this example shows that the inverse image of an asymptotically conical set under an affine mapping is not necessarily asymptotically conical.

□

Finally, we have the following simple observation.

Lemma 4.2.5 Let (M, E) be an a. c. r. of G. Then, for any u ∈ IR^{m},

inf_{x∈G} ⟨ u, x ⟩ = \{ \begin{matrix}min_{x∈M}⟨ u, x ⟩ & ifu ∈ E^{∗}, \\ −∞ & else. \end{matrix}  (4.16)

4.3 The Problem Addressed and its Dual

The problem addressed in this chapter has the form

inf g(x) := γ(Cx + c) + ⟨ d, x ⟩,
subject to x ∈ G,

(4.17)

where γ is a gauge in the IR^{n}, C is a matrix in IR^{n×m}, c ∈ IR^{n} and d ∈ IR^{m} are vectors, and G = M + E ⊆ IR^{m} is a nonempty asymptotically conical set with asymptotical conical representation (M, E).

Observe that, in particular, problem (4.17) contains as instances all the examples discussed in Section 4.1.2 and most of the scalarizations discussed in Chapter 2 (see especially Section 2.2.3.3). Moreover, the case G = IR^{m} has been addressed in [59], whereas the case d = 0 leads to the so-called gauge- or homogeneous program, addressed, among others in [40, 25, 44, 48].

In these references, duals are derived and Slater-type assumptions are made to link primal and dual optimality. We show below and illustrate by examples how the knowledge of an a. c. r. can be successfully used to address duality questions and to design efficient algorithms as well.

A dual for (4.17) can easily be derived using minimax theorems as basic tool. Indeed, one has

inf_{x∈G} γ(Cx + c) + ⟨ d, x ⟩ = inf_{x∈G} sup_{γ^{∘}(u)≤1} (⟨ u, Cx + c ⟩ + ⟨ d, x ⟩)  (4.18)
= sup_{γ^{∘}(u)≤1} (⟨ u, c ⟩ + inf_{x∈G}⟨ C^{⊤}u + d, x ⟩)  (4.19)
= sup_{{γ^{∘}(u)≤1, \atop C^{⊤}u+d∈E^{∗}}} (⟨ u, c ⟩ + min_{x∈M}⟨ C^{⊤}u + d, x ⟩) , (4.20)

where (4.18) follows from the representation of a gauge as the support of its polar unit ball, see Theorem 14.5 of [77], and (4.19) follows from the Minimax Theorem stated as Corollary 37.3.2 in [77] and the fact that γ^{∘}, the dual gauge of γ, has compact level sets (recall that we are assuming that gauges γ has the origin in its interior, which guarantees the compactness of its dual ball). Finally, (4.20) follows from Property 4.2.5.

Denoting by δ_{S}^{∗} the support of a set S ⊆ IR^{m},

δ_{S}^{∗}(x) = sup {⟨ x, y ⟩ | y ∈ S} ,

the chain of equalities above yields

inf  γ(Cx + c) + ⟨ d, x ⟩ = max  ⟨ u, c ⟩ − δ_{M}^{∗}(−C^{⊤}u − d)
subject to x ∈ M + E  subject to C^{⊤}u + d ∈ E^{∗}
γ^{∘}(u) ≤ 1

(4.21)

From this equivalence, we will call the optimization problem

maximize ⟨ u, c ⟩ − δ_{M}^{∗}(−C^{⊤}u − d)
subject to C^{⊤}u + d ∈ E^{∗},
γ^{∘}(u) ≤ 1

(4.22)

from the right-hand side of (4.21) the dual of problem (4.17) and we have already shown that (4.17) and (4.22) have identical optimal value.

Before exploring further the relations between (4.17) and (4.22) we now present some particular instances of (4.17) whose corresponding dual (4.22) has a simple (explicit) form.

Example 4.3.1 (Compact Feasible Sets) Let γ̃ be a gauge in IR^{m}, let x_{0} ∈ IR^{m} and let M = {x ∈ IR^{m} | γ̃(x − x_{0}) ≤ r} for some constant r ≥ 0. Since, by definition of dual gauges,

min {⟨ C^{⊤}u + d, x ⟩ | γ̃(x − x_{0}) ≤ r} = ⟨ C^{⊤}u + d, x_{0} ⟩ − rγ̃^{∘}(−C^{⊤}u − d),

we get the dual

maximize ⟨ u, c ⟩ + ⟨ C^{⊤}u + d, x_{0} ⟩ − rγ̃^{∘}(−C^{⊤}u − d)
subject to C^{⊤}u + d ∈ E^{∗},
γ^{∘}(u) ≤ 1.

(4.23)

Example 4.3.2 (Unconstrained Problems) Setting G = IR^{m} and M = {0}, we have E = IR^{m} and E^{∗} = {0}. Hence, the dual (4.22) takes the form

maximize ⟨ u, c ⟩
subject to C^{⊤}u + d = 0,
γ^{∘}(u)  ≤ 1.

This dual has been derived in [59] using the same idea but lengthier arguments, see their Theorem 1, Remark 4 and Remark 5. □

Example 4.3.3 (Distance between two Sets) Given asymptotically conical sets G_{1} ⊆ IR^{m} and G_{2} ⊆ IR^{m}, with respective a. c. r. G_{1} = M_{1} + E_{1} and G_{2} = M_{2} + E_{2} and a gauge γ in IR^{m}, formula (4.21) provides an alternative expression for the distance δ_{γ}(G_{1}, G_{2}) between G_{1} and G_{2}. Indeed, since

δ_{γ}(G_{1}, G_{2}) = inf {γ(x_{1} − x_{2}) | x_{1} ∈ G_{1}, x_{2} ∈ G_{2}}
= inf {γ(Cx) | x ∈ G} ,

for G := G_{1} × G_{2} and C = (I_{m}, −I_{m}) and I_{m} the m × m identity matrix, one gets

δ_{γ}(G_{1}, G_{2}) = max_{u∈E_{1}^{∗}∩(−E_{2}^{∗})∩B^{∘}} ( min_{x_{1}∈M_{1},x_{2}∈M_{2}}⟨ u, x_{1} − x_{2} ⟩)
= max_{u∈(E_{1}^{∗}∩(−E_{2}^{∗})∩B^{∘})} min_{x∈M_{1}−M_{2}}⟨ u, x ⟩

where B^{∘} is the unit ball of the gauge γ^{∘}. When G_{1} is an affine manifold, G_{1} = {p_{0}}+E_{1} for some vector space E_{1} ⊆ IR^{m}, and G_{2} is a cone, thus having ({0}, G_{2}) as a. c. r. decomposition, we get

δ_{γ}(G_{1}, G_{2}) = max {⟨ u, p_{0} ⟩ | u ∈ E_{1}^{∗} ∩ (−E_{2}^{∗}) ∩ B^{∘}}

This expression yields a simple characterization for G_{1} ∩ G_{2} ≠ ∅. Indeed, for any γ one has that G_{1} ∩ G_{2} ≠ ∅ if and only if δ_{γ}(G_{1}, G_{2}) ≤ 0 (see Corollary 4.4.10), thus (by a slight abuse of notation)

G_{1} ∩ G_{2} ≠ ∅ if and only if ⟨ u, p_{0} ⟩ ≤ 0  ∀u ∈ E_{1}^{∗} ∩ (−E_{2}^{∗}) ∩ B^{∘}
if and only if ⟨ u, p_{0} ⟩ ≤ 0  ∀u ∈ E_{1}^{∗} ∩ (−E_{2}^{∗})
if and only if p_{0} ∈ − (E_{1}^{∗} ∩ (−E_{2}^{∗}))^{∗}
□

The following result will be useful to rephrase the dual (4.22) if the gauge γ in use is a composite gauge:

Lemma 4.3.1 (Composite Gauges) Let γ_{1}, …, γ_{k} be gauges in IR^{n_{1}}, …, IR^{n_{k}}, and let γ̃ be a gauge in IR^{k}, monotone in IR_{+}^{k}. The gauge γ in IR^{n_{1}} × ⋯ × IR^{n_{k}} defined by

γ(u_{1}, …, u_{k}) = γ̃(γ_{1}(u_{1}), …, γ_{k}(u_{k}))  (4.24)

has as dual the gauge γ^{∘} with

γ^{∘}(x_{1}, …, x_{k}) = γ̃^{∘}(γ_{1}^{∘}(x_{1}), …, γ_{k}^{∘}(x_{k})).

Proof: Let x := (x_{1}, …, x_{k}) ∈ IR^{n_{1}} × … × IR^{n_{k}}. Since γ̃ is monotone in IR_{+}^{k}, by definition of a dual gauge one has

γ̃^{∘}(γ_{1}^{∘}(x_{1}), …, γ_{k}^{∘}(x_{k}))
= max { ∑_{j=1}^{k} α_{j}γ_{j}^{∘}(x_{j}) | α ∈ IR_{+}^{k}, γ(α) = 1}
= max { ∑_{j=1}^{k} α_{j}⟨ u_{j}, x_{j} ⟩ | α ∈ IR_{+}^{k}, γ(α) = 1, γ_{j}(u_{j}) = 1 ∀j}
= max { ∑_{j=1}^{k}⟨ α_{j}u_{j}, x_{j} ⟩ | α ∈ IR_{+}^{k}, γ(α_{1}γ_{1}(u_{1}), …, α_{k}γ_{k}(u_{k})) = 1}
= max { ∑_{j=1}^{k}⟨ α_{j}u_{j}, x_{j} ⟩ | α ∈ IR_{+}^{k}, γ(γ_{1}(α_{1}u_{1}), …, γ_{k}(α_{k}u_{k})) = 1} .

With the change of variables α_{j}u_{j} = ω_{j} (j = 1, …, k), we get (abusing the notation again)

γ̃^{∘}(γ_{1}^{∘}(x_{1}), …, γ_{k}^{∘}(x_{k})) = max { ∑_{j=1}^{k}⟨ ω_{j}, x_{j} ⟩ | γ(γ_{1}(ω_{1}), …, γ_{k}(ω_{k})) = 1}
= (γ̃(γ_{1}(x_{1}), …, γ_{k}(x_{k})))^{∘} .
□

This lemma yields a very simple dual for gauges of the form (4.24):

Corollary 4.3.2 Let C_{i} ∈ IR^{m_{i}×n} (i = 1, …, k) be matrices and define the matrix C by C^{⊤} := (C_{1}^{⊤}, C_{2}^{⊤}, …, C_{k}^{⊤}). Let c_{i} ∈ IR^{m_{i}} (i = 1, …, k) be vectors and set c^{⊤} := (c_{1}^{⊤}, …, c_{k}^{⊤}). Moreover, let γ be defined as in (4.24) and let G be asymptotically conical with a. c. r. (M, E). Then (4.22) admits the form

maximize ∑_{j=1}^{k}⟨ u_{j}, c_{j} ⟩ + min_{x∈M} ⟨x, ∑_{j=1}^{k} C_{j}^{⊤}u_{j} + d⟩
subject to ∑_{j=1}^{k} C_{j}^{⊤}u_{j} + d ∈ E^{∗},
γ̃^{∘}(γ_{1}^{∘}(u_{1}), …, γ_{k}^{∘}(u_{k})) ≤ 1

We illustrate the power of our strategy for deriving the dual by applying it to two problems previously addressed in the literature, as discussed in the following examples.

Example 4.3.4 (Infeasible Flow Problems) We consider the flow problem of [63]. Let (F, E) be a directed graph. (Here, the edge set E is not to be confused with the conic part of an a. c. r. for a given convex subset of the IR^{m}.) Associate with each arc e ∈ E a lower bound l_{e} and an upper bound u_{e} on its capacity, l_{e} ∈ [−∞, +∞[, u_{e} ∈] − ∞, +∞]. Associate with each node f ∈ F its demand d_{f} ∈ IR. Flows on (F, E) are vectors x in IR^{E}; a flow x is called feasible if it satisfies both flow conservation,

∑_{f: (f,g)∈E} x_{(f,g)} − ∑_{f: (g,f)∈E} x_{(g,f)} = d_{g}  ∀ g ∈ F,  (4.25)

and boundedness,

\begin{aligned}x_{e} & ≥ & l_{e} \\ x_{e} & ≤ & u_{e} \end{aligned}  ∀ e ∈ E.  (4.26)

When no feasible flow exists, McCormick proposed in [63] to consider (4.25)-(4.26) as goal constraints and to solve the corresponding problem (4.1) for γ equal to the (weighted) ∥ · ∥_{1}-, ∥ · ∥_{∞}-, and ∥ · ∥_{2}-norm.

We first reformulate (4.26) as distance constraints following (4.2):

\begin{aligned}x_{e} & ∈ & [l_{e}, +∞[ \\ x_{e} & ∈ & ] − ∞, u_{e}] \end{aligned}  ∀e ∈ E.  (4.27)

Then, the problem can be written as

minimize γ (( ∑_{f: (f,g)∈E} x_{(f,g)} − ∑_{f: (g,f)∈E} x_{(g,f)} − d_{g})_{g∈F}, (x_{e} − \underline{u}_{e})_{e∈E} , (x_{e} − \overline{l}_{e})_{e∈E})
subject to \underline{u}_{e} ∈] − ∞, u_{e}]  ∀ e ∈ E,
\overline{l}_{e} ∈ [l_{e}, +∞[ ∀ e ∈ E,

or

minimize γ (( ∑_{f: (f,g)∈E} x_{(f,g)} − ∑_{f: (g,f)∈E} x_{(g,f)} − d_{g})_{g∈F}, (x_{e} − \underline{u}_{e})_{e∈E} , (x_{e} − \overline{l}_{e})_{e∈E})
subject to (x, \underline{u}, \overline{l}) ∈ \{((0)_{e∈E}, (u_{e})_{e∈E}, (l_{e})_{e∈E})\} + (IR^{E} × IR_{−}^{E} × IR_{+}^{E}) .  (4.28)

Associate dual variables π := (π_{g})_{g∈F}, π^{+} := (π_{e}^{+})_{e∈E}, and π^{−} := (π_{e}^{−})_{e∈E} with the three blocks of components in (4.28). Since, in the a. c. r. of (4.28), one has

(IR^{E} × IR_{−}^{E} × IR_{+}^{E})^{∗} = {0}^{E} × IR_{−}^{E} × IR_{+}^{E},

we obtain the dual

maximize ∑_{g∈F} π_{g}(−d_{g}) + min_{x=0,\underline{u}=u,\overline{l}=l} ∑_{e∈E}(−π_{e}^{+})\underline{u}_{e} + ∑_{e∈E}(−π_{e}^{−})\overline{l}_{e}
subject to π_{e}^{+} ≥ 0  ∀e ∈ E,
π_{e}^{−} ≤ 0 ∀e ∈ E,
π_{e}^{+} + π_{e}^{−} = π_{f} − π_{g} ∀ e = (f, g) ∈ E,
γ^{∘}(π, π^{+}, π^{−}) ≤ 1,

i. e.

maximize − ∑_{g∈F} π_{g}d_{g} − ∑_{e∈E} π_{e}^{+}u_{e} − ∑_{e∈E} π_{e}^{−}l_{e}
subject to π_{e}^{+} ≥ 0  ∀ e ∈ E,
π_{e}^{−} ≤ 0 ∀ e ∈ E,
π_{e}^{+} + π_{e}^{−} = π_{f} − π_{g}, ∀ e = (f, g) ∈ E,
γ^{∘}(π, π^{+}, π^{−}) ≤ 1,

an expression which includes the particular cases derived in [63]. □

As another application, we derive the dual of the quite general unconstrained multifacility location problem (4.10) introduced in Section 4.1.2.2.

Example 4.3.5 (Facility Location) With the notation as used in (4.10), for G defined in (4.11), one has that

G_{∞} = ∏_{f∈F} G^{f},

where G^{f} = {0} ⊂ IR^{N} if f ∈ A and G^{f} = IR^{N} if f ∉ A. Then, one obtains the dual

maximize ∑_{f∈A} ⟨a_{f}, ∑_{h: (f,h)∈E} u_{(f,h)} − ∑_{h: (h,f)∈E} u_{(h,f)}⟩
subject to ∑_{g: (f,g)∈E} u_{(f,g)} − ∑_{g: (g,f)∈E} u_{(g,f)} = 0  ∀f ∉ A,
γ^{∘}((γ_{e}^{∘}(u_{e}))_{e∈E}) ≤ 1,

which covers most of the instances previously addressed in the literature, e. g. [39, 44]. □

Remark 4.3.3 The minimax argument of (4.19) can also be used re-derive (4.21) from the basic theory of Fenchel conjugacy. Indeed, denoting by g the function

g(x) = γ(Cx + c) + ⟨ d, x ⟩,

and by δ_{G} the indicator function of G,

δ_{G}(x) = \{ \begin{matrix}0 & if x ∈ G, \\ +∞, & else, \end{matrix}

the conjugate of g + δ_{G} is given by

(g + δ_{G})^{∗} (p) = inf_{γ^{0}(u)≤1} {−⟨ u, c ⟩ + sup_{x∈G} ⟨x, −C^{⊤}u + p − d⟩} .  (4.29)

In particular, by inspecting −(g + δ_{G})^{∗} at 0, one obtains again the dual (4.22). □

Example 4.3.6 (Positive Polynomials) In algebra, questions of the following type frequently arise. Suppose that a subset S ⊂ P_{n} of the set P_{n} of univariate polynomials with maximal degree n is given. The question is if at least one of the polynomials in S is nonnegative on the whole IR. Now let E_{n} ⊂ P_{n} be the cone of nonnegative polynomials with maximal degree n. A natural reformulation of the question is then the optimization problem

minimize γ(p − q)
subject to p ∈ S,
q ∈ E_{n},

where γ is an arbitrary gauge on P_{n}. Let (M, R) be an a. c. r. of S and let p̂ ∈ int(M). The space P_{n} is (n + 1)-dimensional, and we choose as a basis the monomials 1, t, t^{2}, …, t^{n}. Each polynomial p(t) = ∑_{i=0}^{n} p_{i}t^{i} can then be identified with the vector (p_{0}, …, p_{n})^{⊤} ∈ IR^{n+1}. We interpret all sets and variables in the problem above as sets and vectors in the IR^{n+1} and use the standard scalar product on IR^{n+1}. The dual problem can then be written as

minimize γM̂^{∘}(−u) − ⟨ p̂, u ⟩
subject to u ∈ R^{∗} ∩ (−E_{n}^{∗}),

γ^{0}(u) ≤ 1, where M̂ = −p̂+ M and γM̂^{∘} is the gauge with unit ball M̂^{∘}. Note that we can describe the set E_{n} as the image of the cone of all positive semidefinite matrices in IR^{(n+1)×(n+1)} under a simple linear operator Λ : IR^{(n+1)×(n+1)} ⟶ IR^{n+1}, see Nesterov [66]. Similar results hold if the domain of the polynomials considered is restricted to a (possibly unbounded) interval [66]. While the primal problem can be used to construct a certificate that at least one of the polynomials in S has a strictly negative function value (by constructing such a polynomial and the corresponding function value, i. e. by solving the primal problem), the dual problem can be used to construct a certificate that no polynomial in S has nonpositive function values (by constructing a strictly positive lower bound for all function values of all polynomials in S, i. e. by solving the dual problem). Polynomial-time algorithms for the primal and the dual problem are discussed in Section 4.5. □

4.4 Existence of Primal and Dual Solutions

In this section we study the finiteness and the attainment of the optimal value of the problems (4.17) and (4.22).

Theorem 4.4.1 (Finite Optimal Values) Let B^{∘} be the polar of the unit ball B of γ. The following statements are equivalent.

1. (4.17) (and (4.22)) have finite optimal value.

2. For all y ∈ G_{∞} one has γ(Cy) + ⟨ d, y ⟩ ≥ 0.

3. d ∈ (G_{∞})^{∗} − C^{⊤}B^{∘}.

Proof: Denote by g_{∞} the recession function of g and let y ∈ IR^{m}. For arbitrary x ∈ IR^{m} we have that

g_{∞}(y) = sup_{λ>0} \frac{g(x + λy) − g(x)}{λ}
= lim_{λ→+∞} \frac{g(x + λy) − g(x)}{λ}
= lim_{λ→+∞} \frac{γ(C(x + λy) + c) + ⟨ d, x + λy ⟩ − γ(Cx + c) − ⟨ d, x ⟩}{λ}
= γ(Cy) + ⟨ d, y ⟩,  (4.30)

where the second equation is due to Theorem 8.5 of [77] and the last equation follows because of the homogeneity of γ. If (4.17) has a finite optimal value, then Part 2 is a consequence of (4.30) and Theorem 27.1 (parts (a) and (i)) of [77]. Conversely, if Condition 2 holds, we have for any a. c. r. (M, E) of G and for any x ∈ G, x = x_{M} +x_{E} with x_{M} ∈ M, x_{E} ∈ E = G_{∞}, that

γ(Cx + c) + ⟨ d, x ⟩ = γ(Cx_{M} + Cx_{E} + c) + ⟨ d, x_{M} ⟩ + ⟨ d, x_{E} ⟩
≥ γ(Cx_{E}) + ⟨ d, x_{E} ⟩ − γ(−Cx_{M} − c) + ⟨ d, x_{M} ⟩
≥ min_{x_{M} ∈M} (−γ(−Cx_{M} − c) + ⟨ d, x_{M} ⟩) > −∞,

showing that Condition 1 holds.

The equivalence between Condition 1 and Condition 3 follows from (4.21) and (4.15). □

Remark 4.4.2 Conditions 2 and 3 do not imply Condition 1 for sets G which are not asymptotically conical. As a simple counterexample, take G = {x = (x_{1}, x_{2})^{⊤} ∈ IR^{2} | x_{2} ≥ x_{1}^{2}}, let C be the 1 × 2 matrix C = (1, 0), d := (2, 0)^{⊤}, c := (0), and let γ(s) = |s| for all s ∈ IR. Then G_{∞} is the ray expanded by the vector (0, 1)^{⊤}, and thus γ(Cx) + ⟨ d, x ⟩ = 0 for all x ∈ K_{∞}. Hence, Condition 2 holds. Moreover,

(G_{∞})^{∗} = {(x_{1}, x_{2})^{⊤} | x_{2} ≥ 0},

thus, taking u = 0 ∈ B^{∘}, one obtains d ∈ (G_{∞})^{∗} − C^{⊤}B^{∘}, and Condition 3 holds also. However, g((−k, k^{2})^{⊤}) = γ(C(−k, k^{2})^{⊤} + c) + ⟨ d, (−k, k^{2})^{⊤} ⟩ = −k for every integer k, thus the optimal value of (4.17) is −∞. □

The duality scheme previously described enables us to easily characterize the (possibly empty) set of optimal solutions of (4.17) in terms of any optimal solution ū of (4.22). See also Theorem 1 of [40] for the case of a polyhedral feasible set G, or Theorem 1.1 of [44] for related constraint qualification assumptions.

Theorem 4.4.3 (Saddlepoints) Let one of the equivalent conditions of Theorem 4.4.1 hold. Then,

1. The set of optimal solutions of (4.22) is not empty.

2. Let x̄ be feasible for (4.17) and ū feasible for (4.22). Then x̄ is optimal for (4.17) and ū is optimal for (4.22) if and only if the pair (x̄, ū) is a saddle-point for the problem

inf_{x∈G} sup_{γ^{∘}(u)≤1} ⟨u, Cx + c⟩ + ⟨ d, x ⟩.  (4.31)

Proof: By (4.21), under the assumptions of Theorem 4.4.1, the dual (4.22) consists of the maximization of the continuous function

u ↦ ⟨ u, c ⟩ + min_{x∈M} ⟨x, C^{⊤}u + d⟩

over the nonempty compact set {u | C^{⊤}u + d ∈ E^{∗}, γ^{∘}(u) ≤ 1}. Thus an optimal solution ū for (4.22) always exists.

For Part 2, observe that, due to Theorem 4.4.1 and (4.21), the saddle value exists and is finite. Hence, saddle points exist, cmp. e. g. Theorem 4.2.5 in [50]. Moreover, the set of saddle points coincides with the Cartesian product of the set of optimal solutions for (4.17) and (4.22), as asserted. □

The characterization of optimal solutions of (4.17) as part of saddle points yields the following result.

Theorem 4.4.4 (Optimality Conditions) Let (M, E) be an a. c. r. of G. The feasible point x̄ := x̄_{M} + x̄_{E}, x̄_{M} ∈ M, x̄_{E} ∈ E, is optimal for (4.17) if and only if there exists a point ū ∈ IR^{n} satisfying

γ^{∘}(ū) ≤ 1,
C^{⊤}ū + d ∈ E^{∗} ∩ −N_{M}(x̄_{M}),
⟨ x̄_{E}, C^{⊤}ū + d ⟩ = 0,
Cx̄ + c ∈ N_{B^{∘}}(ū).

In that case, such a point ū is an optimal solution for (4.22).

Proof: By Theorem 4.4.3, x̄ is optimal for (4.17) if and only if there exists ū optimal for (4.22) such that the pair (x̄, ū) is a saddle point. In other words, x̄ is optimal for (4.17) if and only if there exists some ū ∈ IR^{n} satisfying

C^{⊤}ū + d ∈ E^{∗},  (4.32)
γ^{∘}(ū) ≤ 1,  (4.33)
γ(Cx̄ + c) + ⟨ d, x̄ ⟩ = ⟨ ū, Cx̄ + c ⟩ + ⟨ d, x̄ ⟩  (4.34)
= inf_{x_{M} ∈M,x_{E}∈E}⟨ ū, C(x_{M} + x_{E}) ⟩ + ⟨ d, x_{M} + x_{E} ⟩ (4.35)
+⟨ ū, c ⟩.

But (4.34) holds if and only if

Cx̄ + c ∈ N_{B^{∘}}(ū).

Moreover, for vectors u satisfying (4.32), it follows from Lemma 4.2.5 that

inf_{x∈G}⟨ u, Cx ⟩ + ⟨ d, x ⟩ = min_{x_{M} ∈M}⟨ u, Cx_{M} ⟩ + ⟨ d, x_{M} ⟩,

and thus condition (4.35) is equivalent (for vectors ū satisfying (4.32)) to

⟨ ū, Cx̄_{M} ⟩ + ⟨ d, x̄_{M} ⟩ = min_{x_{M} ∈M}⟨ ū, Cx_{M} ⟩ + ⟨ d, x_{M} ⟩,  (4.36)
⟨ ū, Cx̄_{E} ⟩ + ⟨ d, x̄_{E} ⟩ = 0.

Since (4.36) is equivalent to

C^{⊤}ū + d ∈ −N_{M}(x̄_{M}),

the result follows. □

Note that the conditions derived in [23] and [24] for the single-facility location model (see Subsection 4.1.2) are special cases of the ones derived in the last theorem.

Remark 4.4.5 (Dual Gap) Let x = x_{M} +x_{E} (x_{M} ∈ M, x_{E} ∈ E) be primal feasible and let u be dual feasible. Additionally, let γ(Cx + c) = ⟨ u, Cx + c ⟩ and let δ_{M}^{∗}(−C^{⊤}u − d) = −⟨ x_{M}, C^{⊤}u + d ⟩. (Note that these two additional conditions are equivalent to u ∈ ∂γ(Cx + c) and x_{M} ∈ ∂δ_{M}^{∗}(−C^{⊤}u − d).) A simple calculation then shows that ⟨ x_{E}, C^{⊤}u + d ⟩ is the dual gap with respect to the feasible points x and u.

Remark 4.4.6 If G_{∞} is a linear subspace then G_{∞}^{∗} = G_{∞}^{⊥}, and thus the complementary condition ⟨ C\overline{u} + d, \overline{x}_{E} ⟩ = 0 in Theorem 4.4.4 is redundant.

We have shown in Theorem 4.4.3 that primal and dual optimal solutions (when they exist!) are related with each other as saddle point solutions of (4.31). However, the existence of optimal solutions for (4.17) is not guaranteed when (4.17) has a finite optimal value. Since in applications E = {0} usually does not hold (see Section 4.1), a deeper analysis is required. This is the purpose of the rest of the section.

For certain instances of (4.17), the non-emptiness and compactness of the set of optimal primal solutions can be derived by ad-hoc procedures, as done, e.g., in [10, 14, 70]. For the general situation we have the following.

Theorem 4.4.7 (Bounds on Primal Solutions) If all nonzero y ∈ G_{∞} satisfy γ(Cy) + ⟨ d, y ⟩ > 0, then

1. The set of optimal solutions is nonempty, convex, and compact.

2. Let (M, E) be an a. c. r. of G and suppose E ≠ {0}. Let

L ≥ max_{x∈M} (γ(−Cx − c) − ⟨ d, x ⟩) ,
\overline{v} ≥ min_{x∈G} (γ(Cx + c) + ⟨ d, x ⟩) ,
r ∈ ]0, min_{e∈E,∥e∥_{2}=1} γ(Ce) + ⟨ d, e ⟩] .

Then, any optimal solution x_{M} +x_{E} for (4.17) with x_{M} ∈ M and x_{E} ∈ E satisfies

∥x_{E}∥_{2} ≤ \frac{\overline{v} + L}{r}

Proof: Part 1 follows from the fact that, under these assumptions, the function g + δ_{K} is coercive and thus level bounded. See, e. g., [4].

To show Part 2, let x_{E} ∈ E be given with ∥x_{E}∥_{2} > (\overline{v} + L)/r and let x_{M} ∈ M. The triangle inequality now shows that

γ(Cx_{M} + Cx_{E} + c) + ⟨ d, x_{M} ⟩ + ⟨ d, x_{E} ⟩ ≥
≥ γ(Cx_{E}) + ⟨ d, x_{E} ⟩ − γ(−Cx_{M} − c) + ⟨ d, x_{M} ⟩
= ∥x_{E}∥_{2} (γ (\frac{Cx_{E}}{∥x_{E}∥_{2}} ) + ⟨d, \frac{x_{E}}{∥x_{E}∥_{2}} ⟩) + ⟨ d, x_{M} ⟩ − γ(−Cx_{M} − c)
> \frac{\overline{v} + L}{r} (γ (\frac{Cx_{E}}{∥x_{E}∥_{2}} ) + ⟨d, \frac{x_{E}}{∥x_{E}∥_{2}} ⟩) − L
≥ \overline{v},

contradicting the optimality of x_{M} + x_{E}. □

Remark 4.4.8 In Theorem 1 of [59], it is assumed that C is a p × q matrix (q < p) with rank q, G = IR^{q}, and 0 ∈ C^{⊤}int(B^{∘}) + d. This is clearly stronger than the assumption in Theorem 4.4.7. Indeed, for any e ∈ G_{∞}, one of the two following conditions hold:

max_{u∈B^{∘}}⟨ e, C^{⊤}u + d ⟩ > 0,  (4.37)
⟨ e, C^{⊤}u + d ⟩ = 0  ∀ u ∈ B^{∘}.  (4.38)

If (4.37) holds, then

γ(Ce) + ⟨ d, e ⟩ > 0,

whereas if (4.38) holds, then

d + C^{⊤}B^{∘} ∈ {e}^{⊥},

which, since C^{⊤}B^{∘} has full dimension, implies that e = 0. Hence, γ(Ce) + ⟨ d, e ⟩ > 0 for all nonzero e ∈ G_{∞}. □

When the assumption in Theorem 4.4.7 is invalid, non-emptiness of the set of optimal solutions cannot be guaranteed in general, even for the case of polyhedral recession cone G_{∞}. This is shown in the following example.

Example 4.4.1 Let C = (1, 0)^{⊤}, d = (−1), c = (0, 1)^{⊤}, the feasible set G = [0, +∞[= G_{∞}, and let γ be the Euclidean norm in the plane. Then,

γ(Cx + c) + ⟨ d, x ⟩ = \sqrt{x^{2} + 1} − x,

which is always non-negative, but tends to zero when x grows to infinity. Hence, no optimal solution exists. □

The case d = 0 (in fact the common one in applications) simplifies the analysis since then the objective function of (4.17) is bounded below. However, this does not guarantee the attainment of the optimal value, as shown in the following example.

Example 4.4.2 Let G = {(x_{1}, x_{2}, x_{3}) : x_{3}^{2} ≥ x_{1}^{2} + x_{2}^{2}, x_{i} ≥ 0, i = 1, 2, 3} and let C be the matrix

( \begin{matrix}1 & 0 & −1 \\ 0 & −1 & 0 \end{matrix} ) .

Let c = (0, 1)^{⊤}, d = 0, and let γ be the Euclidean norm. Then,

γ(Cx + c) + ⟨ d, x ⟩ = γ(x_{1} − x_{3}, 1 − x_{2}).

Since the system

\{ \begin{matrix}x_{1} & = & x_{3} \\ x_{2} & = & 1 \end{matrix}

has no solution on G, the objective function is strictly positive on G. However, for the feasible sequence ((k, 1, \sqrt{1 + k^{2}})^{⊤})_{k∈IN} the objective value tends to zero, showing that the infimum (zero) is not attained. □

In spite of this negative result, a geometrical condition can be given to guarantee the attainment of the optimal value for d = 0:

Theorem 4.4.9 (Attainment of Optimal Values) Let d = 0. The following conditions are equivalent.

1. (4.17) attains its optimal value for each c ∈ IR^{n}.

2. The set CG is closed.

Proof: Assume Condition 1 and suppose CG is not closed. Then there exist a vector v^{∗} and a sequence (x_{j})_{j∈IN} with x_{j} ∈ G for all j ∈ IN such that

lim_{j→+∞} Cx_{j} = v^{∗} ∉ CG.  (4.39)

However, this would imply that, for c = −v^{∗}, the objective function value is always strictly positive but converging to zero for the feasible sequence (x_{j})_{j∈IN}. This contradicts the assumption.

Conversely, if Condition 2 holds, then, formulating (4.17) as the problem

minimize γ(y)
subject to y ∈ CG + c

one immediately obtains that (4.17) amounts to finding the point in the closed set CG + c closest to the origin with respect to the gauge γ, which always admits an optimal solution.

□ Corollary 4.4.10 Let G_{1}, G_{2} be asymptotically conical sets such that inf{γ(x − y) | x ∈ G_{1}, y ∈ G_{2}} = 0. Then, G_{1} ∩ G_{2} ≠ ∅.

Proof: The set G := G_{1} × G_{2} is asymptotically conical, see Property 4.2.2. Let C = (I_{n}, −I_{n}), where I_{n} is the n × n identity matrix. Then CG = G_{1} − G_{2}, which is asymptotically conical as well, see Property 4.2.2, thus it is closed. By Theorem 4.4.9, the problem

inf_{x∈G} γ(Cx)

attains its infimal value (zero), which means that some x exists in G_{1} ∩ G_{2}. □

4.5 Solving Goal Programs Efficiently

The aim of this section is to show how the structure of (4.17) can be exploited to derive polynomial time interior-point schemes for solving the problem at hand. Note that we are not considering warm-start techniques to compute several (or all) efficient points of a given multicriteria optimization problem: this was the purpose of Chapter 3. When the number of criteria is large, the strategy outlined in Chapter 3 might become to unwieldy, and a decision-maker might be contempt with the construction of just a few efficient points. In this case, the technique outlined in this section might be useful.

The discussion below will show that particular instances of the general problem can be solved by interior-point methods of various types (e. g. primal, primal-dual, short-step, long-step, etc.), provided that self-concordant barriers for the unit ball of the gauge γ and M as well as for the cone E are given.

For instance, consider the primal problem (4.17) restated in the form

minimize s + t,
subject to (Cx + c, s) ∈ epi(γ),
⟨ d, x ⟩ ≤ t,
x ∈ G

and let B be the unit ball of the gauge γ. Then epi(γ) is the conic hull of this unit ball, i. e.

epi(γ) = {(y, λ) ∈ IR^{n+1} | y ∈ λB; λ ≥ 0}.

The reformulated problem has therefore a linear objective function, a conic constraint, and a convex constraint of a rather special structure, which makes it easily exploitable for interior-point methods.

As a simple example, we might use the standard primal path-following algorithm from Nesterov and Nemirovskii [67], already outlined in short in Chapter 3. For this, we need not only a starting point in the interior of a compact set of feasible points, but also a lower bound on the asymmetry coefficient of this starting point. Recall from Chapter 3 that the asymmetry coefficient a(x, G) of a point x lying in the interior of the convex compact set G is defined by

a(x, G) := sup{α ≥ 0 | x + α(x − G) ⊆ G}.

Denoting by γ_{G−x} the gauge with unit ball G − x, one immediately obtains from the definition that

a(x, G) = sup{1/β > 0 | (x − G) ⊆ β(G − x)}
= (inf{β > 0 | (x − G) ⊆ β(G − x)})^{−1}
= (inf{β > 0 | ∀ y ∈ G − x : −y ∈ β(G − x)})^{−1}
= (inf{β > 0 | ∀ y ∈ G − x : γ_{G−x}(−y) ≤ β})^{−1}
= ( max_{y∈G−x} γ_{G−x}(−y))^{−1} .  (4.40)

We have the following estimate for a point x ∈ int(G).

Lemma 4.5.1 Let G ⊆ IR^{n} be a convex and compact set, let x ∈ int(G) and let B_{1}, B_{2} ⊆ IR^{n} be convex compact sets with the origin in their interior, satisfying

B_{1} ⊆ G − x ⊆ B_{2}.

Denote by γ_{B_{1}} the gauge with unit ball B_{1}. Then,

a(x, G) ≥ (max_{y∈B_{2}} γ_{B_{1}}(−y))^{−1}

Proof: By assumption and the definition of gauges,

γ_{G−x} ≤ γ_{B_{1}}.

Hence, it follows from (4.40) that

a(x, G) ≥ ( max_{y∈G−x} γ_{B_{1}}(−y))^{−1}
≥ (max_{y∈B_{2}} γ_{B_{1}}(−y))^{−1} .
□

Corollary 4.5.2 Let G ⊂ IR^{n} be a convex compact set and let x ∈ int(G). Let B_{1} and B_{2} be the unit balls of the ∥ · ∥_{1}- and ∥ · ∥_{2}-norm, respectively. Moreover, let r_{1}, r_{2} > 0 be two constants such that x + r_{1}B_{1} ⊆ G and G ⊆ x + r_{2}B_{2}. Then

a(x, G) ≥ \frac{r_{1}}{\sqrt{n}r_{2}} .

Proof: By Lemma 4.5.1,

a(x, G) ≥ ( max_{y∈r_{2}B_{2}} γ_{r_{1}B_{1}}(−y))^{−1}
= \frac{r_{1}}{r_{2}} (max_{y∈B_{2}} γ_{B_{1}}(−y))^{−1}
= \frac{r_{1}}{r_{2}\sqrt{n}}.
□

Suppose now that we are given a self-concordant barrier b_{B} for the unit ball B with self-concordancy parameter ϑ_{B} ≥ 1 and a self-concordant barrier b_{G} for the set of feasible points G with self-concordancy parameter ϑ_{G} ≥ 1. If an a. c. r. G = M + E for G is known, the latter barrier will usually be written as b_{G} = b_{M} + b_{E}, where b_{M} is a barrier for the compact set M, while b_{E} is the corresponding barrier for the cone E.

Theorem 4 from [41] tells us that it is relatively easy to construct a self-concordant barrier for the epigraph of γ explicitly. Indeed, such a barrier takes the form

b^{+}(x, t) = βb_{B}(x/t) − αϑ_{B} ln t,  (4.41)

where α, β > 0 are explicitly given constants, depending only on ϑ_{B}. Both constants are of magnitude O(1). Moreover, b^{+} has a self-concordancy parameter of order O(ϑ_{B}).

Let there be given a point y ∈ int(G). In the next step, we have to assume that we know a bound r̃ > 0 such that for every solution x ∈ G of our primal problem the relation ∥x∥_{2} ≤ r̃ holds. See Theorem 4.4.7 for methods for constructing such r̃. Moreover, we assume that r̃ is chosen in such a way that ∥y∥_{2} + 1 ≤ r̃. Define now

ϱ_{P} := γ(Cy + c) + ⟨ d, y ⟩ + 3

and

G_{P} := {(x, s, t) ∈ IR^{m+2} | x ∈ G, ∥x∥_{2} ≤ r̃, (Cx + c, s) ∈ epi(γ),
⟨ d, x ⟩ ≤ t, s + t ≤ ϱ_{P} }.

Obviously, G_{P} is a convex compact set. With u := γ(Cy + c) + 1 and v := ⟨ d, y ⟩ + 1 we have that the point ŷ:= (y, u, v) is in the strict interior of G_{P}. Moreover,

b_{P} (x, t) := b^{+}(Cx + c, t) + b_{K}(x) − ln(r̃^{2} − ∥x∥_{2}^{2}) − ln(t − ⟨ d, x ⟩) − ln(ϱ_{P} − s − t)

is a self-concordant barrier for G_{P} with self-concordancy parameter

ϑ_{P} := O(1)ϑ_{B} + ϑ_{K} + 3 = O(ϑ_{B} + ϑ_{K})

(see [67, Proposition 5.1.1 and 5.1.2]). This means that we can opt to solve the problem

minimize s + t
subject to (x, s, t) ∈ G_{P}

with an interior-point method, using ŷ as a starting point.

Lemma 4.5.3 (Bounds on the Asymmetry Coefficient) Denote by e_{i} be the ith Euclidean unit vector (i = 1, …, m) in the IR^{m}. Let δ_{G} > 0 be such that y+δ_{G}e_{i}, y−δ_{G}e_{i} ∈ G for all i = 1, …, m, and denote by c_{i} ∈ IR^{n}, i = 1, …, m, the columns of the matrix C. Define

r_{P} := min{1, δ_{G}, 1/γ(c_{i}), 1/γ(−c_{i}), 1/|d_{i}| : i = 1, …, n}

and

R := (r̃^{2} + (ϱ_{P} + r̃∥d∥_{2})^{2} + (max{|ϱ_{P} |, r̃∥d∥_{2}})^{2})^{1/2} .

It then follows that

a(ŷ, G_{P}) ≥ \frac{r_{P}}{2\sqrt{n + 2}R}.

Proof: First, let (x, s, t) ∈ G_{P}. It then follows that ∥x∥_{2} ≤ r̃, 0 ≤ s, and −r̃∥d∥_{2} ≤ t ≤ ϱ_{P}. As a consequence, we have s ≤ ϱ_{P} + r̃∥d∥_{2}. This means that

∥(x, s, t)∥_{2}^{2} = ∥x∥_{2}^{2} + s^{2} + t^{2} ≤ r̃^{2} + (ϱ_{P} + r̃∥d∥_{2})^{2} + (max{|ϱ_{P} |, r̃∥d∥_{2}})^{2}

Therefore,

G_{P} ⊆ ŷ + 2RB_{2},

where B_{2} is the unit ball of the ∥ · ∥_{2}-norm in IR^{m+2}.

Second, with the ith Euclidean unit vector e_{i} ∈ IR^{m} we have for all z ∈ IR^{m} that γ(C(z + r_{P} e_{i}) + c) ≤ γ(Cz + c) + 1 and γ(C(z − r_{P} e_{i}) + c) ≤ γ(Cz + c) + 1 for i = 1, …, m. Moreover, ⟨ d, y +r_{P} e_{i} ⟩ ≤ v as well as ⟨ d, y −r_{P} e_{i} ⟩ ≤ v (i = 1, …, m). This means that ŷ + r_{P} B_{1} ⊆ G_{P}, where B_{1} is the unit ball of the ∥ · ∥_{1}-norm in IR^{m+2}. The result follows with Corollary 4.5.3. □

With the last lemma, it is easy to see that Stage 1 of the standard primal path-following algorithm from [67] takes

O (\sqrt{ϑ_{P}} (ln ϑ_{P} + ln n + ln R − ln r_{P}))

iterations, while Stage 2 of this method takes

O (\sqrt{ϑ_{P}} (ln ϑ_{P} + ln(1/ε) + ln(ϱ_{P} + 2r̃∥d∥_{2})))

iterations to achieve ε-accuracy.

Bounding r_{P}, ϑ_{B}, and ϑ_{G} depends, of course, on the nature of the actual data at hand. We will consider typical examples of goal programming problems in Section 4.5.1.

Note that the dual problem in the formulation (4.23) can be treated in the same way as the primal one, provided that there are known self-concordant barriers for M^{∘}, B^{∘}, and E^{∗}. A latter one can, at least in principle, be constructed from a self-concordant barrier for E, see Section 2.4.1 in [67]. One just has to change the maximization problem into an equivalent minimization problem and add one slack variable s ∈ IR for the constraint (−C^{⊤}u − d, s) ∈ epi(γ̃^{∘}).

Another possibility is to consider a primal-dual reformulation of the problem. Assume that the set M has non-empty interior, and let x_{0} ∈ int(M) be given. Define the gauge γ̃ by its unit ball: γ̃ := γ_{M−x_{0}} . Example 4.3.1 now shows that

minimize γ (Cx + c) + γ̃^{∘} (−C^{⊤}u − d) + ⟨ d, x ⟩ − ⟨ Cx_{0} + c, u ⟩ − ⟨ d, x_{0} ⟩
subject to x ∈ M + E,
C^{⊤}u + d ∈ E^{∗}
γ^{∘}(u) ≤ 1.

is the primal-dual reformulation of problem (4.17) with objective function value 0. Note that this again is a problem of the same type as (4.17).

Other algorithms, especially long-step methods, can be derived when more knowledge is available about the problem structure. As a trivial example, if K as well as B is polyhedral, the problem reduces to a linear one, for which methods of higher efficiency then the one depicted above are readily available. (See also Examples 4.5.2, 4.5.3, and 4.5.4 in the next subsection.) Other possibilities include the cases in which the cones E and epi(γ) appearing in the problem formulation are direct sums of cones of positive semidefinite symmetric real matrices and cones of the form epi(∥ · ∥_{2}). (See Example 4.5.1 in the next subsection.) The standard barriers for these cones are self-scaled, allowing for especially efficient algorithms, cmp. [68]. Note, however, that the use of a p-norm with p ≠ 2 (an important case in applications [72]) does not allow for a self-scaled cone, and that a self-dual formulation for the corresponding problem is not readily at hand. Indeed, interior-point methods proposed up to now for this class of problems do not use a self-dual formulation, see [96] and Example 4.5.5 in the next subsection.

4.5.1 Particular Cases

In this section we take a quick look at how self-concordant barriers for the unit balls of typical gauges encountered in applications can be easily derived. Some of these cases have already been discussed in [33], in the context of interior point algorithms applied to specific location problems similar to the one discussed in Subsection 4.1.2.2.

Example 4.5.1 (Euclidean Norm) (cp. Proposition 5.4.2 and 5.4.3 in [67]) If γ is the Euclidean norm, a self-concordant barrier with self-concordancy parameter ϑ_{B} = 1 for the unit ball B of γ is given by b_{B}(x) := − ln(1 − ∥x∥_{2}^{2}). Likewise, if the unit ball of γ is an ellipsoid, B = {x ∈ IR^{n} | ⟨ x, Qx ⟩ ≤ 1}, where Q is a positive definite n × n-matrix, a self-concordant barrier with self-concordancy parameter ϑ_{B} = 1 is given by b_{B}(x) := − ln(1 − ⟨ x, Qx ⟩). In both cases, the tedious general construction of a selfconcordant barrier for epi(γ) can be avoided by noting that b^{+}(x, t) = − ln(t^{2}−⟨ x, Qx ⟩) is a self-concordant barrier for this cone with self-concordancy parameter ϑ_{B}^{+} = 2. Note that epi(∥ · ∥_{2}) is just the standard second-order cone, while b^{+} is the corresponding self-scaled barrier for this cone, see [68]. □

Example 4.5.2 (Polyhedral Gauges) Let γ be a polyhedral gauge whose unit ball is given by a set of k linear inequalities: B = {x ∈ IR^{n} | Ax ≤ g}, A ∈ IR^{k×n}, g ∈ IR^{k}. Of course, the standard logarithmic barrier b_{B}(x) = − ∑_{i=1}^{k} ln(g_{i} −a_{i}^{⊤}x) for the polytope B can be used to define b_{B}^{+}(x, t) = − ∑_{i=1}^{k} ln(g_{i}t − a_{i}^{⊤}x), a self-concordant barrier for the epigraph of γ with self-concordancy parameter ϑ_{B}^{+} = k. □

Example 4.5.3 (Skewed Gauges) Let γ be a gauge, A ∈ IR^{n×n} be a regular matrix and c ∈ IR^{n} be a vector with γ^{∘}(A^{⊤}c) < 1. Then γ̃(x) := γ(Ax) + ⟨ c, x ⟩ defines a gauge [71]. Gauges defined like this have important applications in location science, see, e. g., [72, 36], and are called skewed gauges. It is easy to see that the unit ball of γ̃ is given by (c + A^{⊤}(B^{∘}))^{∘}. However, finding a barrier for this unit ball does not seem to be so easy. On the other hand, finding a self-concordant barrier for epi(γ̃) is simple, as long as such a barrier for the unit ball B of γ is given. Let b be such a barrier with self-concordancy parameter ϑ_{B} and let b^{+} be a barrier for epi(γ) with self-concordancy parameter ϑ_{B}^{+} (note that we have ϑ_{B}^{+} = O(ϑ_{B}), according to [41]). Then,

b̃^{+}(x, t) := b^{+}(Ax, t − ⟨ c, x ⟩) − ln t

is a barrier for the epigraph of γ̃ with self-concordancy parameter ϑ̃^{+} = ϑ_{B}^{+} + 1. This is a simple application of Proposition 5.1.1 and 5.2.5 from [67]. □

Example 4.5.4 Suppose that we are given k gauges γ_{i} (i = 1, …, k) and we want to use the gauge γ defined by γ(x) := max_{i=1}^{k} γ_{i}(x). If B_{i} is the unit ball of γ_{i} (i = 1, …, k), it is easy to see that B = ⋂_{i=1}^{k} B_{i} is the unit ball of γ. Given self-concordant barriers b_{i} for the unit balls B_{i} with self-concordancy parameter ϑ_{i}, we have that b_{B} := ∑_{i=1}^{k} b_{i} is a self-concordant barrier for B with self-concordancy parameter ϑ_{B} = ∑_{i=1}^{k} ϑ_{i}. □

Example 4.5.5 (p-Norms) We consider the case that γ is a p-norm with p ∈]1, ∞[. By introducing slack variables for the inequalities describing the unit ball B of γ, we can consider the set

B̂_{p} := {(x, y, z) ∈ IR^{3n} | −y_{i} ≤ x_{i} ≤ y_{i}, 0 ≤ y_{i}, y_{i}^{p} ≤ z_{i} (i = 1, …, n),

∑_{i=1}^{n} z_{i} ≤ 1}.

Using Proposition 5.3.1 from [67], we see that

b̂(x, y, z) := − ∑_{i=1}^{n} (ln(y_{i} − x_{i}) + ln(y_{i} + x_{i}) + ln y_{i} + ln z_{i} + ln (z_{i}^{1/p} − y_{i}))

− ln (1 − ∑_{i=1}^{n} z_{i})

is self-concordant barrier for B̂ with self-concordancy parameter 6n + 1. Moreover, (4.41) tells us how to construct a self-concordant barrier with self-concordancy parameter ϑ^{+} = O(n) for the conic hull

K(B̂_{p}) = {(x, y, z, t) ∈ IR^{3n+1} | (x, y, z) ∈ tB̂, t ≥ 0}
= {(x, y, z, t) | −y_{i} ≤ x_{i} ≤ y_{i}, 0 ≤ y_{i}, y_{i}^{p} ≤ z_{i}t^{p−1} (i = 1, …, n),

∑_{i=1}^{n} z_{i} ≤ t, 0 ≤ t}.

This set can be used to replace the epigraph of γ in the definition of G_{P}, see also [33].

Note that, using the result of [41], it is easy to refine the bound on ϑ^{+} and show that ϑ^{+} ≤ 86.43n+14.41 holds. In contrast to this, Xue and Ye [96] use different embeddings to derive barriers for the corresponding conic hulls with self-concordancy parameters of 200n and 200p^{2}n + 100n + 200, respectively.

The construction of a starting point lying in the strict interior of the set of feasible points and the estimate of the asymmetry of this point can be done as shown in Section 4.5 and is discussed in more detail in [33].

(The rest of this example has been developed in collaboration with Yinyu Ye from the University of Iowa.) Note that, depending on the formulation one is using, either the dual cone (K(B̂_{p}))^{∗} or K(B̂_{q}) (1/p + 1/q = 1) will appear in the dual problem. But with the embedding above, (K(B̂_{p}))^{∗} ≠ K(B̂_{q}) holds. Therefore, an explicit construction of (K(B̂_{p}))^{∗} is in order. To fix notation, recall that

(K(B̂_{p}))^{∗} = {(α, β, γ, τ) ∈ IR^{3n+1} |⟨(α, β, γ, τ), (x, y, z, t)⟩ ≥ 0 ∀(x, y, z, t) ∈ K(B̂_{p})}.

As a first step, write

K(B̂_{p}) = K̂_{p}^{1} ∩ K̂_{p}^{2} ∩ K̂_{p}^{3}

with

K̂_{p}^{1} := {(x, y, z, t) ∈ IR^{3n+1} | −y ≤ x ≤ y, y ≥ 0},
K̂_{p}^{2} := {(x, y, z, t) ∈ IR^{3n+1} | ∑_{i=1}^{n} z_{i} ≤ t} ,
K̂_{p}^{3} := {(x, y, z, t) ∈ IR^{3n+1} | y ≥ 0, z ≥ 0, t ≥ 0, y_{i}^{p} ≤ z_{i}t^{p−1} (i = 1, …, n)}.

Using [77, Corollary 16.4.2] we have that

(K(B̂_{p}))^{∗} = (K̂_{p}^{1})^{∗} + (K̂_{p}^{2})^{∗} + (K̂_{p}^{3})^{∗},

so we just have to compute the cones on the right-hand side. First,

(K̂_{p}^{1})^{∗} = {(α, β, γ, τ) ∈ IR^{3n+1} | ∀(x, y, z, t) ∈ K̂_{p}^{1} : ⟨(α, β, γ, τ), (x, y, z, t)⟩ ≥ 0}
= {(α, β, γ, τ) ∈ IR^{3n+1} | ∀(x, y, z, t) ∈ K̂_{p}^{1} :
⟨ α, x ⟩ + ⟨ β, y ⟩ + ⟨ γ, z ⟩ + τt ≥ 0}
= {(α, β, 0, 0) ∈ IR^{3n+1} | ∀(x, y, z, t) ∈ K̂_{p}^{1} : ⟨ α, x ⟩ + ⟨ β, y ⟩ ≥ 0}
= {(α, β, 0, 0) ∈ IR^{3n+1} | β ≥ 0, |α_{i}| ≤ β_{i} (i = 1, …, n)}.

Second, since K̂_{p}^{2} is a halfspace,

(K̂_{p}^{2})^{∗} = {(0, 0, γ, τ) ∈ IR^{3n+1} | 0 ≤ τ = −γ_{1} = ⋯ = −γ_{n}}
= {τ(0, 0, −e, 1) ∈ IR^{3n+1} | τ ≥ 0}

(with e = (1, …, 1) ∈ IR^{n}). Third, let there be given (α, β, γ, τ) ∈ (K̂_{p}^{3})^{∗}, i. e. ⟨ α, x ⟩ + ⟨ β, y ⟩ + ⟨ γ, z ⟩ + τt ≥ 0 for all (x, y, z, t) ∈ K̂_{p}^{3}. Since we can choose x independently (and arbitrarily) of y, z, t we immediately get α = 0.

Let (x, y, z, t) ∈ K̂_{p}^{3}. If t = 0 we have y = 0, and

⟨ β, y ⟩ + ⟨ γ, z ⟩ + τt ≥ 0  (4.42)

holds as long as γ ≥ 0 (note z ≥ 0). Now consider values of t with t > 0. Under this assumption, the inequality (4.42) is equivalent to

\frac{1}{t}⟨ β, y ⟩ + \frac{1}{t}⟨ γ, z ⟩ + τ ≥ 0.

Moreover, y_{i}^{p} ≤ z_{i}t^{p−1} is equivalent to (y_{i}/t)^{p} ≤ z_{i}/t. This allows a change of variables according to ŷ := (1/t)y and ẑ := (1/t)z. With this, (0, β, γ, τ) ∈ (K̂_{p}^{3})^{∗} is equivalent to

⟨ β, ŷ⟩ + ⟨ γ, ẑ⟩ + τ ≥ 0

for all ŷ, ẑ ≥ 0 with ŷ_{i}^{p} ≤ ẑ_{i} (i = 1, …, n). As a consequence, we have necessarily γ ≥ 0 (take ẑ_{i} large) and τ ≥ 0 (take ŷ = ẑ = 0). Moreover,

⟨ β, ŷ⟩ + ⟨ γ, ẑ⟩ + τ = ∑_{i=1}^{n}(β_{i}ŷ_{i} + γẑ_{i}) + τ,

and since γ ≥ 0 we conclude that only the case ẑ_{i} = ŷ_{i}^{p} (i = 1, …, n) needs to be considered. But then

⟨ β, ŷ⟩ + ⟨ γ, ẑ⟩ + τ = ∑_{i=1}^{n} f_{i}(ŷ_{i}),

where f_{i}(ŷ_{i}) := β_{i}ŷ_{i} + γ_{i}ŷ_{i}^{p} + τ/n. Since only the constraints ŷ ≥ 0 have to be fulfilled, we can minimize each function f_{i} separately. If γ_{i} = 0, we have f_{i}(ŷ_{i}) ≥ 0 for all ŷ_{i} ≥ 0 if and only if β_{i} ≥ 0. Consider now the case γ_{i} > 0. If β_{i} ≥ 0, then f_{i}(ŷ_{i}) ≥ 0 for all ŷ_{i} ≥ 0 follows. Suppose therefore that β_{i} < 0 and γ_{i} > 0. Then the function f_{i} takes on its minimal value at

ŷ_{i}^{∗} := (\frac{−β_{i}}{pγ_{i}} )^{\frac{1}{p−1}} ,

and we arrive at

f_{i}(ŷ_{i}^{∗}) = β_{i}ŷ_{i}^{∗} + γ_{i}(ŷ_{i}^{∗})^{p} + τ/n
= −|β_{i}|^{1+ \frac{1}{p−1}} p^{−1/(p−1)}γ_{i}^{−1/(p−1)} + |β_{i}|^{p/(p−1)}p^{−p/(p−1)}γ_{i}^{1− \frac{p}{p−1}} + τ/n
= −|β_{i}|^{p/(p−1)}γ_{i}^{−1/(p−1)} (p^{−1/(p−1)} − p^{−p/(p−1)}) + τ/n.

Define now q ∈]1, +∞[ by 1/p + 1/q = 1 and

c_{p} := p^{−1/(p−1)} − p^{−p/(p−1)} = p^{−q/p} − p^{−q} = p^{−(q−1)} − p^{−q} = p^{−q}(p − 1) = \frac{p^{1−q}}{q} .

Then f_{i}(ŷ_{i}^{∗}) ≥ 0 if and only if

nc_{p}|β_{i}|^{q} ≤ τγ_{i}^{q/p}.

With q/p = q − 1, the inequality above is equivalent to

nc_{p}|β_{i}|^{q} ≤ τγ_{i}^{q−1}

i. e.

(nc_{p})^{p/q}|β_{i}|^{p} ≤ τ^{p−1}γ_{i}

(use p/q = p − 1). In this formulation, β plays the same role as the primal variable y, while γ takes over the role of z. Of course, τ has the same meaning as t.

As a consequence of the discussion above, we arrive at

(K̂_{p}^{3})^{∗} = {(0, β, γ, τ) ∈ IR^{3n+1} | γ ≥ 0, τ ≥ 0,
(nc_{p})^{p−1}(max(−β_{i}, 0))^{p} ≤ γ_{i}τ^{p−1} (i = 1, …, n)}
= {(0, β, γ, τ) ∈ IR^{3n+1} | γ ≥ 0, τ ≥ 0,
(nc_{p})^{1/q}(max(−β_{i}, 0)) ≤ γ_{i}^{1/p}τ^{1/q} (i = 1, …, n)}
= {(0, β, γ, τ) ∈ IR^{3n+1} | γ ≥ 0, τ ≥ 0,
nc_{p}(max(−β_{i}, 0))^{q} ≤ γ_{i}^{q−1}τ (i = 1, …, n)}

which is the cone induced by the set

B̂_{p}^{3}{(0, β, γ) ∈ IR^{3n} | (nc_{p})^{p−1}(max(−β_{i}, 0))^{p} ≤ γ_{i} (i = 1, …, n)}.

Taking the results for (K̂_{p}^{1})^{∗}, (K̂_{p}^{2})^{∗}, and (K̂_{p}^{3})^{∗} together, we get

(K(B̂_{p}))^{∗} = \{( \begin{matrix}α^{[1]} \\ β^{[1]} + β^{[3]} \\ −τ^{[2]}e + γ^{[3]} \\ τ^{[2]} + τ^{[3]} \end{matrix} ) ∈ IR^{3n+1} |\begin{matrix}−β^{[1]} ≤ α^{[1]} ≤ β^{[1]}, \\ γ^{[3]} ≥ 0, τ^{[2]} ≥ 0, τ^{[3]} ≥ 0, \\ (nc_{p})^{p−1}(max(−β^{[3]}, 0))^{p} ≤ (τ^{[3]})^{p−1}γ^{[3]} \end{matrix} \} ,

where the exponents, the max-operator as well as the inequalities are understood to be taken componentwise.

This representation of the dual cone is slightly less economical than the representation of the primal one. More precisely, 4n + 2 real-valued variables are needed instead of 3n + 1. Moreover, another n variables are used when the max-operator is resolved.

It seems also that there is no symmetry at all in this formulation. Take p = 2. Then K(B̂_{2}) ≠ (K(B̂_{2}))^{∗}. (Just use α^{[1]} = β^{[1]} = β^{[3]} = γ^{[3]} = 0. Then (0, 0, −τ^{[1]}e, τ^{[1]} + τ^{[3]}) ∈ (K(B̂_{2}))^{∗}, in contrast to z ≥ 0 if (x, y, z, t) ∈ K(B̂_{2}).)

But fortunately, an explicit description of a self-concordant barrier for (K(B̂_{p}))^{∗} poses no problem. By adding variables δ ∈ IR^{n} to get rid of the max-operator, we can describe the nonlinear constraints defining B_{p}^{3} by the barrier function

b_{3}(β, δ, γ) := − ∑_{i=1}^{n}(ln δ_{i} + ln(δ_{i} + β_{i}) + ln γ_{i} + ln(γ_{i}^{1/p} − (nc_{p})^{1/q}δ_{i})),

which is self-concordant with self-concordancy parameter 5n. From this, we can build a barrier for the cone K̂_{p}^{3} with self-concordancy parameter O(n). All other constraints in the description of (K(B̂_{p}))^{∗} are linear and and can be taken care of in the usual way. □

Example 4.5.6 (Composite Gauges) Let γ be a gauge as in (4.24). Using the fact that γ̃ is monotonic, it is sufficient to consider the set

B̂ = {(u_{1}, t_{1}, …, u_{k}, t_{k}) ∈ IR^{n_{1}+1} × ⋯ × IR^{n_{k}+1} | γ_{i}(u_{i}) ≤ t_{i} (i = 1, …, k),
γ̃(t_{1}, …, t_{k}) ≤ 1}.

Obviously, to construct a self-concordant barrier for the set B̃, one can use self-concor-dant barriers b_{i}^{+} with self-concordancy parameter ϑ_{i}^{+} for the cones epi(γ_{i}) and a selfconcordant barrier b̃ with self-concordancy parameter ϑ̃ for the unit ball of γ̃ to define

b̂(u_{1}, t_{1}, …, u_{k}, t_{k}) := b̃(t_{1}, …, t_{k}) + ∑_{i=1}^{k} b_{i}^{+}(u_{i}, t_{i}),

a self-concordant barrier for B̂ with self-concordancy parameter ϑ̃+∑_{i=1}^{k} ϑ_{i}^{+}. Moreover, Lemma 4.3.1 can be used to construct barriers for the dual gauge γ^{∘} and its unit ball. □

Nonparametric Descent Methods for Nonconvex Multicriteria Problems

In case the number of criteria of a given multicriteria optimization problem is large or the criteria functions are nonconvex, the strategy outlined in Chapter 2 and Chapter 3 might not be feasible. Moreover, if convexity is missing, the duality results from Chapter 4 do not hold, and consequently all that we can hope for is to efficiently compute locally efficient points. It is therefore important that we are able to compute such points for problems of the type mentioned above.

In this chapter, we propose a parameter-free optimization strategy for computing a point satisfying the first-order necessary conditions for multicriteria optimization. This strategy generalizes descent direction methods from single-criterion optimization, methods known for their computationally efficiency when solving large-scale nonconvex problems. Neither ordering information nor weighting factors or similar value function parameters for the different objective functions given are assumed to be known. Moreover, no value function will be used to scalarize the multicriteria problem at hand.

5.1 The Problem Considered

As in Chapter 2, let K ⊆ IR^{n} be a cone and let the IR^{n} be equipped with the ordering ≤_{K} induced by K. Suppose that we are in search for a (local) minimum of a given function

F : IR^{m} ⟶ IR^{n},

not necessarily convex. A (local) minimum x ∈ IR^{m} of F corresponds to a (local) efficient point F (x) ∈ IR^{n} of the manifold F (IR^{m}). Constrained problems will be considered in Section 5.7. In what follows, we assume that the function F is continuously differentiable. For x ∈ IR^{m} the Jacobian of F at x is denoted by JF (x), an n × m matrix with entries

(JF (x))_{i,j} = \frac{∂F_{i}}{∂x_{j}} (x).

A necessary condition for a point x ∈ IR^{m} to be a local minimum is

range(JF (x)) ∩ (−int(K)) = ∅.  (5.1)

Therefore, we will call points that satisfy the above condition (5.1) Pareto-critical points. If a point x is not Pareto-critical then there exists a direction v ∈ IR^{m} satisfying

JF (x)v ∈ −int(K),  (5.2)

i. e. a descent direction for the vector-valued objective function F. We have the following simple observation with respect to (5.2).

Lemma 5.1.1 Let K ⊆ IR^{n} be convex cone with 0 ∈ K and A ∈ IR^{n×m} be a matrix. If

Av ∈ −int(K)  (5.3)

holds, then it follows that

v ∈ (A^{⊤}(K^{∗} \ {0}))^{+}  (5.4)

holds. Moreover, if int(K) ≠ ∅ then (5.4) is sufficient for (5.3) if, additionally,

(K^{∗} \ {0}) ∩ kern(A^{⊤}) = ∅.  (5.5)

Proof: Let Av ∈ −int(K). Then int(K) ≠ ∅ and therefore

Av ∈ int(K) = {x ∈ IR^{n} | ∀ ω ∈ K^{∗} \ {0} : ⟨ ω, x ⟩ > 0},

which is equivalent to

⟨ ω, v ⟩ > 0 for all ω ∈ A^{⊤}(K^{∗} \ {0}).

From this, v ∈ (A^{⊤}(K^{∗} \ {0}))^{+} follows. Now let int(K) ≠ ∅. Then (5.5) means that

A^{⊤}(K^{∗} \ {0}) = (A^{⊤}(K^{∗} \ {0})) \ {0} = (A^{⊤}K^{∗}) \ {0}

and the result follows immediately. □

The idea of the general algorithm proposed in this chapter is straightforward: choose an x and check if (5.1) holds. If not, compute a direction v for which (5.2) holds and make a step with a suitably chosen steplength from x along v. This results in a new point, and the scheme can be repeated.

In what follows, we will consider only the case of

K = IR_{+}^{n}.

Due to this, we can abbreviate ≤:= ≤_{IR}_{+}^{n}, and we understand the inequality x < y for two vectors x, y ∈ IR^{n} componentwise, that is x_{i} < y_{i} for all i = 1, …, n. Moreover, (5.1) just means

range(JF (x)) ∩ −int(IR_{+}^{n}) = ∅.  (5.6)

5.2 Computing the Search Direction

Suppose that we have given a point x ∈ IR^{m}. Define

A := JF (x)

and the function f : IR^{m} ⟶ IR by

f(v) := max{(Av)_{i} | i = 1, …, n}.

Of course, f is convex (as the maximum of linear functions) and positive homogeneous:

f(λv) = λf(v) ∀ λ ∈ IR_{+}.

Next, consider the unconstrained minimization problem

minimize f(v) + \frac{1}{2}∥v∥_{2}^{2}
subject to v ∈ IR^{m}

(5.7)

Since the objective function is proper, closed, and strongly convex, it has always a unique solution. Note that a simple reformulation of problem (5.7) to get rid of the nondifferentiabilities would be

minimize α + \frac{1}{2}∥v∥_{2}^{2}
subject to (Av)_{i} ≤ α,  i = 1, …, n,

(5.8)

which is a convex quadratic problem with linear inequality constraints, readily solvable by a variety of methods. The next lemma shows the interconnection between Paretocriticality and the program (5.7).

Lemma 5.2.1 Let v(x) be the solution of problem (5.7) and α(x) the optimum value.

1. If x is Pareto-critical, then v(x) = 0 ∈ IR^{m} and α(x) = 0.

2. If x is not Pareto-critical, then α(x) < 0 and

f(v(x)) ≤ −\frac{1}{2}∥v(x)∥_{2}^{2} < 0,
(JF (x)v(x))_{i} ≤ f(v(x)),  i = 1, …, n.

3. The mappings x ↦ v(x) and x ↦ α(x) are continuous.

Proof: All claims are immediately clear from the definition of problem (5.7). □

Remark 5.2.2 If n = 1, we fall back to the single-objective minimization case. In this case, f(v) = ⟨ ∇F (x)), v ⟩ and

v(x) = −∇F (x).

So, we retrieve the steepest descent direction for the case n = 1. We might therefore call the descent direction computed by solving (5.7) multicriteria steepest descent direction or simply steepest descent direction.

Instead of solving problem (5.7) exactly, it is interesting for algorithmic purposes to deal with inexact solutions. So, if x is not Pareto-critical, we say that v is an approximated solution of (5.7) with tolerance σ ∈]0, 1] if

f(v) + \frac{1}{2}∥v∥_{2}^{2} ≤ σα(x)  (5.9)

where again f(v) := max_{i=1}^{n}(JF (x)v)_{i} and α(x) is the optimum value of problem (5.7). Observe that for σ = 1 only the exact solution satisfies the above inequality.

Define, for A ∈ IR^{n×m},

∥A∥_{∞,2} := max_{x≠0} \frac{∥Ax∥_{∞}}{∥x∥_{2}}

Then ∥ · ∥_{∞,2} is a norm in IR^{n×m}. It is easy to see that

∥A∥_{∞,2} := max_{i=1}^{n} ∥A_{i,·}∥_{2}
= max_{i=1}^{n} \sqrt{∑_{j=1}^{m} A_{i,j}^{2}}

Lemma 5.2.3 Suppose x is not Pareto-critical and that v is an approximated solution of (5.7) with tolerance σ ∈]0, 1]. Then

∥v∥_{2} ≤ 2∥A∥_{2,∞}.

Proof: evident □

5.2.1 Other Possibilities for the Search Direction

Of course, there is no need for the specific choice of (1/2)∥ · ∥_{2}^{2} as the function which is added to f to define problem (5.7). In fact, any proper closed strictly convex function g : IR^{m} → IR with g(0) = 0, g(v) > 0 (v ≠ 0) and g(v) = o(∥v∥) (v → 0) can be used to define the minimization problem min_{v∈IR^{m}} f(v) + g(v). The results of Lemma 5.2.1 carry over to this generalized case.

Now let us discuss another possibility for choosing descent search directions at x. Instead of solving (5.7) we could take v as a solution of

minimize f(v)
subject to ∥v∥_{∞} ≤ 1

(5.10)

This problem can also be reformulated as

minimize h
subject to (Av)_{i} ≤ h,  i = 1, …, m,
∥v∥_{∞} ≤ 1,

(5.11)

which is a linear problem. Let h(x) be the optimal function value of (5.11). We have h(x) ≤ 0, because (0, 0) ∈ IR^{m+1} is feasible for (5.11). As above, if h(x) = 0, the point x is Pareto-critical. If h(x) < 0, it is possible to compute a descent direction by solving (5.11).

Lemma 5.2.4 Let V (x) be the solution-set and h(x) be optimum value of problem (5.10). 1. If x is Pareto-critical, then 0 ∈ V (x) and h(x) = 0.

2. If x is not Pareto-critical, then h(x) < 0 and for any v ∈ V (x) we have

(JF (x)v)_{i} ≤ f(v) < 0, i = 1, …, n.

3. The mapping x ↦ h(x) is continuous.

4. Let the sequence (x^{(k)})_{k∈IN} converge to x̄, let the sequence (v^{(k)})_{k∈IN} converge to

v̄, and let v^{(k)} ∈ V (x^{(k)}) hold for all k ∈ IN. Then v̄ ∈ V (x̄) Proof: All the results are obvious. □

The choice of the constraint ∥v∥ ≤ 1 can be replaced by the more general constraint v ∈ B where B ∈ IR^{m} is a closed compact set with 0 ∈ int(B). The obvious advantage of the choice made is, of course, the simple linear structure of problem (5.11). On the other hand, every norm or gauge γ defines the set B considered above by its unit ball. We will return to this observation in Section 5.6.

No matter which auxiliary problem is used to compute a descent direction, note that there is no need to solve these auxiliary problems exactly. A descent direction is found as soon as we have a function value in the auxiliary problem of less than zero.

5.3 Computing the Step Length

Suppose that we have a direction v ∈ IR^{m} with JF (x)v < 0, i. e. (JF (x)v)_{i} < 0 for all i = 1, …, n. To compute a steplength t > 0 for the direction v, we use an Armijo-like rule. Let β ∈]0, 1[ be a prespecified constant. The condition to accept t is

F (x + tv) ≤ F (x) + βtJF (x)v.  (5.12)

We start with t = 1 and while (5.12) is not satisfied, we set

t := t/2.

Finiteness of this procedure follows from the fact that (5.12) holds strictly for t > 0 small enough.

Lemma 5.3.1 (Multicriteria Armijo Condition) If the function F is differentiable and JF (x)v < 0 holds, then there exist some ε > 0 (which may depend on x, v and β) such that

F (x + tv) < F (x) + βtJF (x)v

for all t ∈]0, ε].

Proof: Since we assume that F is differentiable, we have that

F (x + h) = F (x) + JF (x)h + R(h)

holds with

lim_{h→0} \frac{|R(h)_{i}|}{∥h∥_{2}} = 0, i = 1, …, n.

Define

α := max {(JF (x)v)_{i} | i = 1, …, n}

and observe that α < 0 and v ≠ 0. Since β < 1, there exist an ε > 0 such that

0 < t ≤ ε ⟹ \frac{|R(tv)_{i}|}{∥tv∥_{2}} < \frac{(1 − β)|α|}{∥v∥_{2}}, i = 1, …, n.

Hence, for 0 < t ≤ ε we get

|R_{i}(tv)| < t(1 − β)|α|, i = 1, …, n,

and, since |α| = −α = min_{i=1}^{n} −(JF (x)v)_{i},

R(tv) < −t(1 − β)JF (x)v.

Therefore

F (x + tv) = F (x) + t(JF (x)v) + R(tv)
< F (x) + tβ(JF (x)v) + (−t(1 − β)(JF (x)v))
= F (x) + tβ(JF (x)v)

for 0 < t ≤ ε. □

Remark 5.3.2 Let γ be an arbitrary norm in IR^{n} and γ^{∘} be the dual norm. Suppose that JF is Lipshitz-continuous and

γ^{∘} (∇F_{i}(y) − ∇F_{i}(ỹ)) ≤ Lγ (y − ỹ) , i = 1, …, n,

for all y, ỹ ∈ IR^{m}.

If JF (x)v < 0, define

ε̃:= −\frac{2(1 − β)}{L} \frac{max_{i}(JF (x)v)_{i}}{(γ(v))^{2}} > 0.

Then

F (x + tv) ≤ F (x) + βtJF (x)v

for all t ∈ [0, ε̃].

Under the conditions of the above remark, if γ is the Euclidean norm and v = v(x), then ε̃ ≥ (1 − β)/L.

5.4 The Complete Algorithm

We are now able to state the complete parameter-free descent algorithm for differentiable unconstrained multicriteria optimization problems.

Algorithm 5.4.1 (Steepest Descent for Multicriteria Problems)

1. Choose β ∈]0, 1[, σ ∈]0, 1], and x^{(1)} ∈ IR^{m}, set k := 1.

2. If x^{(k)} is Pareto-critical, stop.

3. Compute v^{(k)}, an approximative solution of (5.7) at x = x^{(k)} with tolerance σ.

4. Compute a steplength t_{k} ∈]0, 1] as the maximum of

T_{k} := {t = 1/2^{j} | j ∈ IN, F (x^{(k)} + tv^{(k)}) ≤ F (x^{(k)}) + βtJF (x^{(k)})v^{(k)}} .

5. Set x^{(k+1)} := x^{(k)} + t_{k}v^{(k)}, k := k + 1, and goto Step 2.

Observe that if Step 3 is reached at iteration k, then x^{(k)} is not Pareto-critical and

max_{i=1}^{n}(JF (x^{(k)})v^{(k)})_{i} ≤ max_{i=1}^{n}(JF (x^{(k)})v^{(k)})_{i} + \frac{1}{2}∥v^{(k)}∥_{2}^{2} ≤ σα(x^{(k)}) < 0.

So, JF (x^{(k)})v^{(k)} < 0. Hence at this iteration the set T_{k} is not empty and Step 4 is well defined. Indeed, note that the computation of the steplength in Step 4 can be done by the Armijo-like rule outlined above.

Of course, ”parameter-free” refers to the fact that the algorithm is meant to solve multicriteria optimization problems, but does not use some sort of scalarization. Any numerical method has to use some parameters, at the very least for checking the numerical accuracy of the computed approximated solution. This is what the parameters β and σ in Algorithm 5.4.1 are for.

5.5 Convergence of the Steepest Descent Method

Observe that if Algorithm 5.4.1 terminates after finitely many iterations, then it will terminate at a Pareto-critical point. From now on we suppose that an infinite sequence is generated, so α(x^{(k)}) ≠ 0 for all k ∈ IN.

We present a simple application of the standard convergence argument for the steepest descent method.

Theorem 5.5.1 Every accumulation point of the sequence (x^{(k)})_{k∈IN} produced by Algorithm 5.4.1 is a Pareto-critical point. If the function F has bounded level sets in the sense that the set {x ∈ IR^{m} | F (x) ≤ F (x^{(1)})} is bounded, then the sequence (x^{(k)})_{k∈IN} is bounded and has at last one accumulation point.

Proof. Let y be an accumulation point of the sequence (x^{(k)})_{k∈IN} and let v(y) and α(y) be the solution and the optimum value of (5.7) at y:

v(y) := arg min_{v∈IR^{m}} (f_{y}(v) + \frac{1}{2}∥v∥_{2}^{2})

and

α(y) := inf_{v∈IR^{m}} (f_{y}(v) + \frac{1}{2}∥v∥^{2})

where

f_{y}(v) := max {(JF (y)v)_{i} | i = 1, …, n} .

According to Lemma 5.2.1 it suffices to prove that α(y) = 0.

Clearly, the sequence (F (x^{(k)}))_{k∈IN} is componentwise strictly decreasing and we have

lim_{k→+∞} F (x^{(k)}) = F (y).

Therefore

lim_{k→+∞} ∥F (x^{(k)}) − F (x^{(k+1)})∥_{2} = 0

But

F (x^{(k)}) − F (x^{(k+1)}) ≥ −t_{k}βJF (x^{(k)}) v^{(k)} ≥ 0,

and therefore

lim_{k→+∞} t_{k}JF (x^{(k)}) v^{(k)} = 0  (5.13)

Observe that t_{k} ∈]0, 1] for all k. Now take a subsequence (x^{(k_{u})})_{u∈IN} converging to y. We will consider two possibilities:

lim sup_{u→+∞} t_{k_{u}} > 0

and

lim sup_{u→+∞} t_{k_{u}} = 0.

First case. Suppose that there exist a subsequence (x^{(k_{ℓ})})_{ℓ∈IN} converging to y and satisfying

lim_{ℓ→+∞} t_{k_{ℓ}} = t̄> 0

Using (5.13) we conclude that

lim_{ℓ→+∞} JF (x^{(k_{ℓ})}) v^{(k_{ℓ})} = 0,

which also implies

lim_{ℓ→+∞} α (x^{(k_{ℓ})}) = 0

Since x ↦ α(x) is continuous, we conclude that α(y) = 0, so y is Pareto-critical. Second case. Using Lemma 5.2.3 we conclude that the sequence (v^{(k_{u})})_{u∈IN} is bounded. So we can take a subsequence (x^{(k_{r})})_{r∈IN} of (x^{(k_{u})})_{u∈IN} such that the sequence (v^{(k_{r})})_{r∈IN} also converges to some v̄. Note that for any r we have

max_{i=1}^{n} (JF (x^{(k_{r})})v^{(k_{r})})_{i} ≤ σα (x^{(k_{r})}) < 0.

Therefore, passing to the limit for r → +∞ we get

\frac{1}{σ} max_{i=1}^{n} (JF (y)v̄)_{i} ≤ α(y) ≤ 0.

Take some q ∈ IN. For r large enough,

t_{k_{r}} < \frac{1}{2^{q}} ,

which means that the Armijo condition (5.12) is not satisfied for t = 1/2^{q}, i. e.

F (x^{(k_{r})} + \frac{1}{2^{q}} v^{(k_{r})}) ≰ F (x^{(k_{r})}) + β \frac{1}{2^{q}} JF (x^{(k_{r})}) v^{(k_{r})}

(for r large enough). Passing onto the limit for r → +∞ and using the fact that there are only finitely many components in F, we get

F_{j} (y + \frac{1}{2^{q}} v̄) ≥ F_{j}(y) + β \frac{1}{2^{q}} (JF (y)v̄)_{j}

for at least one j ∈ {1, …, n}. Note that this inequality holds for all q ∈ IN. From Lemma 5.3.1 it follows that

max_{i=1}^{n} (F (y)v̄)_{i} ≥ 0

which implies α(y) = 0. So, again we conclude that y is Pareto-critical. □

5.6 Towards an implementable Method

Theorem 5.5.1 guarantees that under suitable conditions the accumulation points of the sequence (x^{(k)})_{k∈IN} are Pareto-critical. But in practice, every iterative algorithm has to terminate in finite time. So, a stopping criteria other than the theoretical one used in Step 2 of Algorithm 5.4.1 should be used.

In single objective optimization, when the optimal value is neither known nor estimated a priori, the norm of the gradient of the objective function is frequently used in a stopping criterion. If this norm is small enough, the algorithm in use terminates. This, however, is not possible in multicriteria optimization since the gradients of the objective functions usually do not vanish at an efficient point.

But since the direction provided by (5.7) generalizes the steepest descent direction, the optimal value of the direction search problem (5.7) can be used as a stopping criterion in an implementation of the algorithm. In this case, we might replace Steps 2 and 3 of the algorithm by

2. Compute v^{(k)}, an approximative solution of (5.7) at x = x^{(k)} with tolerance σ. Let α(x^{(k)}) be the function value of (5.7) at v^{(k)}.

3. If α(x^{(k)}) > −τ, stop.

Here, τ > 0 is a prespecified accuracy value. As it can be seen, we do not need to compute the exact optimal function value of the direction search program, but only an upper bound for it. This is exactly what a generic minimization technique employed to solve (5.7) provides.

Nevertheless, the knowledge of the optimal value α(x^{(k)}) is useful to accept a direction v as an approximated solution. This can be seen in inequality (5.9), where the exact optimal value explicitly appears. The use of this inequality in an iterative algorithm solving the direction search program makes it necessary to compute good bounds for α(x^{(k)}).

Dual (and primal dual) algorithms (and this includes log-barrier methods) also provides lower bounds for the optimum value, which obviates the calculation of the exact optimum value of problem (5.7). Furthermore, to work with a primal-dual reformulation of problem is quite natural, because the dual of (5.8) is

maximize −\frac{1}{2} ∥∑_{i=1}^{n} θ_{i}A_{i,·}∥_{2}^{2}
subject to ∑_{i=1}^{n} θ_{i} = 1,
θ ≥ 0.

The same consideration applies to search directions based on the linear programming formulation (5.11). Here, the dual is

maximize − ∥∑_{i=1}^{n} θ_{i}A_{i,·}∥_{1}
subject to ∑_{i=1}^{n} θ_{i} = 1,
θ ≥ 0.

We can see clearly that both duals to the direction search programs proposed here are equivalent to standard goal programming problems, as considered in Chapter 4. It has already been noted in Chapter 4 that goal programming problems are usually used only to generate (some) efficient points, while they appear here as duals of ”linearized” nonlinear multicriteria problems, which might even be nonconvex. Indeed, both duals above are special cases of the general problem

minimize γ(u)
subject to u ∈ conv{A_{i,·} | i = 1, …, n},

where γ : IR^{m} ⟶ IR is a gauge as in (4.9). Using the results from Section 4.3, it is easy to see that

minimize max{−(Av)_{i} | i = 1, …, n}
subject to γ^{∘}(v) ≤ 1

is the corresponding primal problem. Clearly, both direction search programs described in Section 5.2 are special cases of this generalization.

5.7 The Constrained Case

Suppose now that a set of feasible points G ⊆ IR^{m} is given. Given z ∈ G, a vector v ∈ IR^{m} is a tangent direction of G at z if there exist a sequence (z^{(k)})_{k∈IN} with z^{(k)} ∈ G and a scalar λ ≥ 0 such that

lim_{k→+∞} z^{(k)} = z,  lim_{k→+∞} λ\frac{z^{(k)} − z}{∥z^{(k)} − z∥_{2}} = v

The set of all tangent direction of G at z is called the tangent cone of G at z and will be denoted by T (G, z).

Observe that if z is locally Pareto-optimal and v ∈ T (G, z), then we have that ⟨ ∇F_{i}(z), v ⟩ < 0 can not hold for all i = 1, …, n. So, a necessary condition for local Pareto-optimality is:

(JF (z) · T (G, z)) ∩ (−int(IR_{+}^{n})) = ∅  (5.14)

which is equivalent to saying that the system

⟨ ∇F_{i}(z), v ⟩ < 0, v ∈ T (G, z)  (i = 1, …, n)  (5.15)

has no solutions. Now suppose that G is given by

G = {x ∈ IR^{m} | g_{j}(x) ≤ 0, j = 1, …, ℓ}

with continuously differentiable functions g_{j} : IR^{m} → IR. Define I := {1, …, ℓ}. The set of active constraints at a feasible point z will be denoted by I_{0}(z), i. e.

I_{0}(z) := {i ∈ I | g_{i}(z) = 0},

and the linearized cone at z, denoted by C(z) is defined as

C(z) := {v ∈ IR^{m} | ⟨ ∇g_{i}(z), v ⟩ ≤ 0, ∀i ∈ I_{0}(z)}.

We will suppose, as a constraint qualification, that

int(C(z)) = T (G, z)

holds for all feasible z. Under this assumption, condition (5.14) (resp. (5.15)) becomes equivalent to saying that the system

\begin{aligned}⟨ ∇F_{i}(x), v ⟩ & < & 0  ∀ i = 1, …, n, \\ ⟨ ∇g_{j}(x), v ⟩ & ≤ & 0  ∀ j ∈ I_{0}(z) \end{aligned}  (5.16)

has no solution v. Let Â be the matrix obtained by adding rows (∇g_{j}(x))^{⊤} for all active constraints j ∈ I_{0}(z) to the Jacobian JF (z). Necessary for the condition above is

range(Â) ∩ (−int(IR_{+}^{n+|I_{0}(z)|})) = ∅.

We see that we can use the same steepest descent algorithm as above, by adding some rows in each step according to the active constraints, and use a step length control to enforce feasibility as well as sufficient decrease. The auxiliary program is virtually the same as (5.7), the only difference being that the matrix A is replaced by the matrix Â.

To use the above idea in a convergent algorithm, we replace (5.7) by the following program. For a given ε ∈ IR_{+} and x ∈ IR^{m} define the index set I_{ε}(x) by

I_{ε}(x) := {j ∈ I | g_{j}(x) ≥ −ε}.

We then use the direction search program

minimize α
subject to ⟨ ∇F_{i}(x), v ⟩ ≤ α,  i = 1, …, n,
⟨ ∇g_{j}(x), v ⟩ ≤ α, j ∈ I_{ε}(x),
∥v∥_{∞} ≤ 1,

(5.17)

and denote the optimal function value of this problem by α(x, ε). Following the results found in Zoutendijk [99], we have the following lemmas.

Lemma 5.7.1 Let (x^{(k)})_{k∈IN} be a sequence in G with lim_{k→∞} x^{(k)} = y and let (ε_{k})_{k∈IN} be a sequence in IR_{+} with lim_{k→∞} ε_{k} = 0. Then we have

lim sup_{k→∞} α(x^{(k)}, ε_{k}) ≤ α(y, 0).

Lemma 5.7.2 Let ∇g_{i} be Lipshitz-continuous with Lipshitz-constant L for all i ∈ I. Let x ∈ G, v ∈ IR^{m}, ε ∈ IR_{+}, and ϱ < 0 be given such that

⟨ ∇g_{j}(x), v ⟩ ≤ ρ for all j ∈ I_{ε}(x).

It then follows that

x + tv ∈ G for all t ∈ [0, min \{1, \frac{ε}{L}, −2\frac{ρ}{L}\}] .

These preliminary results enable us to prove the next theorem.

Theorem 5.7.3 Let JF and ∇g_{j} (j ∈ I) be Lipshitz-continuous and let x^{(1)} ∈ G. For a fixed ε > 0, let the sequence (x^{(k)})_{k∈IN} be generated according to the direction search program (5.17), where the step length is chosen according to the Armijo criterion (5.12). Then we either have

lim_{k→+∞} F_{i}(x^{(k)}) = −∞

for at least one index i ∈ {1, …, n}, or we have

lim_{k→+∞} α(x^{(k)}, ε) = 0.  (5.18)

Proof: We have

F (x^{(k+1)}) ≤ F (x^{(k)}) + t_{k}JF (x^{(k)})v^{(k)} + t_{k}^{2} \frac{L}{2}e,

where e = (1, …, 1)^{⊤} ∈ IR^{n} and L is the Lipshitz constant of JF and the g_{i}. Suppose that (5.18) does not hold. But we always have α(x^{(k)}, ε) ≤ 0 for all k, which means that there exists a number ϱ < 0 and a subsequence (x^{(k_{j})})_{j∈IN} with α(x^{(k_{j})}, ε) < ϱ. Then,

F (x^{(k_{j}+1)}) ≤ F (x^{(k_{j})}) + t_{k_{j}}ϱe + t_{k_{j}}^{2} \frac{L}{2}e  (5.19)

for all j ∈ IN. Without loss of generality, assume −2ϱ ≤ ε and −2ϱ ≤ L. According to Lemma 5.7.2, x^{(k)} + tv^{(k)} ∈ G holds as long as t ≤ −2ϱ/L. But Lemma 5.3.1 and the accompanying remark tell us that the Armijo criterion holds for 0 ≤ t ≤ −2ϱ(1 − β)/L, and so we have

t_{k} ∈ [−ϱ\frac{1 − β}{L} , −2ϱ\frac{1 − β}{L} ] .

Using this in (5.19), we see that

F (x^{(k_{j}+1)}) < F (x^{(k_{j})}) − ϱ\frac{1 − β}{L} (ϱ − ϱ(1 − β)) e = F (x^{(k_{j})}) − ϱ^{2}β \frac{1 − β}{L} e

holds for all j ∈ IN. As a result, lim_{k→+∞} F_{i}(x^{(k)}) = −∞. □

We can now adapt Zoutendijk’s method of feasible directions to the multicriteria case considered here.

Algorithm 5.7.4 (Method of Feasible Directions for Multicriteria Problems)

1. Choose β ∈]0, 1[, ε_{1} > 0, and x^{(1)} ∈ G, set k := 1.

2. Solve (5.17), thereby computing a direction v^{(k)} and an optimal function value α(x^{(k)}, ε_{k}) for the auxiliary problem.

3. If α(x^{(k)}, ε_{k}) ≥ −ε_{k}, set ε_{k+1} := ε_{k}/2, x^{(k+1)} := x^{(k)}, k := k + 1, and goto Step 2.

4. Compute a steplength t_{k} ∈]0, 1] such that

F (x^{(k)} + t_{k}v^{(k)}) ≤ F (x^{(k)}) + βt_{k}JF (x^{(k)})v^{(k)}

and x^{(k)} + t_{k}v^{(k)} ∈ G.

5. Set x^{(k+1)} := x^{(k)} + t_{k}v^{(k)}, k := k + 1, and goto Step 2.

Theorem 5.7.5 Let JF and ∇g_{j} (j ∈ I) be Lipshitz-continuous and let the sequence (x^{(k)})_{k∈IN} generated by Algorithm 5.7.4 be bounded. It then follows that the sequence (x^{(k)})_{k∈IN} has an accumulation point y with α(y, 0) = 0.

Proof: The proof is identical to the convergence theorem of Zoutendijk’s P1-method [99] and is therefore omitted here. □

Note that when F has bounded level sets, the sequence (x^{(k)})_{k∈IN} is necessarily bounded.

Similar to the unconstrained case, it is not necessary to solve (5.17) exactly. Instead, it is sufficient to stop as soon as a function value α has been computed with α < −ε_{k}. In case α(x^{(k)}, ε_{k}) ≥ −ε_{k} holds for the optimal function value α(x^{(k)}, ε_{k}), dual bounds can be used to stop prematurely.

5.8 A Quadratic Approximation Strategy

Let the function F be twice continuously differentiable and let the problem be unconstrained. If we use a quadratic approximation

q_{i}(v) := ⟨ ∇F_{i}(x), v ⟩ + \frac{1}{2}⟨ v, ∇^{2}F_{i}(x)v ⟩

at the point x ∈ IR^{m} for each component function F_{i} and define the function q : IR^{m} ⟶ IR^{n} by q(v) := (q_{1}(v), …, q_{n}(v))^{⊤}, we have the following necessary condition for x being Pareto-optimal:

{q(v) | v ∈ IR^{m}} ∩ (int(−IR_{+}^{n})) = ∅  (5.20)

i. e. q(v) ≮ 0 for all v ∈ IR^{m}. This condition can be checked by solving the auxiliary problem

minimize max {q_{i}(v) | i = 1, …, n}  (5.21)
subject to v ∈ IR^{m},

which can be rewritten as

minimize t  (5.22)
subject to q_{i}(v) ≤ t,  i = 1, …, n.

If the optimal function value of this auxiliary program is −∞, we can compute a common descent direction v. Otherwise, let (v, t) be a solution of the auxiliary problem. As before, if t = 0, the point x is critical. Else, we have t < 0, and v is a descent direction.

In case the Hessians ∇^{2}F_{i}(x) are all positive semidefinite, we have as an auxiliary problem a problem with convex quadratic constraints, readily solvable.

The next results show that solving the auxiliary problem is, in fact, a Newton step for a disguised scalarized problem.

Lemma 5.8.1 Let x ∈ IR^{m} and let v ∈ IR^{m} be a solution to problem (5.21). Then there exists a vector ω ∈ IR_{+}^{n} = (IR_{+}^{n})^{∗} such that v is a Newton-step for the function ⟨ ω, F (·) ⟩ from the point x.

Proof: Problem (5.21) is equivalent to Problem (5.22). The Karush-Kuhn-Tucker conditions for (5.22) can be written as

∑_{i=1}^{n} ω_{i} (∇F_{i}(x) + ∇^{2}F_{i}(x)v) = 0,  (5.23)
−1 + ∑_{i=1}^{n} ω_{i} = 0,  (5.24)
ω_{i} ≥ 0,  i = 1, …, n,  (5.25)
q_{i}(v) − t ≤ 0,  i = 1, …, n,  (5.26)

where the ω_{i} (i = 1, …, n) are the Lagrange-multiplicators. Equation (5.23) can be rewritten as

(∑_{i=1}^{n} ω_{i}∇^{2}F_{i}(x)) v = − ∑_{i=1}^{n} ω_{i}∇F_{i}(x)  (5.27)

from which the result follows by noting (5.25). □

Note that (5.24) just normalizes the vector ω into the standard base of IR_{+}^{n}. The corresponding algorithm making use of the quadratic approximation of the function F is now immediately at hand.

Algorithm 5.8.2

1. Choose x^{(1)} ∈ IR^{m}, set k := 1.

2. If (5.21) is unbounded below, stop.

3. Compute v^{(k)} ∈ IR^{m}, a solution of (5.21) at x = x^{(k)} with optimal function value t_{k}.

4. If t_{k} = 0 then stop.

5. Set x^{(k+1)} := x^{(k)} + v^{(k)}, k := k + 1, and goto Step 2.

Theorem 5.8.3 Algorithm 5.8.2 either stops after a finite number of iterations at a point for which (5.20) holds, or generates an infinite sequence (x^{(k)})_{k∈IN}. In this case, let the sequence (x^{(k)})_{k∈IN} converge to a point z such that the matrices ∇^{2}F_{i}(z) (i = 1, …, n) are positive definite. Then z is weakly Pareto-efficient.

Proof: Only the second part of the theorem has to be proven. As it has been shown in the proof of Lemma 5.8.1, for each point v^{(k)} there exist multiplicators ω^{(k)} ∈ IR_{+}^{n} such that (5.23)–(5.26) holds. Since all ω^{(k)} belong to the compact standard base of IR_{+}^{n}, there exists an accumulation point ω of this sequence. Now ω ≠ 0, and

lim_{k→+∞} v^{(k)} = lim_{k→+∞} x^{(k+1)} − x^{(k)} = 0.

Equation (5.27) together with the assumed positive definiteness of the matrices involved leads to ∇_{x}⟨ ω, F (z) ⟩ = ⟨ ω, ∇_{x}F (z) ⟩ = 0, which proves the weak Pareto-optimality of z. □

104

Bibliography

[1] Warren Adams, Marie Coffin, Peter Kiessler, K. B. Kulasekera, Matthew J. Saltzmann, and Margaret Wiecek. An optimization-based framework for affordable science and technology decisions at the office of naval research. Clemson report, Clemson University, Clemson, South Carolina, October 16 1999.

[2] G[autam] Appa and C[yril] Smith. On l_{1} and chebyshev estimation. Mathematical Programming, 5:73–87, 1973.

[3] K. J. Arrow, E. W. Barankin, and D. Blackwell. Admissible points of convex sets. In H. W. Kuhn and A. W. Tucker, editors, Contributions to the Theory of Games, pages 481–492. Princeton University Press, Princeton, New Jersey, 1953.

[4] A. Auslender. How to deal with the unbounded in optimization: Theory and algorithms. Mathematical Programming, 79:3–18, October 1997.

[5] U. Bagchi. Simultaneous minimization of mean and variation of flow time and waiting time in single machine systems. Operations Research, 37:118–125, 1989.

[6] F. L. Bauer, J. Stoer, and C. Witzgall. Absolute and monotonic norms. Numerische Mathematik, 3:257–264, 1961.

[7] Harold P. Benson and Serpil Sayin. Towards Finding Global Representations of the Efficient Set in Multiple Objective Mathematical Programming. Naval Research Logistics, 44:47–67, 1997.

[8] Å. Björck. Numerical Methods for Least Squares Problems. SIAM, Philadelphia, 1996.

[9] J[oseph] Frédéric Bonnans and Alexander Shapiro. Perturbation Analysis of Optimization Problems. Springer Series in Operations Research. Springer Verlag, New York, Berlin, Heidelberg, 2000.

[10] E[milio] Carrizosa, E. Conde, F. R. Fernández, Manuel Muñoz, and J[usto] Puerto. Pareto-optimality in linear regression. Journal of Mathematical Analysis and Applications, 190:129–141, 1995.

[11] E[milio] Carrizosa, E. Conde, M[anuel] Munoz, and J[usto] Puerto. Planar pointobjective location problems with nonconvex constraints:a geometrical construction. Journal of Global Optimization, 6:77–86, January 1995.

[12] Emilio Carrizosa and Jörg Fliege. Generalized goal programming: Polynomial methods and applications. Report BEIF/108, Vakgroep Beleidsinformatica, Department of Management Informatics, Vrije Universiteit Brussel, Pleinlaan 2, B-1050 Brussel, Belgium, January 1999. Submitted.

[13] E[milio] Carrizosa and J. B. G. Frenk. Dominating sets for convex functions with some applications. Preprint TI 96-156/9, Tinbergen Institute, The Netherlands, 1996. Published as [14].

[14] E[milio] Carrizosa and J. B. G. Frenk. Dominating sets for convex functions with some applications. Journal of Optimization Theory and Applications, 96:281–295, 1998.

[15] A. Charnes, W. W. Cooper, and R. Ferguson. Optimal estimation of executive compensation by linear programming. Management Science, 1:138–151, 1955.

[16] Andrew R. Conn, Nicholas I. M. Gould, and Philippe L. Toint. Trust-Region Methods. MPS/SIAM Series on Optimization. Society for Industrial and Applied Mathematics, Philadelphia, 2000.

[17] Indraneel Das and J. E. Dennis. Normal-Boundary Intersection: A New Method for Generating Pareto Optimal Points in Nonlinear Multicriteria Optimization Problems. SIAM Journal on Optimization, 8(3):631–657, 1998.

[18] Indraneel Das and John Dennis. A closer look at drawbacks of minimizing weighted sums of objectives for pareto set generation in multicriteria optimization problems. Structural Optimization, 14(1):63–69, 1997.

[19] J. Dauer and Y. H. Liu. Multi-criteria and goal programming. In T. Gal and H. J. Greenberg, editors, Advances in Sensitivity Analysis and Parametric Programming, Boston, 1997. Kluwer.

[20] Zvi Drezner, editor. Facility Location: A Survey of Applications and Methods. Springer Series in Operations Research. Springer Verlag, New York, 1995.

[21] R[oland] Durier. Weighting Factor Results in Vector Optimization. Journal of Optimization Theory and Applications, 58(3):411–430, September 1988.

[22] R[oland] Durier. A general framework for the one center location problem. In W. Oettli and D. Pallaschke, editors, Advances in Optimization, volume 382, pages 441–457. Springer Verlag, 1992.

[23] R[oland] Durier. The general one center location problem. Mathematics of Operations Research, 20:400–414, May 1995.

[24] Roland Durier and Christian Michelot. Geometrical properties of the Fermat-Weber problem. European Journal of Operation Research, 20:332–343, 1985.

[25] E. Eisenberg. Duality in homogeneous programming. Proceedings of the American Mathematical Society, 12:783–787, 1961.

[26] H. Eschenauer, J. Koski, and A. Osyczka. Multicriteria Design Optimization. Springer, Berlin, 1990.

[27] G. Evans. Overview of techniques for solving multiobjective mathematical programs. Management Science, 30:1268–1282, 1984.

[28] Peter C. Fishburn. Mathematics of decision theory. methods and models in the social siences. Mouton, The Hague, Paris, 1972.

[29] R. Fletcher and S. Leyffer. Nonlinear programming without penalty functions. Numerical Analysis Report NA/171, Department of Mathematics, University of Dundee, Scotland, 1997.

[30] R. Fletcher and S. Leyffer. User manual for filterSQP. Numerical Analysis Report NA/181, Department of Mathematics, University of Dundee, Scotland, 1998.

[31] Jörg Fliege. Nondifferentiability Detection and Dimensionality Reduction in Minisum Multifacility Location Problems. Journal of Optimization Theory and Applications, 94(2):365–380, August 1997.

[32] Jörg Fliege. Project report Research in Brussels: Atmospheric Dispersion Simulation. Report BEIF/116, Vakgroep Beleidsinformatica, Department of Management Informatics, Vrije Universiteit Brussel, Pleinlaan 2, B-1050 Brussel, Belgium, October 1999. Also published by the Region Bruxelles-Capital, Belgium. Available electronically via http://www.mathematik.uni-

dortmund.de/user/lsx/fliege/olaf/.

[33] Jörg Fliege. Solving Convex Location Problems with Gauges in Polynomial Time. Studies in Locational Analysis, 14:153–172, June 2000.

[34] Jörg Fliege. OLAF — a general modeling system to evaluate and optimize the location of an air polluting facility. ORSpektrum, 2001. To appear.

[35] Jörg Fliege and Andree Heseler. Constructing approximations to the efficient set of convex quadratic multiobjective problems. Unpublished Manuscript.

[36] Jörg Fliege and Stefan Nickel. An Interior Point Mehod for Multifacility Location Problems with Forbidden Regions. Ergebnisberichte Angewandte Mathematik 142, Universität Dortmund, Fachbereich Mathematik, Universität Dortmund, 44221 Dortmund, Germany, March 1997. Published as [37].

[37] Jörg Fliege and Stefan Nickel. An Interior Point Mehod for Multifacility Location Problems with Forbidden Regions. Studies in Locational Analysis, 14:23–46, 2000.

[38] Jörg Fliege and Benar Fux Svaiter. Steepest descent methods for multicriteria optimization. Mathematical Methods of Operations Research, 51:479–494, 2000.

[39] R[ichard] L. Francis and A. V. Cabot. Properties of a Multifacility Location Problem Involving Euclidean Distances. Naval Research Logistics Quarterly, 19:335–353, 1972.

[40] R. M. Freund. Dual gauge programs, with applications to quadratic programming and the minimum-norm problem. Mathematical Programming, 38:47–67, 1987.

[41] Roland W. Freund, Florian Jarre, and Siegfried Schaible. On interior-point methods for fractional programs and their convex reformulation. Technical Report 535, The Institute of Statistical Mathematics, Tokyo, Japan, 1994. Electronically

distributed via http://cm.bell-labs.com/cm/cs/doc/nam.html. Pub-

lished as [42].

[42] R[oland] W. Freund, F[lorian] Jarre, and S[iegfried] Schaible. On self-concordant barrier functions for conic hulls and fractional programming. Mathematical Programming, 74:237–246, 1996.

[43] Arthur M. Geoffrion. Proper Efficiency and the Theory of Vector Maximization. Journal of Optimization Theory and Applications, 22:618–630, 1968.

[44] C. R. Glassey. Explicit duality for convex homogeneous programs. Mathematical Programming, 10:176–191, 1976.

[45] A[lfred] Göpfert and R[einhard] Nehse. Vektoroptimierung — Theorie, Verfahren und Anwendungen, volume 74 of Mathematisch-Naturwissenschaftliche Bibliothek. BSB B. G. Teubner Verlagsgesellschaft, Leipzig, 1990.

[46] M. Gravel, I. M. Martel, R. Madeau, W. Price, and R. Tremblay. A multicriterion view of optimal ressource allocation in job-shop production. European Journal of Operation Research, 61:230–244, 1992.

[47] J. Gross and J. Talavage. A multiple-objective planning methodology for information service managers. Information Processing and Management, 15:155–167, 1979.

[48] J. Gwinner. An extension lemma and homogeneous programming. Journal of Optimization Theory and Applications, 47:321–335, 1985.

[49] Andree Heseler. Bestimmung effizienter portefeuilles durch warmstart-innere-punkt-methoden. Diplomarbeit, Lehrstuhl für Wissenschaftliches Rechnen, Fachbereich Mathematik, Universität Dortmund, 44221 Dortmund, Germany, 2000.

[50] J. B. Hiriart-Urruty and C. Lemaréchal. Convex Analysis and Minimization Algorithms, volume 1. Springer Verlag, Berlin, 1993.

[51] R. Horst and P. Pardalos, editors. Handbook of Global Optimization. Kluwer Academic Publishers, Dordrecht, 1995.

[52] R. Horst and H. Tuy. Global Optimization: Deterministic Approaches. Springer, Berlin, second edition, 1993.

[53] Angelika Hutterer and Johannes Jahn. Optimization of the location of antennas for treatment planning in hyperthermia. Preprint 265, Institut für Angewandte Mathematik, Universität Erlangen-Nürnberg, Martensstraße 3, D-91058 Erlangen, June 15 2000.

[54] H. F. Idrissi, O[dile] Lefebvre, and C[hristian] Michelot. Solving constrained multifacility minimax location problems. Working paper, Centre de Recherches de Mathématiques Statistiques et Économie Mathématique, Université de Paris 1, Panthéon-Sorbonne, Paris, France, 1991.

[55] Johannes Jahn. Mathematical Vector Optimization in Partially Ordered Linear Spaces. Verlag Peter Lang, Frankfurt, 1986.

[56] J[ohannes] Jahn. Parametric Approximation Problems Arising in Vector Optimization. Journal of Optimization Theory and Applications, 54:503–516, September 1987.

[57] C. R. Johnson and P. Nylen. Monotonicity of norms. Linear Algebra and Applications, 148:43–58, 1991.

[58] A. Jüschke, J. Jahn, and A. Kirsch. A bicriterial optimization problem of antenna design. Computational Optimization and Applications, 7:261–276, 1997.

[59] W. Kaplan and W. H. Yang. Duality theorem for a generalized Fermat-Weber problem. Mathematical Programming, 76:285–297, 1997.

[60] Helmut Laux. Entscheidungstheorie. Springer, Berlin, Heidelberg, 1998. Vierte, neubearbeitete und erweiterte Auflage.

[61] O[dile] Lefebvre, C[hristian] Michelot, and F[rank] Plastria. Geometric Interpretation of the Optimality Conditions in Multifacility Location and Applications. Journal of Optimization Theory and Applications, 65(1):85–101, April 1990.

[62] T. M. Leschine, H. Wallenius, and W. A. Verdini. Interactive multiobjective analysis and assimilative capacity-based ocean disposal decisions. European Journal of Operation Research, 56:278–289, 1992.

[63] S. T. McCormick. How to compute least infeasible flows. Mathematical Programming, 78:179–194, 1997.

[64] C[hristian] Michelot. The mathematics of continuous location. Special Isolde of Studies in Locational Analysis, J. Karkasis & B. Boffey (Eds), 5, 1993.

[65] K[aisa] Miettinen and M[arko] M. Mäkelä. Interactive bundle-based method for nondifferentiable multiobjective optimization: Nimbus. Optimization, 34:231–246, 1995.

[66] Yu[rii] Nesterov. Squared functional systems and optimization problems. Unpublished manuscript, March 4 1998.

[67] Yurii Nesterov and Arkadii Nemirowskii. Interior-Point Polynomial Algorithms in Convex Programming, volume 13 of SIAM Studies in Applied Mathematics. Society for Industrial and Applied Mathematics, Philadelphia, 1994.

[68] Yu[rii] E. Nesterov and M[ichael] J. Todd. Self scaled barrier and interior-point methods for convex programming. Mathematics of Operations Research, 22:1–42, 1997.

[69] R. W. Owens and V. P. Sreedharan. Least squares methods to minimize errors in a smooth, strictly convex norm on IR^{m}. Journal of Approximation Theory, 73:180– 198, 1993.

[70] F[rank] Plastria. Localization in single facility location. European Journal of Operational Research, 18:215–219, 1984.

[71] Frank Plastria. Private communication, May 7 1993.

[72] Frank Plastria. Continuous Location Problems. In Drezner [20], chapter 11, pages 225–262.

[73] D. Prabuddha, J. B. Gosh, and C. E. Wells. On the minimization of completion time variance with a bicriteria extension. Operations Research, 40:1148–1155, 1992.

[74] James Renegar. Linear programming, complexity theory and elementary functional analysis. Mathematical Programming, 70A:279–351, 1995.

[75] James Renegar. A Mathematical View of Interior-Point Methods for Convex Optimizations. Unpublished lecture notes., July 28, 1998.

[76] James Renegar. Linear Programming, Complexity Theory and Elementary Functional Analysis. Published as [74], June 1994.

[77] R[alph] Tyrrell Rockafellar. Convex Analysis, volume 28 of Princeton Mathematical Series. Princeton University Press, Princeton, New Jersey, 1970.

[78] C. Romero. Multiobjective and Goal Programming approaches as a distance function model. Journal of the Operational Research Society, 36:249–251, 1985.

[79] C. Romero. Handbook of Critical Issues in Goal Programming. Pergamon Press, Oxford, 1991.

[80] B. Roy. Méthodologie Multicritère d’Aide à la Décision. Economica, 1985.

[81] B. Roy and D. Bouyssou. Aide Multicritère à la Décision: Méthodes et cas. Economica, 1993.

[82] T. L. Saaty. The Analytical Hierarchy Process. McGraw-Hill, New York, 1980.

[83] Yoshikazu Sawaragi, Hirotaka Nakayama, and Tetsuzo Tanino. Theory of Multiobjective Optimization, volume 176 of Mathematics in Science and Engineering. Academic Press, 1985.

[84] M. J. Schniederjans. Goal Programming. Methodology and Applications. Kluwer, Boston, 1995.

[85] H. Späth. Mathematical Algorithms for Linear Regression. Academic Press, Boston, 1992.

[86] V. P. Sreedharan. Solutions of overdetermined linear equations which minimize error in an abstract norm. Numerische Mathematik, 13:146–151, 1969.

[87] R. E. Steuer. Multiple Criteria Optimization: Theory, Computations, and Application. Wiley, New York, 1986.

[88] R. E. Steuer. Adbase multiple objective linear programming package. University of Georgia, Athens, Georgia, 1989.

[89] M. Tamiz, D. Jones, and C. Romero. Goal programming for decision making: An overview of the current state-of-the-art. European Journal of Operational Research, 11:569–581, 1998.

[90] M. Tamiz, D. F. Jones, and E. El-Darzi. A review of Goal Programming and its applications. Annals of Operations Research, 58:39–53, 1993.

[91] Walter Vogel. Vektoroptimierung in Produkträumen, volume 35 of Systems in Economicss. Verlag Anton Hain, Meisenheim am Glan, 1977.

[92] D[ouglas] J[ohn] White. Optimality and Efficiency. John Wiley & Sons, Chichester, New York, Brisbane, Toronto, Singapore, 1982.

[93] D[ouglas] J[ohn] White. Epsilon-dominating solutions in mean-variance portfolio analysis. European Journal of Operational Research, 105:457–466, 1998.

[94] Wohlmuth, Barbara I. and Hoppe, Ronald H. W. A comparison of a posteriori error estimators for mixed finite element discretizations by Raviart-Thomas elements. Mathematical Computation, 68:1347–1378, 1999.

[95] Guoliang Xue and Yinyu Ye. An Efficient Algorithm for Minimizing a Sum of Euclidean Norms with Applications. SIAM Journal on Optimization, 7(4):1017– 1036, November 1997.

[96] Guoliang Xue and Yinyu Ye. An efficient algorithm for minimizing a sum of pnorms. SIAM Journal on Optimization, 10:551–579, 1997.

[97] E. Alper Yildirim and S. J. Wright. Warm-start strategies in interior-point methods for linear programming. Technical Report 1258, School of Operations Research and Industrial Engineering, Cornell University, 2000.

[98] P. L. Yu and M. Zeleny. The set of all nondominated solutions in linear cases and a multicriteria simplex method. Journal of Mathematical Analysis and Applications, 49:430–468, 1975.

[99] G. Zoutendijk. Methods of feasible directions. Elsevier, Amsterdam, 1960.