Nonconvex min–max fractional quadratic problems under quadratic constraints: copositive relaxations

In this paper we address a min–max problem of fractional quadratic (not necessarily convex) over linear functions on a feasible set described by linear and (not necessarily convex) quadratic functions. We propose a conic reformulation on the cone of completely positive matrices. By relaxation, a doubly nonnegative conic formulation is used to provide lower bounds with evidence of very small gaps. It is known that in many solvers using Branch and Bound the optimal solution is obtained in early stages and a heavy computational price is paid in the next iterations to obtain the optimality certiﬁcate. To reduce this effort tight lower bounds are crucial. We will show empirical evidence that lower bounds provided by the copositive relaxation are able to substantially speed up a well known solver in obtaining the optimality certiﬁcate.


Introduction and motivation
In optimization, the objective function of the underlying mathematical model represents a certain performance measure to be optimized. In many situations this performance is described as the ratio of two functions to accurately reflect the trade-off between two, sometimes competing, aspects of the system, like profit versus the number of employees or production versus consumption, mean value and variance, cost versus time or volume. Problems of this kind arise in many different activities such as financial planning (where the ratio of debt/equity must be considered), in production planning (inventory/sales or revenue/employee), health care (cost/patient), blending problems (income/quantity of raw material), in the domains of optics, engineering, or information retrieval, to name just a few. More application examples can be found in the first chapter of [25]. Sometimes the fractional relation is ignored for the sake of model simplicity and replaced by a parameterized linear function. An alternative approach introduces bounds on one of the functions as a constraint in the model. However, recent interest has focused on more realistic mathematical models, and increased research effort focuses on advanced methodologies to directly address the problem of ratio optimization, also known as fractional programming.
Consider the so called single-ratio fractional problem, where ⊆ R n is a bounded convex set and it is assumed that g(x) > 0 for x ∈ . If f and g are affine functions and is a polyhedron, problem (1.1) is called a linear fractional optimization problem. If the functions f and g are quadratic (one of them can be affine) and is a polyhedron, we will call this problem as a quadratic fractional optimization problem on a polyhedron. If is not a polyhedron it will just be referred as a quadratic fractional optimization problem.
Fractional problems (1.1) are nonconvex in general. Under convexity and concavity assumptions for f and g respectively, and nonnegativity of f if g is not affine, (1.1) is termed a convex-concave fractional program [22], and the objective function is strictly quasi-convex. Any local solution to a convex-concave fractional program is also a global solution, as it happens for general convex programs, and in addition a solution of the Karush-Kuhn-Tucker optimality conditions is a global solution if the numerator and denominator functions are differentiable on an open set containing . But the convexity and concavity assumptions are restrictive, and in many practical problems they are not valid.
As a generalization of (1.1) we can find the generalized fractional program or min-max fractional optimization problem min x∈ max y∈U f (x, y) g(x, y) , (1.2) and as a special case the discrete case, where we have a finite set I = {1, . . . , m} instead of U : In the sequel we will consider the discrete case and assume that: • for all i ∈ I , f i and g i are continuous on an open set˜ ⊆ R n containing ; (1.5) • g i (x) > 0 for all x ∈ and all i ∈ I ; (1.6) • f i (x) > 0 for all x ∈ and all i ∈ I . (1.7) Note that condition (1.6) aims at ensuring that g i (x) = 0 and so the problem is wellposed. The positivity condition is not restrictive because we can multiply both numerator and denominator by (−1) ensuring this condition even if g i (x) < 0. We can assume (1.7) without loss of generality since we can always consider an equivalent problem adding L to g i (x) in the objective function for i ∈ I with L > 0 large enough such that the new numerator In this paper we are particularly interested in problems (1.3) since they are relevant for performing worst-case analysis as mentioned in [15,20]. Indeed, suppose that the uncertainty set U ⊂ R p is a polytope with known (few) vertices v 1 , . . . v m , so U = conv {v i : i ∈ I }, and further assume that the objective f (x, y)/g(x, y) depends in a quasiconvex way on y for any fixed x ∈ . Now putting f i (x) = f (x, v i ) and g i (x) = g(x, v i ), it is obvious that problems (1.2) and (1.3) are the same, so that the robust formulation can be phrased in terms of finitely many scenarios which do not affect the set representing precisely known (not uncertain, global) constraints. Closely related applications arise from goal programming and multicriteria optimization when Chebychev's norm is used to evaluate the goals or criteria, leading to The equivalence of (1.3) and (1.8) follows as above, or else directly from the fact that the maximum of m numbers is equal to the maximum of their convex combinations [8,13]. Problem (1.8) is particularly interesting when pooling of the rival objective functions is considered, in which case the components w i are the pooling weights. Then the min-max solution provides a robust pooling by suggesting a strategy that does not depend on the prevailing scenario.
The solution of (1.3) allows the decision maker to compute the best strategy under the worst-case scenario. In practical terms, solving the min-max problem ensures that the solution value will not deteriorate whichever scenario turns out to be the true one. The real solution evaluation will be at least as good as the min-max value. In this context the above mentioned applications of (1.1) can be directly extended to (1.3) in cases where there is some uncertainty related to outside factors. The maximum operator can be avoided by introducing another variable; indeed, problem (1.3) is equivalent to the following problem We will investigate (1.3) in case of quadratic functions f i and affine g i . Problems of this type have applications in maximal predictability portfolio [11], cluster analysis [19], transportation problems with several objective functions [24], and investment allocation problems [27], to name but a few. This paper is organized in five sections. After this introduction, Sect. 2 presents a review on some methods for addressing general min-max problems. In Sect. 3 the conic reformulation and a tractable relaxation is constructed. This section also incorporates a discussion regarding the extension of these results to the quadratic over quadratic case. Sect. 4 reports computational experience, followed by concluding remarks in the last section.

