Constructing Approximations to the Efficient Set

of Convex Quadratic Multiobjective Problems

Jörg Fliege & Andree Heseler

January 7, 2002

Abstract

In multicriteria optimization, several objective functions have to be minimized simultaneously. For this kind of problem, no single solution can adequately represent the whole set of optimal points. We propose a new efficient method for approximating the solution set of a convex quadratic multiobjective programming problem. The method is based on a warm-start interior point algorithm for which we derive complexity results, thereby extending previous results by Yildirim & Wright. Numerical results on bicriteria problems from power plant optimization and portfolio optimization show that the method is an order of magnitude faster than standard methods applied to the problems considered.

1 Introduction

In multicriteria optimization, several objective functions have to be minimized simultaneously. Applications for these problems can be found in engineering design [14] (especially truss optimization [12] and antenna design [27]) location science [8], statistics [10], management science [15] (e. g. portfolio analysis [35], production & project planning [22, 1], and scheduling [3, 31]), environmental analysis [28, 16, 18], cancer treatment planning [25], etc. Usually, no single point will minimize all of the several objective functions given at once. As a consequence, a feasible point is defined as optimal if there does not exist a different feasible point with the same or smaller objective function values such that there is a strict decrease in at least one objective function value. However, this optimality definition is not as simple as the one from classical single-criteria optimization, because different optimal points of a multicriteria problem will be incomparable to each other in terms of the order defined above. As a consequence, no single point can represent the whole solution set, which will usually be quite large. Therefore, it is important to gain as much information as possible about the solution set of a given problem, preferably by constructing a well-defined approximation to it. This is the subject of this paper.

The necessity for a strategy like the one outlined above was also observed by Benson and Sayin [5], by Das and Dennis [13, 12], and by Hillermeier [24]. 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 [13] generates a discrete approximation to the set of efficient points, a counterexample by the authors themselves [13, Figure 3] shows that this is not the case. Berkelaar et al. [6, 7] propose a basis identification algorithm based on interior-point methodology for convex quadratic bicriteria optimization problems. However, since their algorithm computes combinatorial information of all optimal bases of the solutions, it is not clear if their algorithm is computationally or theoretically efficient.

The rest of this paper is as follows. In Section 2, the problem considered is defined, and the basic solution methodology is discussed. We will recall that the construction of the set of solutions to a multicriteria optimization problem involves solving a family of single-criteria problems, each of which can be viewed as a perturbation of every other one. Section 3 considers solving a single-criteria convex quadratic problem with an interior point method, an approach which will be put to heavy use in Section 4. In this Section 4, a warm-start strategy is derived with which several single-criteria optimization problems can be solved with less effort than previously thought. Using the results discussed in Section 2 and 4, we propose in Section 5 an efficient warm-start interior point algorithm for approximating the solution set of a convex quadratic multiobjective problem. Moreover, numerical results on several different problems types are presented and discussed.

2 Convex Quadratic Multiobjective Optimization

2.1 The Problem

Let there be given p > 1 convex quadratic objective functions of the form

f_{i}(x) = \frac{1}{2}x^{T} Q_{i}x + c_{i}^{T} x,  i = 1, …, p,  (1)

with positive semidefinite matrices Q_{i} ∈ ℝ^{n×n} and vectors c_{i} ∈ ℝ^{n} for all i. Moreover, let

G := {x ∈ ℝ^{n} | Ax = b , x ≥ 0}  (2)

be the set of feasible points. Here, A ∈ ℝ^{m×n} is a matrix with full row rank and b ∈ ℝ^{m}. We are interested in minimizing simultaneously the functions

f_{1}, …, f_{p} : G ⟶ ℝ  (3)

on the set G in a sense to be specified as follows. We have to define the notion of a solution to the given multicriteria problem. We do this by defining minimal elements of an arbitrary set in the image space.

Definition 2.1 Let M ⊆ ℝ^{p} be a nonempty set. The element y^{∗} ∈ M is called efficient, if and only if there is no other y ∈ M with

y_{i} ≤ y_{i}^{∗} ∀ i ∈ {1, 2, …, p}

and

y_{k} < y_{k}^{∗} for at least one k ∈ {1, 2, …, p}.

The set of all efficient points of the set M is called the efficient set E(M).

Define the function f : G ⟶ ℝ^{p} by f = (f_{1}, …, f_{p})^{T}. With the definition above, it becomes clear that in multicriteria optimization we are in search for the whole set E(f(G)) and, obviously, for the corresponding set of optimal decision variables f^{−1}(E(f(G))). For typical examples for this type of problem we refer to Section 5. Clearly, two efficient points f(x^{(1)}), f(x^{(2)}) ∈ E(f(G)) (x^{(1)}, x^{(2)} ∈ G) with f(x^{(1)}) ≠ f(x^{(2)}) are incomparable to each other. By their very nature of being efficient, there exist two indices i, j ∈ {1, …, p} such that f_{i}(x^{(1)}) < f_{i}(x^{(2)}) and f_{j}(x^{(2)}) < f_{j}(x^{(1)}). Therefore, just one efficient point can hardly capture the possible optimal alternatives we face when solving a multicriteria optimization problem. This shows that human decision makers need information about the whole set E(f(G)).

2.2 Scalarization

It is well-known that we can find a point close to E(f(G)) of the problem specified by (3) by solving the single-objective optimization problem

min z^{T} f(x)
subject to x ∈ G,  (4)

with z ∈ Z an arbitrary weight vector from the set

Z := {z ∈ ℝ^{p} | ∑_{i=1}^{p} z_{i} = 1, z_{i} > 0 ∀ i ∈ {1, 2, …, p}} .  (5)

This approach is often called scalarization. (For a discussion of this and other scalarization techniques see e. g. [21, 26, 24].) Indeed, defining the set of properly efficient points P by