Short review of some traditional methods
Before we proceed we will present in this section two traditional approaches for addressing the min-max problem. Both suffer from the same drawback, namely for general matrices the subproblems are nonconvex, so most of them cannot be solved quickly and the first approach does not even guarantee optimality. Throughout this section, we assume for simplicity of exposition that is compact and all involved functions are continuous.

The minimum-regret approach
We start by introducing a naive minimax strategy based on the regret concept [20]. First, an optimal solution (which generically is unique) for each individual objective function is obtained: Next, the impact of solution x * k , k ∈ I , on the other objective functions is investigated, tantamount to the loss incurred by adopting x * k if a scenario different from k occurs: .
The best strategy x * will be the one corresponding to the solution that would be the least damaging if a different scenario replaces the one for which the solution is optimal, i.e., select i ∈ I such that α i = min k∈I α k and put x * = x * i .

A generalization of the Dinkelbach approach
Denote by In [9] a method for solving (1.3) based on the Dinkelbach approach for a single ratio is presented. considering a parameterized version of problem (2.2): and a solution of the equation F(λ) = 0 is sought. It is known that F(λ) has some important properties, namely: Algorithm 1 describes the adaptation of the single ratio Dinkelbach method to multiple ratios. The difference consists mainly in the update of the parameter λ in each iteration, steps (1) and (2). The convergence of this method is linear, in contrast to the superlinear convergence in case of the Dinkelbach method for a single ratio (which follows by analogy to the Newton method). A better convergence is obtained by modifying the update of x k to by approximating a Newton step, and thus increasing the convergence rate to superlinearity. This method is also equivalent to the partial linearization procedure in [3]. See [9].
Data: , f i , g i for i ∈ I Result: Solution of problem (2.2) x 0 ∈ ; STOP=false; k ← 1; = 10 −6 ; STOP=true; output the optimal solution x k and optimal value λ k ; else Algorithm 1: Dinkelbach for multiple ratios

A copositive approach
Let us recall the problem (1.3) under investigation .
In this work we are considering with A a k × n matrix while a q ∈ R n and A q ∈ M n for all q∈ [1 : p]. We write 2r i in the linear term later in (3.2), to keep the analogy with the numerator and to allow for a simpler notation when discussing extensions to the quadratic over quadratic case. Observe that we have implicitly introduced slack variables for the linear constraints but keep quadratic inequalities as constraints. No definiteness assumptions are made on A q . Likewise, in the sequel we will make no assumptions about matrices Q i in (3.1), other than symmetry. Our proposal is to use copositive optimization to construct good lower bounds that can be incorporated in a global optimization method for solving this class of min-max problems. Many solvers return a feasible solution at early stages of the branch and bound tree which turns out to be optimal after a following, in many cases long, sequence of iterations, incurring heavy computational burden for obtaining an optimality certificate. Good lower bounds are crucial at this stage. We will show how the lower bounds provided by the copositive relaxation can speed up BARON to obtain this optimality certificate earlier. But before that we briefly recapitulate some essential concepts about copositive optimization.

Quick facts about copositive optimization
A conic optimization problem in the space S n of symmetric matrices of order n consists in optimizing a linear function with linear constraints over a cone K ⊆ S n : is also a conic optimization problem of the form which involves the dual cone As usual, we drop sign constraints u i ≥ 0 on the dual variables in case of linear equality constraints A i • X = b i on the primal side.
Well-investigated special matrix cones include the semidefinite cone S + n ⊂ S n of all positive-semidefinite symmetric n × n matrices (i.e. those X with no negative eigenvalue, denoted by X O), and the cone N n ⊂ S n of symmetric n × n matrices X with no negative entry (denoted by X ≥ O). These two cones are selfdual, S + n = S + * n and N n = N * n . Combination of both properties leads to the cone of so-called doubly nonnegative matrices The dual cone D * n = S + n + N n , so D n is no longer self-dual. While all of aforementioned cones are tractable in a certain sense, there is an important intractable subcone of D n , namely the cone of all completely positive symmetric n × n matrices, with its dual cone, that of copositive matrices which are symmetric of order n: We have CP 4 = D 4 but CP n D n for n ≥ 5. A copositive optimization problem is a conic problem (3.4) with K = COP n , and most authors include in the definition of copositive optimization also the case K = CP n . This problem class has recognized merits since it admits reformulations of hard optimization problems, such as continuous nonconvex quadratic [4,5,18], mixed-integer quadratic [7], continuous [2] and mixed-integer [1] fractional quadratic problems. Basically, by lifting quadratic expressions, all constraints can be linearized, pushing the hardness entirely into the cone description. Indeed, checking that a matrix is in CP n is NP-hard [10], and it is co-NP-complete to check whether a matrix is copositive [14]. Nevertheless it is an important feature of completely positive formulations that relaxations giving tight lower bounds can be obtained by replacing CP n by D n , for instance. Even tighter relaxation bounds can be obtained by the so-called approximation hierarchies which replace CP n by other tractable cones (polyhedral or expressible by linear matrix inequalities). These approximation hierarchies can be constructed using simplices [6,23], sums-of-squares [16] or polynomials [4,17].

Copositive reformulation
We now present a reformulation based on copositive optimization of problem (1.3) written now in a more descriptive way to facilitate further reading.
Using a Shor lifting and considering y = 1 , Next we square the homogenized linear constraints: Likewise, homogenize the quadratic constraints by introducing, for all q∈ [1 : p], we arrive at The linear constraints (3.11) are homogenized as inequalities instead of equivalent equalities (given that A 0 is pd) to simplify notation, defining all the constraints, quadratic and linear, with the same structure.
Next, borrowing from (1.9), we introduce another variable v ∈ R and using hypothesis (1.6), we obtain we obtain Now for X = zz , we obtain the following reformulation of (3.14), more precisely, of (λ * ) 2 , in view of z * n+1 ≥ 0: where CP rk 1 n+2 := {X ∈ CP n+2 : rank X = 1}. Dropping the rank constraint leads to the completely positive relaxation where E k • X = X kk . Weak duality and rank relaxation yield immediately the relations γ * C O P ≤ γ * C P ≤ γ * = (λ * ) 2 . If an optimal solution of (3.16) has rank one then we have γ * C P = γ * and vice versa.
Further, tractable relaxations of above conic problems may use approximation hierarchies for CP n+2 , for instance by replacing CP rk1 n+2 by D n+2 , the cone of doubly nonnegative matrices. In this case we obtain the following problem: Example 3.1 Consider the following example: Upon 0 ≤ x 2 = 1 − x 1 2 , the quadratic constraints reduces to x 1 ≥ 0 (which is redundant) and we obtain the following one-dimensional problem min The next example considers a convex instance where the objective h(x) is smooth at the optimal solution. γ * C O P = γ * C P = √ γ * = λ * , in this example. We will specify a more general condition for zero relaxation gap below.