P (f(G)) := {f(x^{∗}) | z ∈ Z, x^{∗} ∈ G, f(x^{∗}) = min_{x∈G} z^{T} f(x)} ,

it can be shown [20, 33] that

P (f(G)) ⊆ E(f(G)) ⊆ cl(P (f(G)))  (6)

holds. Here, cl(·) is the closure operator. In fact, this result holds for arbitrary functions f : G ⟶ ℝ^{p} as long as f(G) + ℝ_{+}^{p} is closed and convex. Since we can not distinguish numerically between a set and it’s closure, we can therefore replace E by P in all applications involving convex functions. Turning our attention to (3), (2), and (1), we see that we have to consider several scalar problems of the form

min \frac{1}{2}x^{T} Qx + c^{T} x
subject to Ax = b,  (7)
x ≥ 0,

where Q = ∑_{i=1}^{p} z_{i}Q_{i}, c = ∑_{i=1}^{p} z_{i}c_{i}, and z = (z_{1}, …, z_{p})^{T} ∈ Z is a given parameter or weight vector. Of course, in applications we have to content ourselves with a finite subset Z̃ ⊂ Z from which we have to choose our weight vectors. The corresponding objective functions are once again convex quadratic, while the set of feasible points is polyhedral.

Note that there exists other scalarization techniques, notably the so-called maxscalarization, in which P is replaced by

{f(x^{∗}) | z ∈ Z, x^{∗} ∈ G, f(x^{∗}) = inf_{x∈G} max_{i=1}^{p} z_{i}f_{i}(x)} .

For this and other scalarizations, a result similar to (6) holds, too. Since we do not want to discuss the peculiarities of different scalarization techniques here, we stick for the sake of a simplified exposition to the weighted-sum scalarization outlined above. Note, however, that all results presented here hold for different scalarizations, too. For more details about different scalarization techniques, the reader is referred to [24, 9, 17].

As a consequence of the discussion above, we are able to approximate the set E(f(G)) as well as f^{−1}(E(f(G))) by solving optimization problems of the form (7). These ersatz problems are defined by choosing different weights z ∈ Z, see (4).

The basic idea is now as follows. Our aim is to compute a discrete approximation of the set of efficient points. We have to solve a standard scalar optimization problem for each efficient point we want to compute. The different optimization problems we have to consider can be viewed as perturbations of each other, with vectors of weights z ∈ Z serving as parameters defining the perturbations. Unfortunately, it is in general not true that a small perturbation in the weight vector leads to a small perturbation of a corresponding efficient point. As a consequence, we have to be careful when discretizing the set Z. A uniform a priori discretization of the set of weights will usually result in a waste of a large amount of computational effort, since it is not a priori clear in which part of the weight space a finer discretization (i. e. a finer grid) is necessary. On the other hand, it has to be expected that information gained during the computation of some efficient points can be put to use while computing other efficient points, e. g. by choosing the discretization of the set of weights in an adaptive way. Therefore, we propose to use an adaptive discretization technique for the set of weights Z. Basically, we want to use more parameters in those regions of the parameter space where weight vectors which are close together result in efficient points in the image space ℝ^{p} which are far apart from each other. Furthermore, to save work when computing the new efficient points (i. e. when solving the new optimization problems), we propose to use a warm-start strategy. With such a strategy, points from the iteration history of scalar problems already solved are used as starting points for the optimization problems currently under consideration.

3 An Interior-Point Algorithm for Convex Quadratic

Optimization Problems

This section will review in short an interior-point method for solving convex quadratic optimization problems of the form (7). The method we will use is an adapted version of Wright’s primal-dual long step algorithm for linear programs, explained in detail in [36, Chapter 5]. We will consider here only the major differences between the linear case described in [36] and the case of a quadratic objective function needed in what follows.

3.1 The KKT System

The dual problem to (7) is

max −\frac{1}{2}x^{T} Qx + b^{T} λ
subject to −Qx + A^{T} λ + s = c,  (8)
s ≥ 0,

where the dual variables are s ∈ ℝ^{n} and λ ∈ ℝ^{m}. For a primal-dual feasible point (x, λ, s) the inner product s^{T} x is the duality gap between x and (λ, s). Moreover, we have strong duality, i. e. the duality gap is zero in a primal-dual solution, and for arbitrary primal-dual feasible points (x, λ, s) the chain of equalities/inequalities

s^{T} x = c^{T} x − b^{T} λ + x^{T} Qx ≥ 0

holds. The corresponding KKT-system can be written as

F (x, λ, s) := ( \begin{matrix}−Qx + A^{T} λ + s − c \\ Ax − b \\ SXe \end{matrix} ) = ( \begin{matrix}0 \\ 0 \\ 0 \end{matrix} ) ,  (9)
(x, s) ≥ 0,  (10)

with e := (1, 1, …, 1)^{T} ∈ ℝ^{n}. As usual, the matrices X and S are defined by X := diag(x_{1}, x_{2}, …, x_{n}) and S := diag(s_{1}, s_{2}, …, s_{n}). The set of primal-dual feasible points is given by

Ω := {(x, λ, s) ∈ ℝ^{n} × ℝ^{m} × ℝ^{n} | Ax = b, −Qx + A^{T} λ + s = c, (x, s) ≥ 0},

and we denote the set of strictly feasible points by

Ω^{0} := {(x, λ, s) ∈ Ω | (x, s) > 0}.

Note that the Jacobian of F is regular at points lying in Ω^{0}, and that therefore a Newton step for F is well defined for these points. Perturbing the last block of the KKT-system in the usual way with some τ > 0 leads to

( \begin{matrix}−Qx + A^{T} λ + s − c \\ Ax − b \\ SXe \end{matrix} ) = ( \begin{matrix}0 \\ 0 \\ τe \end{matrix} ) ,  (x, s) > 0.  (11)