Theorem 3.3 (a) For any
consider the following symmetric matrix of order n + 2:

19)
where If S u is copositive, then √ u 0 ≤ λ * , so we get a valid lower bound for the problem (1.3) with affine-linear denominators g i (x). (b) Now suppose that x * ∈ and put √ u 0 = h(x * ). If for this value of u 0 , above S u is copositive, then x * solves the problem (1.3) with affine-linear denominators g i (x), and both the conic duality and relaxation gaps are zero: Proof (a) The matrix S u as defined in (3.19) coincides with the slack matrix of the dual problem (3.17), as one can easily verify. By assumption, this S u is feasible to (3.17), hence we obtain , and all assertions are established.

Some considerations about the constrained min-max quadratic over quadratic case
Before presenting the computational experience we discuss some limitations regarding extensions of the previous work, on a conic reformulation, to the quadratic over quadratic case, when: The quadratic over affine case is obtained considering T i = O for i ∈ I . Knowing that it is possible to replace a quadratic expression in the objective function by an additional variable introducing this relation as a constraint of the problem, we immediately looked for the possibility of applying the same conic reformulation to this more general case but we faced strong limitations regarding the conic relaxation. First we show that that this case can be converted to a quadratic over affine problem, as (3.10) with the introduction of additional variables, x n+i replacing the quadratic expressions x T i x in g i (x) and adding the quadratic constraints x T i x ≥ x n+i for i ∈ I . A problem with an affine denominator is obtained subject to m + p + 1 quadratic constraints plus (linear) sign constraints. It is worth noting that additional quadratic constraints should not pose further difficulties concerning the conic reformulation. An equivalent problem is obtained as follows: denoting by O zero matrices of various but suitable size where we abbreviated with e i the ith column of the identity matrix I m of order m.

Example 3.4 Consider the following example
with = x ∈ R 2 + : x 1 + 2x 2 = 2 . Upon replacing x 2 = 1 − x 1 2 we obtain the following one-dimensional problem min Figure 3 shows both quadratic ratios, and the optimal value 0.4381 over [0, 2] attained at the optimal solution x * = [0.9122, 0.5439] is designated by a circle. Applying above reasoning, problem (3.23) can be rewritten as When observing in detail the structure of the constraints of the conic optimization relaxation for the quadratic over quadratic case we have (abbreviating Z = X n+1:n+m,n+1:n+m ) ⎡ The occurrence of a sub-matrix of zeros (in the framed box) in constraints (3.24)-(3.26) contributes to a severe deterioration of the quality of the lower bounds. Further research on this topic is left for future investigation.

Computational experience
The numerical experiments were performed on a PC, Intel(R) Core(TM) i7-2640M, 2.80 Ghz,400 GB RAM. The software Matlab 2013Ra was used to run the global optimization solver BARON [21] and SDPT3(4.0)/Octave [26], a Matlab package for solving convex optimization problems involving linear equations and inequalities, secondorder cone constraints, and semidefinite constraints (linear matrix inequalities). The interface YALMIP [12] was used to call SDPT3.
To study the empirical quality of the lower bounds provided by the relaxation (3.18) we used the global optimization solver BARON to obtain, when possible, the optimal solution to define the gaps and to provide a lower bound for comparison. It would be useless to use the last lower bound provided by BARON since in most problems it would coincide with the optimal value and it is the result of a chain of different methods (which are unknown to the user). To compare with the first lower bound of BARON would be a fair choice, but even better is to compare the performance of BARON with itself when the lower bound of relaxation δ * is provided to this solver.
The formulation defined for BARON input was the following: The software BARON does not accept, as input, an initial lower bound to start the search, but in this case, since the objective function is a single variable v and BARON accepts box bounds (lower and upper) for variables it was in fact possible to indirectly impose a lower bound by introducing the lower box bound v l for v. This allowed to study the quality of the lower bound by solving (3.18), obtaining as optimal value X * n+1,n+1 and running BARON twice, with v l = −∞ and v l = X * n+1,n+1 . The limit time ( in seconds) for the first run of BARON was set to MaxTime = 250 and the gap tolerance EpsA = 10 −4 . On the second run of BARON, the time limit was set differently to each instance to 250 minus the run time of SDPT3 for that instance. The settings for SDPT3 were sdpt3.maxit = 500 for the maximal number of iterations and sdpt3.gaptol = 10 −2 for the gap tolerance.

Empirical study on the quality of lower bounds
In the analysis of the results we focus to answer the question whether the relaxation is more efficient than the battery of techniques used in BARON to generate lower bounds for this problem, using a first run of BARON with v ≥ −∞ (RUN1) and a second run with v ≥ X * n+1,n+1 (RUN2). We compare to the result of the SDPT3 solver for solving (3.18), designated by RUN0. Table 1 presents results from RUN0, RUN1 and RUN2 where the columns represent: • Problem-Problem name 'INSTaMbNc' where 'a' is the instance number (0-9) for the same specifications: 'b' number of ratios and 'c' number of variables.  • GAP0-100 BARON UB−YLB BARON UB for RUN0 • Time-CPU time of RUN0 Looking at Table 1, merits of the lower bounds provided by (3.18) are clear and more pronounced for larger instances. First we should observe the resulting gaps and see that for the smaller problems the lower bound coincides with the optimal value in most instances. As the size increases, it is evident that the lower bound from (3.18) outperforms BARON's lower bounding techniques. In some cases (INST8M5N25, INST8M5N50) only by providing the lower bound from (3.18) it was possible to obtain the optimal solution. The caveats occur for some numerical problems (E2 = 4) encountered in some instances, and the lower bound cannot be trusted. In fact significantly negative gap values (< −10 −2 ) represent these cases while small negative gaps are just the result of error propagation. Figure 4 gives a graphical representation of these gaps. It is clear that the line for GAP1 is mostly above that of GAP2 and GAP0 due to usage of (3.18), and that the line of GAP2 lies close to that of GAP0, meaning that the proposed lower bound represents an important contribution. The decrease in each instance from GAP1 to GAP2 is the direct result of the a priori knowledge of this lower bound.
The next observation is related to CPU time in columns BTime and Time. Most of the time the best upper bound is achieved at an early stage and so BARON puts considerably larger effort in improving the lower bound. So, even discounting the fact that BARON time is not only related to lower bound calculations, Table 1 shows that it clearly pays to have an initial lower bound calculated by (3.18).
In just a few cases (Btime < 0) it seems that it does not compensate the a priori knowledge of the lower bound given by (3.18), even if the bound in the first node is really weak. We must bear in mind that it is not clear at which stage of the procedure (if at all) BARON interprets the lower bound box constraint for variable v also as a lower bound for the optimal value. Nevertheless the overall conclusion is that it is a benefit and sometimes even crucial to have the lower bound provided by (3.18).

Conclusion
In this paper we make a contribution to extend the class of hard optimization problems for which a conic reformulation, based on the cone of copositive or completely positive matrices, is known. We address min-max problems of fractional quadratic (not necessarily convex) over affine functions on a feasible set described by linear and (not necessarily convex) quadratic functions. Relaxations of copositive reformulations have proved to produce tight lower bounds and empirical evidence obtained by the computational experience on this problem corroborates this statement. Relaxation of the completely positive cone by the doubly nonnegative cone led to a tractable relaxation problem providing lower bounds with very small gaps, enabling us to globally solve a difficult optimization problem. It is known that in branch and bound methods, very good, sometimes even optimal, solutions are obtained at early stages of traversing the problem tree. The effort to prove optimality or to state a safe termination criterion for the search depends on closing the gap between the lower and upper bound, so the quality of the first are crucial. In this paper we present a lower bounding procedure that justifies once more the high reputation of relaxations gained from conic reformulations.