For fixed τ > 0, the (unique) point in Ω^{0} solving (11) is denoted by (x_{τ}, λ_{τ}, s_{τ}). The set {(x_{τ}, λ_{τ}, s_{τ}) ∈ Ω^{0} | τ > 0} is the central path.

Denote by I ∈ ℝ^{n×n} the identity matrix. In each iteration of the long-step interior-point algorithm presented below, we do a Newton step by solving a linear system of the form

( \begin{matrix}−Q & A^{T} & I \\ A & 0 & 0 \\ S & 0 & X \end{matrix} ) ( \begin{matrix}Δx \\ Δλ \\ Δs \end{matrix} ) = ( \begin{matrix}0 \\ 0 \\ −XSe + σμe \end{matrix} ) ,  (12)

where

μ = μ(x, λ, s) := x^{T} s/n  (13)

describes the duality gap of the point (x, λ, s) ∈ Ω^{0}. The centering parameter σ controls the bias of the step made towards the central path.

3.2 The Long-Step Algorithm

Before presenting the long-step algorithm for convex-quadratic problems, we define the wide neighborhood of the central path in the usual way, i. e.

N_{−∞}(γ) := {(x, λ, s) ∈ Ω^{0} | x_{i}s_{i} ≥ γμ, i = 1, 2, …, n}.

Given (x, λ, s) ∈ Ω^{0} and a step vector (Δx, Δλ, Δs) satisfying the system (12), we define

(x(α), λ(α), s(α)) := (x, λ, s) + α(Δx, Δλ, Δs)  (14)

as well as

μ(α) := \frac{x(α)^{T} s(α)}{n}.  (15)

The following algorithm is based on Wright’s long-step algorithm for linear programs, see [36, p. 96f].

Algorithm I: Long-step path-following method

Initialization: Define k := 0 and choose ϵ > 0, σ_{min}, σ_{max} with 0 < σ_{min} < σ_{max} < 1, γ ∈] 0, 1[, as well as (x^{0}, λ^{0}, s^{0}) ∈ N_{−∞}(γ).

Loop: While

μ_{k} := \frac{(x^{k})^{T} s^{k}}{n} > ϵ  (16)

do the following:

Compute Step: choose σ_{k} ∈ [σ_{min}, σ_{max}] and solve

( \begin{matrix}−Q & A^{T} & I \\ A & 0 & 0 \\ S^{k} & 0 & X^{k} \end{matrix} ) ( \begin{matrix}Δx^{k} \\ Δλ^{k} \\ Δs^{k} \end{matrix} ) = ( \begin{matrix}0 \\ 0 \\ −X^{k}S^{k}e + σ_{k}μ_{k}e \end{matrix} ) . (17)

Compute step length: set

α_{k}^{∗} := \{ \begin{matrix}\frac{n(1 − σ_{k})}{2(Δx^{k})^{T} Δs^{k}} μ_{k} & if (Δx^{k})^{T} Δs^{k} ≥ n(1 − σ_{k})μ_{k}/2, \\ 1 & else. \end{matrix}  (18)
choose α_{k} as the largest α ∈ [0, α_{k}^{∗}] such that
(x^{k}(α), λ^{k}(α), s^{k}(α)) ∈ N_{−∞}(γ).

Update: Set (x^{k+1}, λ^{k+1}, s^{k+1}) := (x^{k}(α_{k}), λ^{k}(α_{k}), s^{k}(α_{k})) and k := k + 1.

End (Loop)

Comparing the algorithm described above with the one for the linear case, we see that there are two differences. In the linear case we have Q = 0 and due to this, Δx^{T} Δs = 0 as well as μ(α) = (1 − α(1 − σ))μ holds (see [36, Lemma 5.1]). But both equalities break down in the more general case considered here, and we only have

Δx^{T} Δs ≥ 0  (19)

as well as

μ(α) = (1 − α(1 − σ))μ + α^{2} \frac{Δx^{T} Δs}{n}.  (20)

This is the reason why the computation of the step length α has to be modified from the one presented in [36] to (18). A constant step length of α = 1 works only in the linear case, while in the quadratic case the important inequality μ(1) < μ(0) which is used when proving convergence of the method might break down. Fortunately, if μ(α) is nonlinear, we can choose α^{∗} > 0 as the minimizer of μ(α) (cmp. (20) with (18)). By this, we get a large reduction of the duality gap. It is easy to show that μ(α) remains strictly positive for all α ∈] 0, 1] and that μ(α) < μ(0) holds for all α ∈] 0, α^{∗}].

Similar to the linear case, there exists a constant δ > 0 depending only on γ, σ_{min}, and σ_{max} (but not on n) such that

μ_{k+1} ≤ (1 − \frac{δ}{n}) μ_{k}  (21)

holds for all k ≥ 0. This can be shown in the same way as Theorem 5.11 of [36].

Clearly, the iterates produced by this algorithm are strictly feasible. The convergence properties as well as the computational complexity are discussed in Section 4.3, p. 16.

4 A Warm-Start Strategy

Suppose now that two optimization problems of the form (7) are given. If we have solved one of them by an interior-point algorithm, we might view the second problem as a perturbation of the first one and try to use one of the iterates calculated while solving the first problem as a starting point for the second one. It is exactly this situation that we are in if we want to construct an approximation to the set of efficient points of a given convex quadratic multicriteria optimization problem. Of course, the smaller the perturbation of the first problem is, the higher are our hopes that an iteration point generated rather late in the algorithm (i. e. close to a solution of the first problem) is close to a suitable starting point for the second problem, thereby saving a large amount of computational effort. Such a suitable starting point for the second problem is called a warm-start point. We now turn our attention to generating a strictly feasible point for the second optimization problem out of a strictly feasible point for the first optimization problem.

For the warm-start strategy outlined in this section we introduce notation and ideas adapted from the ones given by Nunez and Freund [34] and Yildirim and Wright [39].

4.1 A Warm-Start Step

A data instance

d = (A, b, Q, c)

represents the optimization programs (7) and (8). We define the norm of d by

∥d∥_{2} := max{∥A∥_{2}, ∥b∥_{2}, ∥Q∥_{2}, ∥c∥_{2}}.

Perturbations in the objective function are represented by a data instance like Δd = (0, 0, ΔQ, Δc) with a symmetric matrix ΔQ. Of course, if the problem changes, the feasible set and the wide neighborhood as defined in the last section changes, too. Therefore, we define Q̃ := Q + ΔQ,

Ω̃(Δd) := {(x, λ, s) | Ax = b, A^{T} λ + s = Q̃x + c + Δc, (x, s) ≥ 0},
Ω̃^{0}(Δd) := {(x, λ, s) ∈ Ω̃(Δd) | (x, s) > 0},

and

Ñ_{−∞}(γ, Δd) := {(x, λ, s) ∈ Ω̃^{0}(Δd) | ∀i = 1, 2, …, n : x_{i}s_{i} ≥ γ(x^{T} s/n) = γμ}.

Suppose now that we are given a point w := (x, λ, s) ∈ N_{−∞}(γ). Then we are in search for a corrector step Δw := (Δx, Δλ, Δs) such that

w + Δw =: w̃ = (x̃, λ̃, s̃) ∈ Ñ_{−∞}(γ, Δd)

holds.

A corrector step will be obtained by solving the following system of equations:

( \begin{matrix}−Q̃ & A^{T} & I \\ A & 0 & 0 \\ S & 0 & X \end{matrix} ) ( \begin{matrix}Δx \\ Δλ \\ Δs \end{matrix} ) = ( \begin{matrix}Δc + ΔQx \\ 0 \\ 0 \end{matrix} ) .  (22)

This is a generalization of the system of Yildirim and Todd [37], who considered linear programs only. The first and second block of (22) ensures that w̃ fulfills the primal-dual equality constraints of the perturbed problem represented by d + Δd. The third block ensures that the duality gap of w̃ is at least as small as the duality gap of w, see Theorem 4.1 below.

The corresponding augmented system is

Δx = −D^{2}Δs,  (23)
Δs = B(ΔQx + Δc − A^{T} Δλ),  (24)
Δλ = (AD^{2}BA^{T})^{−1}(AD^{2}B)(Δc + ΔQx),  (25)

where D^{2} := S^{−1}X and

B := (I + Q̃D^{2})^{−1}.  (26)

4.2 Feasibility of the Warm-Start Step

The new point w̃ = w + Δw is not necessarily strictly feasible with respect to the new problem represented by d+Δd. In fact, it might happen that some components of x̃ = x + Δx and s̃ = s + Δs are nonpositive. Changing the step size in the update w̃ = w + Δw to some value not equal to 1 does not help, since then w̃ does not necessarily fulfill the dual equality constraints. (Cmp. with the first block of (22).) However, by taking a closer look at the third block of (22), we can derive a condition equivalent to strict positivity of (x̃, s̃).

Clearly, x̃_{i} > 0 resp. s̃_{i} > 0 is equivalent to Δx_{i}/x_{i} > −1 resp. Δs_{i}/s_{i} > −1. This in turn is equivalent to

∥X^{−1}Δs∥_{∞} = ∥S^{−1}Δs∥_{∞} < 1.  (27)

(Cmp., e. g., the proof of Proposition 5.1 in [39].) But from a practical point of view, checking (27) is quite expensive. This is due to the fact that the definition of Δs and Δx involves the matrix B defined in (26). Moreover, the warm-start step has to be computed before we can check if this step leads to a feasible point. Hence, some simple sufficient conditions, more easier to handle, are needed. The next theorem shows such a weaker condition.

Theorem 4.1 Let there be given a data instance d = (A, b, Q, c) and a strictly feasible point w = (x, λ, s) for d. Furthermore, let Δd = (0, 0, ΔQ, Δc) be a perturbation and let Δw = (Δx, Δλ, Δs) be a solution to (22). Define

w̃ = (x̃, λ̃, s̃) = w + Δw.  (28)

and

T := I − IA^{T} (AD^{2}BA^{T})^{−1}AD^{2}B.

If

∥( \begin{matrix}Δc \\ ΔQx \end{matrix} )∥_{∞} < ∥S^{−1}[BT, BT]∥_{∞}^{−1}  (29)

holds, then w̃ is strictly feasible for the perturbed problem with data instance d+Δd. Moreover,

x̃^{T} s̃ ≤ x^{T} s  (30)

holds.

Proof: The proof follows the proof of Proposition 5.1 from Yildirim and Wright [39]. We just have to replace the linear augmented system by the system (23)–(25). □

In applications it will usually make sense to replace (27) by

∥S^{−1}Δs∥_{∞} ≤ 1 − θ,

where θ ∈] 0, 1[. Now we are able to describe in which wide neighborhood the warm-start point w + Δw will be.

Theorem 4.2 Let there be given w, Δw, d, Δd, B, and T as in Theorem 4.1 and let θ ∈] 0, 1[. If

∥( \begin{matrix}Δc \\ ΔQx \end{matrix} )∥_{∞} ≤ \frac{1 − θ}{∥S^{−1}[BT, BT]∥_{∞}}  (31)

holds, then w̃ = (x̃, λ̃, s̃) = w + Δw is a strictly feasible point of the optimization problem associated with the perturbed data instance d + Δd and we have x̃^{T} s̃ ≤ x^{T} s. If γ_{0} ∈] 0, 1[ and w ∈ N_{−∞}(γ_{0}), then

w̃ ∈ Ñ_{−∞}(θγ_{0}, Δd).

Proof: The proof follows the lines of the proof of Proposition 5.3 of [39], invoking Theorem 4.1 instead of [39, Proposition 5.1] at the appropriate instances. □

Next, we derive a sufficient condition for (27), leading to a complexity analysis of a corresponding warm-start algorithm.

The set

L := {(A, b, Q, c) | ∃ (x, λ, s) : (x, s) > 0, Ax = b, −Qx + A^{T} λ + s = c} (32)

contains those data instances for which strictly primal-dual feasible points exist. Hence, the complement L^{C} of L contains those instances for which there does not exist a strictly primal-dual feasible point. The intersection of the closure of the two sets,

B := cl(L) ∩ cl(L^{C}).  (33)

is the boundary between those data instances for which strictly feasible points exists and those for which no such point exists. Note that (0, 0, 0, 0) ∈ B. Following Nunez & Freund [34], we define the distance to ill-posedness

ρ(d) := inf{∥Δd∥_{2} | d + Δd ∈ B}

and the condition number of a data instance d ∈ L,

C(d) := \frac{∥d∥_{2}}{ρ(d)}

(we will use C(d) = ∞ if ρ(d) = 0). By noting that Δd = −d leads to d+Δd = 0 ∈ B we see that ρ(d) ≤ ∥d∥_{2} holds. Due to this, we always have C(d) ≥ 1.

Lemma 4.3 Let there be given A ∈ ℝ^{m×n}, Q ∈ ℝ^{n×n}, and c ∈ ℝ^{n}. The following systems can not be solved simultaneously:

1.

A^{T} λ < c + Qx

2.

Ax = 0,
x ≥ 0,
c^{T} x + x^{T} Qx ≤ 0,
x ≠ 0.

Proof: Assume that both systems can be solved simultaneously. Then we have due to System 1

c > A^{T} λ − Qx.

But with System 2 we get

c^{T} x + x^{T} Qx > λ^{T} Ax − x^{T} Qx + x^{T} Qx = 0,

which is a contradiction. □

Lemma 4.4 Let d = (A, b, Q, c) ∈ L be a data instance with ρ(d) > 0 and A ≠ 0. Let w = (x, λ, s) be a primal-dual feasible point for the problem described by d. Then we have that

∥x∥_{2} ≤ C(d) + \frac{|c^{T} x + x^{T} Qx|}{ρ(d)}  (34)

holds.

Proof: (Following an idea of Nunez and Freund [34, p. 7f].) Given the data instance d, let w = (x, λ, s) be a point satisfying the assumptions. Assume that x ≠ 0. Define ΔA := −bx^{T} /∥x∥_{2}^{2}, Δc := (−|c^{T} x + x^{T} Qx|/∥x∥_{2}^{2})x and the problem instance d + Δd by Δd := (ΔA, 0, 0, Δc). With this, we have (A + ΔA)x = 0 and

(c + Δc)^{T} x + x^{T} Qx = c^{T} x + x^{T} Qx − |c^{T} x + x^{T} Qx| ≤ 0.

But since x ≥ 0 as well as x ≠ 0, Lemma 4.3 tells us that there does not exist strictly feasible dual points of d + Δd, i. e. there does not exist a pair (x, λ) with (A + ΔA)^{T} λ < c + Δc + Qx. This leads to

ρ(d) ≤ ∥Δd∥_{2} = max{∥ΔA∥_{2}, ∥Δc∥_{2}}
≤ \frac{max\{∥b∥_{2}, |c^{T} x + x^{T} Qx|\}}{∥x∥_{2}}
≤ \frac{∥d∥_{2} + |c^{T} x + x^{T} Qx|}{∥x∥_{2}}.

From this, the conclusion follows. □

Now we can formulate the following sufficient condition for (27).

Theorem 4.5 Let there be given the data instances d and Δd as in Theorem 4.1 and a parameter γ_{0} ∈] 0, 1[ as well as a constant q > 0. For feasible points w = (x, λ, s) of the instance d define μ = μ(w) = x^{T} s/n. Then, for the right hand side of (29) we have that

∥S^{−1}[BT, BT]∥_{∞}^{−1} ≥ \frac{γ_{0}μ}{2n^{1/2} (C(d) + \frac{q}{ρ(d)})}  (35)

holds for all w = (x, λ, s) ∈ N_{−∞}(γ_{0}) with |c^{T} x + x^{T} Qx| ≤ q. (Note that S, B, and T depend on w.)

Proof: We have the estimate

∥S^{−1}[BT, BT]∥_{∞} ≤ 2∥S^{−1}BT ∥_{∞}.

Since the matrix T = I − A^{T} (AD^{2}BA^{T})^{−1}AD^{2}B is a projection matrix onto the nullspace of AD^{2}B, the upper bound ∥T ∥_{2} ≤ 1 follows immediately. In view of ∥T ∥_{∞} ≤ n^{1/2}∥T ∥_{2}, it follows that

2∥S^{−1}BT ∥_{∞} ≤ 2n^{1/2}∥S^{−1}∥_{∞}∥B∥_{2}  (36)

holds. Next, we have ∥B∥_{2} = 1/|ρ_{min}(B^{−1})|, where ρ_{min}(B^{−1}) denotes the smallest singular value of B^{−1}. Since D^{2}Q̃+Q̃D^{2} as well as D^{2}Q̃Q̃D^{2} is positive semidefinite, we obtain

ρ_{min}(B^{−1})^{2} = min_{x≠0} \frac{x^{T} (I + Q̃D^{2})^{T} (I + Q̃D^{2})x^{T}}{x^{T} x}
= min_{x≠0} \frac{x^{T} (I + D^{2}Q̃ + Q̃D^{2} + D^{2}Q̃Q̃D^{2})x}{x^{T} x}
= 1 + min_{x≠0} \frac{x^{T} (D^{2}Q̃ + Q̃D^{2} + D^{2}Q̃Q̃D^{2})x}{x^{T} x}
≥ 1.

Therefore, we arrive at

∥B∥_{2} ≤ 1,

and (36), x_{i}s_{i} ≥ γ_{0}μ as well as Lemma 4.4 lead to

2n^{1/2}∥S^{−1}∥_{∞}∥B∥_{2} ≤ max_{i=1}^{n} \frac{2n^{1/2}}{s_{i}}
= max_{i=1}^{n} \frac{2n^{1/2}x_{i}}{x_{i}s_{i}}
≤ 2n^{1/2} (C(d) + \frac{q}{ρ(d)}) \frac{1}{γ_{0}μ}.
□

Now we can estimate how small a perturbation and the duality gap of a primaldual feasible point has to be in order to reach a wide neighborhood of the perturbed problem with the corrector step. In the following theorem, we will make use of ∥Δd∥_{∞} = max{∥Δc∥_{∞}, ∥ΔQ∥_{∞}}.

Theorem 4.6 Let there be given γ and γ_{0} with 0 < γ < γ_{0} < 1, θ := γ/γ_{0}, and d, Δd, w as well as w̃ as in Theorem 4.2. Suppose that there are given bounds ξ and q with ξ > 1, ξ C(d) ≥ ∥x∥_{2} and |c^{T} x + x^{T} Qx| ≤ q. If w ∈ N_{−∞}(γ_{0}) and

μ ≥ \frac{2}{γ_{0}(1 − θ)}ξn^{1/2}C(d) (C(d) + \frac{q}{ρ(d)}) ∥Δd∥_{∞}  (37)

holds, then

w̃ ∈ Ñ_{−∞}(γ, Δd).

Proof: The inequality (37) can be written as

∥Δd∥_{∞}ξ C(d) ≤ \frac{γ_{0}(1 − θ)μ}{2n^{1/2} (C(d) + \frac{q}{ρ(d)})}

Moreover, we have

∥( \begin{matrix}Δc \\ ΔQx \end{matrix} )∥_{∞} ≤ max\{∥Δc∥_{∞}, ∥ΔQ∥_{∞}∥x∥_{2}\}
≤ max{∥Δd∥_{∞}, ∥Δd∥_{∞}ξ C(d)}
= ∥Δd∥_{∞}ξ C(d),

and after invoking Theorem 4.5, we see that

∥( \begin{matrix}Δc \\ ΔQx \end{matrix} )∥_{∞} ≤ \frac{γ_{0}μ(1 − θ)}{2n^{1/2} (C(d) + \frac{q}{ρ(d)})} ≤ \frac{1 − θ}{∥S^{−1}[BT, BT]∥_{∞}}

holds. Theorem 4.2 tells us that w̃ = w + Δw ∈ Ñ_{−∞}(θγ_{0}, Δd) = Ñ_{−∞}(γ, Δd), where Δw = (Δx, Δλ, Δs) is defined by (22). □

Note that we need an estimate for the slightly enigmatic bound ξ before we can use Theorem 4.6. This, however, turns out to be not that difficult. Indeed, under the assumptions of Lemma 4.4 we can define q := |c^{T} x + x^{T} Qx| to get

∥x∥_{2} ≤ C(d) + \frac{q}{ρ(d)},

and therefore we can use

ξ := 1 + \frac{q}{∥d∥_{2}}  (38)

as a suitable constant ξ for which

ξ C(d) ≥ ∥x∥_{2}, ξ > 1

holds.

4.3 Complexity

Suppose we have given a problem as a data instance d as well as a primal-dual feasible starting point (x^{(0)}, λ^{(0)}, s^{(0)}) with the duality measure μ_{0} > 0. Moreover, let ϵ > 0 be given. We want to stop iterating as soon as for an iterate (x, λ, s) with duality measure μ = x^{T} s/n the inequality

μ ≤ ϵ∥d∥_{2}.  (39)

holds. Due to (21), Algorithm I stops after

K := O (n log \frac{μ_{0}}{ϵ∥d∥_{2}} )

iterations. Afterwards, we have (39).

Now, we want to a use warm-start point for solving a perturbed problem represented by the data instance d + Δd. We have already generated iterates w^{(k)} = (x^{(k)}, λ^{(k)}, s^{(k)}) (k = 0, …, K), primal-dual feasible for d. Assume we have given a perturbation Δd, an iterate w^{(j)}, j ∈ {0, …, K}, as well as a corresponding duality gap μ_{j}, satisfying the assumptions of Theorem 4.6. Then the interior-point algorithm applied to the problem with data d + Δd, using the starting point w^{(j)} + Δw generated by the warm-start method stops after

O (n log \frac{μ_{j}}{ϵ∥d + Δd∥_{2}} )

iterations. Finally, assume that ∥Δd∥_{2} ≤ ∥d∥_{2}/2. This yields

\frac{1}{∥d + Δd∥_{2}} ≤ \frac{1}{∥d∥_{2} − ∥Δd∥_{2}} ≤ \frac{2}{∥d∥_{2}} ,

and we need

O (n log \frac{μ_{j}}{ϵ∥d∥_{2}} )

iterations. Note, however, that the size of μ_{j} is directly linked with ∥Δd∥_{∞} as well as C(d) via (37). A suitable small perturbation ensures that the allowed μ_{j} will be small, too. Moreover, it can be clearly seen that the warm-start method will be more efficient the smaller the perturbation is.

5 Numerical Results

Based on Theorem 4.6, we can generate an approximation to the set of efficient points to a given convex quadratic multicriteria optimization problem in the following way. (For the sake of a simplified exposition and a less cumbersome notation we consider here the bicriteria case.) The accuracy ε > 0 of the interior-point method is chosen beforehand. Next, we generate a strictly primal-dual feasible εapproximation w to a (weakly) efficient point by solving the scalarized optimization problem d corresponding to the parameter z = (1, 0)^{T} (cmp. Subsection 2.2). Then, we choose a parameter vector z +Δz (and thereby a scalarized problem d+Δd) with Δz sufficiently small such that w can be used to generate a warm-start point w̃ for d + Δd by way of Theorem 4.6. An estimate on the size of Δz can be calculated by using Theorem 4.5. The problem d + Δd is then solved from the starting point w̃, and we proceed inductively until we are sufficiently close to the parameter (0, 1)^{T}.

This strategy can be easily generalized to convex quadratic multicriteria problems with an arbitrary number of criteria.

In what follows, we compare the warm-start strategy for approximating the set of efficient points described in this paper with a so-called ”naive” strategy. This simpler strategy consist of using the scalarization approach from Section 2.2 to solve exactly the same scalar optimization problems with Algorithm I that are considered in the warm-start strategy outlined above. However, the only information reused in this naive strategy consists of the starting points.

The code was implemented in Matlab Version 5. The numerical tests described in what follows were done on a 1GHz Pentium PC with Linux operating system. For the stopping criterion we used ε = \sqrt{eps} throughout all our tests (eps the machine precision).

5.1 Power Plant Optimization

Power plant optimization considers the following problem: given a number of electrical power plants and a prognosis of the electricity demand for a time interval [t_{0}, t_{0} + T] (T > 0), compute a control function for each power plant such that

• the total energy production cost is minimized,

• the wear & tear of the power plants is minimized,

• the prognosis is met, and

• all technical constraints met.

This is a bicriteria optimization problem. Basically, such an optimization problem has to be solved all the time, since the data given (especially the prognosis data) is constantly updated. The time horizon T varies from 60 sec. to 20 minutes. Since the actual control unit of a power plant is comparatively simple, we use piecewise constant control functions. With this discretization we arrive at two convex quadratic objective functions which have to be minimized subject to linear equality and inequality constraints. More details on the model and the corresponding optimization problems can be found in [2, 19].

Using real-world data from Northrhine-Westfalia, a region in Germany containing 14 power plants and about 18 million private electricity consumers, we generated 1440 optimization problems, one per minute t_{0} of a given 24 hour period.

As it can be seen from Figure 1 (p. 19) our warm-start strategy is consistently about 5 times faster than the standard method reusing only the starting point. The number of efficient points generated per problem varies roughly between 200 and 300 points, see Figure 2 (p. 20).

Figure 1: Computation times of the warm-start strategy described here and a ”naive” strategy reusing just the starting points for the power plant problems. The x-axis indexes the different optimization problems according to t_{0}, the start of the control period on the given day. As it can be seen, the warm-start strategy is about 5 times faster than reusing the starting points only.

5.2 Portfolio Optimization

For our second set of numerical tests we consider portfolio problems derived from the standard Markowitz model, see [30]. For bicriteria portfolio problems of this type there exist already efficient methods, notably the critical line algorithm [30, 29, 32]. Note, however, that the applicability of this algorithm relies crucially on the fact Figure 2: Number of efficient points generated with the warm-start strategy for the power plant problems considered. The x-axis indexes the different optimization problems according to t_{0}, the start of the control period on the given day.

that one of the objective functions is linear. This is not the case for the more general warm-start method proposed here.

Table 1 (p. 21) describes the results of our numerical tests on the Beasley portfolio problem test set [4]. The number of assets corresponds to the number of unknowns in the primal problem. As it can be seen by the number of floating point operations needed to generate the approximation to the set of efficient points, the strategy presented here is several orders of magnitude faster than the naive strategy outlined above.

Table 2 (p. 21) shows the results of our numerical tests on historical data, see [23]. As it can be seen, the speed-ups differ widely between a factor of 4 and three orders of magnitude. Our conjecture is that the speed-up depends strongly on the influence of Q on the conditioning of the system matrix of (12). Further theoretical as well as numerical investigations with respect to this and other phenomena are under away.

No. Flops/pt. Flops/pt. Sec/pt. Speed-up
of naive new new factor
Assets strategy strategy strategy (in flops)
31  7.5 × 10^{7} 3.1 × 10^{6} 4.7 × 10^{−3}  240
85  1.8 × 10^{8} 6.9 × 10^{6} 2.2 × 10^{−2}  264
89  2.1 × 10^{9} 5.6 × 10^{6} 2.2 × 10^{−2}  381
98  3.4 × 10^{9} 7.6 × 10^{6} 2.2 × 10^{−2}  445
225 6.5 × 10^{10} 9.2 × 10^{7} 1.9 × 10^{−1}  705

Table 1: Summary of numerical results for the Beasley portfolio test set. Flops/pt.: floating point operations per computed efficient point. Sec/pt.: seconds of runtime per computed efficient point.

No. Flops/pt. Flops/pt. Sec/pt. Speed-up
of naive new new factor
Assets strategy strategy strategy (in flops)
100 6.0 × 10^{8} 1.3 × 10^{7} 1.5 × 10^{−1}  42
100 3.1 × 10^{8} 1.4 × 10^{7} 1.7 × 10^{−1}  36
100 4.6 × 10^{8} 3.8 × 10^{7} 4.1 × 10^{−1}  10
100 4.0 × 10^{8} 3.8 × 10^{7} 4.5 × 10^{−1}  7
150 1.9 × 10^{9} 5.5 × 10^{6} 9.3 × 10^{−2}  346
150 6.7 × 10^{9} 3.4 × 10^{8} 2.6 × 10^{0}  20
150 1.8 × 10^{9} 9.7 × 10^{7} 8.3 × 10^{−1}  4
200 6.3 × 10^{9} 2.1 × 10^{8} 1.5 × 10^{0}  29
200 6.2 × 10^{9} 1.0 × 10^{9} 9.4 × 10^{0}  5
300 2.1 × 10^{10} 1.6 × 10^{9} 1.4 × 10^{1}  12
350 3.8 × 10^{10} 9.0 × 10^{9} 5.7 × 10^{1}  4
400 7.6 × 10^{10} 1.5 × 10^{10} 9.5 × 10^{1}  5

Table 2: Summary of numerical results for the Frankfurt stock exchange test set from [23]. Flops/pt.: floating point operations per computed efficient point. Sec/pt.: seconds of runtime per computed efficient point.

References

[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] Jens Albrecht. Vorausschauende optimale Steuer- und Regelstrategien zur Verbesserung der Kraftwerksführung. PhD thesis, Fakultät für Elektrotechnik, Universität Dortmund, Universität Dortmund, 44221 Dortmund, Germany, 1997. Published in Fortschrittsberichte VDI, Series 8, No. 616.

[3] U. Bagchi. Simultaneous minimization of mean and variation of flow time and waiting time in single machine systems. Operations Research, 37:118–125, 1989.

[4] J. E. Beasley. Or-library: distributing test problems by electronic mail. Journal of the Operational Research Society, 41:1069–1072, 1990. Test problems

electronically available at http://www.ms.ic.ac.uk/info.html.

[5] 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.

[6] A. B. Berkelaar. Complementarity and Interior Point Methodology. Master’s thesis, Department of Statistics, Probability and Operations Research, Faculty of Technical Mathematics and Computer Science, Delft University of Technology, Delft, The Netherlands, October 31 1995.

[7] Arjan B. Berkelaar, Benjamin Jansen, Kees Roos, and Tamás Terlaky. An Interior-Point Approach to Parametric Convex Quadratic Programming. Working paper, Erasmus University Rotterdam, Erasmus University Rotterdam, Rotterdam, The Netherlands, February 29 1996. Revised version from December 23, 1997.

[8] E[milio] Carrizosa, E. Conde, M[anuel] Munoz, and J[usto] Puerto. Planar point-objective location problems with nonconvex constraints: a geometrical construction. Journal of Global Optimization, 6:77–86, January 1995.

[9] 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.

[10] 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 [11].

[11] 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.

[12] 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.

[13] 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.

[14] H. Eschenauer, J. Koski, and A. Osyczka. Multicriteria Design Optimization. Springer, Berlin, 1990.

[15] G. Evans. Overview of techniques for solving multiobjective mathematical programs. Management Science, 30:1268–1282, 1984.

[16] 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/.

[17] Jörg Fliege. Approximation Techniques for the Set of Efficient Points. Professorial dissertation, Fachbereich Mathematik, Universität Dortmund, 44221 Dortmund, Germany, January 31 2001.

[18] Jörg Fliege. OLAF — a general modeling system to evaluate and optimize the location of an air polluting facility. ORSpektrum, 23(1):117–136, 2001.

[19] Jörg Fliege, H. Kiendl, and Daniel Schauten. A multicriteria model for power plant optimization. Technical report, Fachbereich Mathematik, Universität Dortmund, 44221 Dortmund, Germany, 2001. In preparation.

[20] Arthur M. Geoffrion. Proper Efficiency and the Theory of Vector Maximization. Journal of Optimization Theory and Applications, 22:618–630, 1968.

[21] 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.

[22] 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.

[23] 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.

[24] Claus Hillermeier. Nonlinear Multiobjective Optimization: a Generalized Homotopy Approach, volume 25 of International Series of Numerical Mathematics. Birkäuser Verlag, Basel — Boston — Berlin, 2001.

[25] 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.

[26] Johannes Jahn. Mathematical Vector Optimization in Partially Ordered Linear Spaces. Verlag Peter Lang, Frankfurt, 1986.

[27] A. Jüschke, J. Jahn, and A. Kirsch. A bicriterial optimization problem of antenna design. Computational Optimization and Applications, 7:261–276, 1997.

[28] 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.

[29] Harry M. Markowitz. Mean-Variance Analysis in Portfolio Choice and Capital Markets. Basil Blackwell Ltd, Oxford, 1987.

[30] H[arry] M. Markowitz. Portfolio selection, efficient diversification of investments. Basil Blackwell Ltd, Oxford, 2nd edition, 1991.

[31] 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.

[32] Markus Rudolf. Algorithms for Portfolio Optimization and Portfolio Insurance. Paul Haupt, Bern — Stuttgart — Wien, 1994.

[33] Yoshikazu Sawaragi, Hirotaka Nakayama, and Tetsuzo Tanino. Theory of Multiobjective Optimization, volume 176 of Mathematics in Science and Engineering. Academic Press, Inc., Orlando, San Diego, New York, London, 1985.

[34] Manuel A. Nunez und Robert M. Freund. Condition measures and properties of the central trajectory of a linear program. Mathematical Programming, 83:1–28, 1998.

[35] D[ouglas] J[ohn] White. Epsilon-dominating solutions in mean-variance portfolio analysis. European Journal of Operational Research, 105:457–466, 1998.

[36] Stephen J. Wright. Primal-dual interior-point methods. SIAM, Philadelphia, 1997.

[37] E. Alper Yıldırım and Michael J. Todd. Sensitivity analysis in linear programming and semidefinite programming using interior-point methods. Technical Report 1253, School of operations research and industrial engieering, Cornell University, Ithaca, NY 14853-3801, 1999. Published as [38].

[38] E. Alper Yıldırım and Michael J. Todd. Sensitivity analysis in linear programming and semidefinite programming using interior-point method. Mathematical Programming, 90:229–261, 2001.

[39] E. Alper Yıldırım 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.