Computational experiments with a linear approximation of second-order cone optimization

Image Technical Report 0001

François Glineur^{∗}

Service de Mathématique et de Recherche Opérationnelle, Faculté Polytechnique de Mons, Rue de Houdain, 9, B-7000 Mons, Belgium,

                         Francois.Glineur@fpms.ac.be.

January 2000, revised November 2000

∗The author is a research fellow supported by the F.N.R.S. (Belgian National Fund for Scientific Research). This research was partly carried out during a visit to the SSOR department of the Delft University of Technology supported by a grant of the Communauté Française de Belgique.

Abstract

In this article, we present and improve a polyhedral approximation of the second-order cone due to Ben-Tal and Nemirovski [BTN98]. We also discuss several ways of reducing the size of this approximation. This construction allows us to approximate second-order cone optimization problems with linear optimization.

We implement this scheme and conduct computational experiments dealing with two classes of second-order cone problems: the first one involves truss-topology design and uses a large number of second-order cones with relatively small dimensions, while the second one models convex quadratic optimization problems with a single large second-order cone.

Keywords. second-order cone optimization, approximation, linear optimization, quadratic optimization.

Contents

1 Introduction 3

2 Approximating second-order cone optimization 4 2.1 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 A first approximation of 𝕃^{2} . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.4 A better approximation of 𝕃^{2} . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.5 Reducing the approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.6 An approximation of 𝕃^{n} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.7 Optimizing the approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.8 An approximation of second-order cones optimization . . . . . . . . . . . . . 16 2.9 Accuracy of the approximation . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3 Computational experiments 19 3.1 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2 Truss-topology design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3 Quadratic optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4 Concluding remarks 31

1 Introduction

Many useful classes of optimization problems one encounters both in theory and in practice are convex, such as linear optimization, convex quadratic optimization, semidefinite optimization and geometric optimization. This fundamental characteristic guarantees some nice properties: every local minima is a global minima, existence of strongly related dual problem and availability of efficient algorithms to solve these problems.

Conic optimization is a powerful setting that relies on convex cones to formulate convex problems: the standard conic primal-dual pair is the following [SW70, Stu99a]

inf_{x∈ℝ^{n}} c^{T} x s.t. Ax = b and x ∈ C  (1.1)
sup_{y∈ℝ^{m},x^{∗}∈ℝ^{n}} b^{T} y s.t. A^{T} y + x^{∗} = c and x^{∗} ∈ C^{∗} ,  (1.2)

where x and (y, x^{∗}) are the primal and dual variables, A is a m × n matrix, b and c are m and n-dimensional column vectors, C ⊆ ℝ^{n} is a closed convex cone and C^{∗} ⊆ ℝ^{n} is its dual cone, defined by C^{∗} = {x^{∗} ∈ ℝ^{n} | x^{T} x^{∗} ≥ 0 ∀x ∈ C}.

Different types of convex cones lead to different classes of problems: for example, linear optimization uses the nonnegative orthant ℝ_{+}^{n} while semidefinite programming relies on the set of positive semidefinite matrices 𝕊_{+}^{n}. In this article, we will focus the second-order cone, also known as Lorentz cone or ice-cream cone, which leads to second-order cone optimization. It is defined as follows:

Definition 1.1. The second order cone 𝕃^{n} is the subset of ℝ^{n+1} defined by

𝕃^{n} = {(r, x) ∈ ℝ × ℝ^{n} | ∥x∥ ≤ r} ,

where ∥·∥ denotes the usual Euclidean norm on ℝ^{n}.

It is indeed straightforward to check that this is a closed convex cone (it is in fact the epigraph of the Euclidean norm). Another interesting property of 𝕃^{n} is the fact that it is self-dual, i.e. (𝕃^{n})^{∗} = 𝕃^{n}.

The standard second-order cone problems are based on the cartesian product of several second-order cones, which can be formalized using r constants n_{k} ∈ ℕ, 1 ≤ k ≤ r such that ∑_{k=1}^{r}(n_{k} +1) = n and defining C = 𝕃^{n_{1}} ×𝕃^{n_{2}} ×⋯×𝕃^{n_{r}}. This set is obviously also a self-dual closed convex cone, which allows us to rewrite problems (1.1) and (1.2) as

inf c^{T} x s.t. Ax = b and x^{k} ∈ 𝕃^{n_{k}} ∀k = 1, 2, …, r  (1.3)
sup b^{T} y s.t. A^{T} y + x^{∗} = c and x^{k∗} ∈ 𝕃^{n_{k}} ∀k = 1, 2, …, r ,  (1.4)

where vectors x and x^{∗} have been split into r subvectors (x^{1}, x^{2}, …, x^{r}) and (x^{1∗}, x^{2∗}, …, x^{r∗}) with x^{k∗}, x^{k} ∈ ℝ^{n_{k}+1} for all k = 1, …, r. It is usually more practical to pivot out variables x^{∗} in the dual problem (1.4), i.e. write them as a function of vector y. Splitting matrix A into (A^{1}, A^{2}, …, A^{r}) with A^{k} ∈ ℝ^{m×(n_{k}+1)} and vector c into (c^{1}, c^{2}, …, c^{r}) with c^{k} ∈ ℝ^{n_{k}+1}, we have x^{k∗} = c^{k} − A^{k}^{T} y. The last step is to isolate the first column in A^{k} and the first component in c^{k}, i.e. letting A^{k} = (f^{k}, G^{k}) with f^{k} ∈ ℝ^{m} and G^{k} ∈ ℝ^{m×n_{k}} and c^{k} = (d^{k}, h^{k}) with d^{k} ∈ R and h^{k} ∈ ℝ^{n_{k}}, we can rewrite the dual problem (1.4) as

max b^{T} y s.t. ∥G^{k}^{T} y + h^{k}∥ ≤ f^{k}^{T} y + d^{k} ∀k = 1, 2, …, r ,

which is more convenient to formulate real-world problems.

Second-order optimization admits many different well-known classes of optimization problems as special cases, such as linear optimization, linearly and quadratically constrained convex quadratic optimization, robust linear optimization, matrix-fractional problems and problems with hyperbolic constraints (see the survey [LVBL98]). Applications arise in various fields such as engineering (antenna array design, finite impulse response filter design, truss design) and finance (portfolio optimization), see again [LVBL98].

From the computational point of view, second-order cone optimization is a relatively young field if compared to linear and quadratic optimization (for example, the leading commercial linear and quadratic solvers do not yet offer the option of solving second-order cone optimization problems). This observation led Ben-Tal and Nemirovski to develop an interesting alternative approach to solving second-order cone problems: they show in [BTN98] that it is possible to write a polyhedral approximation of the second order cone 𝕃^{n} with a prescribed accuracy ϵ using a number of variables and constraints that is polynomial in n and log \frac{1}{ϵ}. This implies that second-order cone optimization problems can be approximated with an arbitrarily prescribed accuracy by linear optimization problems using this polyhedral approximation.

This potentially allows the approximate resolution of large-scale second-order cone problems using state of the art linear solvers, capable of handling problems with hundreds of thousands of variables and constraints.

This paper is organized as follows: Section 2 presents a polyhedral approximation of the second-order cone. This construction relies on a decomposition scheme based on threedimensional second order cones. We present first an efficient way to approximate these cones and then show how to combine them in order to approximate a second-order cone of higher dimension, which ultimately gives a method to approximate any second-order cone optimization problem with a linear problem. Section 3 reports our computational experiments with this scheme. After a presentation of our implementation and some related issues, we describe two classes of second-order problems: truss-topology design problems and convex quadratic optimization problems. We present and discuss the results of our computational experiments, highlighting when necessary the particular features of each class of problems (guaranteed accuracy, alternative formulations). We conclude this article with a few remarks and suggestions for further research.

2 Approximating second-order cone optimization

In this section, we present a polyhedral approximation of the second-order cone 𝕃^{n} which allows us to derive a linearizing scheme for second-order cone optimization. It is a variation of the construction of Ben-Tal and Nemirovski that features slightly better properties.

2.1 Principle

The principle that lies behind their approximation is twofold:

a. Decomposition. Since the Lorentz cone 𝕃^{n} is a n + 1-dimensional subset, any circumscribed polyhedral cone around 𝕃^{n} is bound to have its number of facets growing exponentially with the dimension n, i.e. will need an exponential number of linear inequalities to be defined. The remedy is to decompose the second-order cone into a polynomial number of smaller second-order cones with fixed dimension, for which a good polyhedral approximation can be found. In the present case, 𝕃^{n} can be decomposed into n − 1 three-dimensional second-order cones 𝕃^{2}, at the price of introducing n − 2 additional variables (see Section 2.2).

b. Projection. Even the three-dimensional second-order cone 𝕃^{2} is not too easy to approximate: the most obvious way to proceed, a regular circumscribed polyhedral cone, requires hundreds of inequalities even for an approximation with modest accuracy (see Section 2.3). The key idea to lower the number of inequalities is to introduce several additional variables, i.e. lift the approximating polyhedron into a higher dimensional space and consider its projection onto a (n + 1)-dimensional subspace as the approximation of 𝕃^{n} (see Section 2.4).

To summarize, the introduction of a certain number of additional variables, combined with a projection, can be traded against a much lower number of inequality constraints defining the polyhedron. We first concentrate on the decomposition of 𝕃^{n} into smaller second-order cones.

2.2 Decomposition

Let us start with the following equivalent definition of 𝕃^{n}

𝕃^{n} = {(r, x_{1}, x_{2}, …, x_{n}) ∈ ℝ_{+} × ℝ^{n} | ∑_{i=1}^{n} x_{i}^{2} ≤ r^{2}} .

Introducing a vector of ⌊\frac{n}{2}⌋ additional variables y = (y_{1}, y_{2}, …, y_{⌊ \frac{n}{2} ⌋}), we consider the set 𝕃^{n}′ defined by

\{(r, x, y) ∈ ℝ_{+} × ℝ^{n+⌊ \frac{n}{2} ⌋} | x_{2i−1}^{2} + x_{2i}^{2} ≤ y_{i}^{2}, 1 ≤ i ≤ ⌊\frac{n}{2}⌋, \begin{cases}∑i=1^{\frac{n}{2}} y_{i}^{2} & ≤ r^{2} (n even) \\ ∑i=1^{\frac{n}{2}} y_{i}^{2} + x_{n}^{2} & ≤ r^{2} (n odd) \end{cases} \} .

It it straightforward to prove that the projection of this set on the subspace of its first n + 1 variables (r, x_{1}, …, x_{n}) is equal to 𝕃^{n}, i.e.

(r, x) ∈ 𝕃^{n} ⇔ ∃y ∈ ℝ^{⌊ \frac{n}{2} ⌋} s.t. (r, x, y) ∈ 𝕃^{n}′ .

It is also worth to point out that all the constraints defining 𝕃^{n}′ are second-order cone constraints, i.e. that 𝕃^{n}′ can also be written as

\{(r, x, y) ∈ ℝ_{+} × ℝ^{n+⌊ \frac{n}{2} ⌋} | (y_{i}, x_{2i−1}, x_{2i}) ∈ 𝕃^{2}, 1 ≤ i ≤ ⌊\frac{n}{2}⌋, \begin{cases}(r, y) & ∈ 𝕃^{⌈ \frac{n}{2} ⌉} (n even) \\ (r, y, x_{n}) & ∈ 𝕃^{⌈ \frac{n}{2} ⌉} (n odd) \end{cases}\} .

(2.1)

This means that 𝕃^{n} can be decomposed into ⌊\frac{n}{2}⌋ small 3-dimensional second-order cones and a single 𝕃^{⌈ \frac{n}{2} ⌉} second-order cone, at the price of introducing ⌊\frac{n}{2}⌋ auxiliary variables. This procedure can be applied recursively to the largest of the remaining second-order cone 𝕃^{⌈ \frac{n}{2} ⌉} until it also becomes equal to 𝕃^{2}.

It is not too difficult to see that there are in the final expression n − 1 second-order cones 𝕃^{2} and n − 2 additional y_{i} variables. Indeed, the addition of each small cone 𝕃^{2} reduces the size of the largest cone by one, since we remove two variables from this cone (the last two variables in 𝕃^{2}) but replace them with a single new variable (the first variable in 𝕃^{2}).

Since we start with this largest cone equal 𝕃^{n} and stop when its size is equal to 2, we need n − 2 small cones along with n − 2 auxiliary variables to reduce the cone to 𝕃^{2}. But this last 𝕃^{2} cone also has to be counted, which gives then a total number of cones equal to n − 1.

The existence of this decomposition implies that any second-order cone optimization problem can be transformed into a problem using only 3-dimensional second-order cones, using the construction above. We note however that strictly speaking, the resulting formulation is not a conic problem, since some variables belong to two different cones at the same time. It is nonetheless possible to add a extra variable for each shared variable, along with a constraint to make them equal on the feasible region, to convert this formulation into the strict conic format (1.1)–(1.2).

2.3 A first approximation of 𝕃^{2}

The previous section has shown that we can focus our attention on approximations of the 3-dimensional second-order. Moreover, it seems reasonable to require this approximation to be a cone itself too. Taking into account this additional assumption, we can take advantage of the homogeneity property of these cones to write

(r, x_{1}, x_{2}) ∈ 𝕃^{2} ⇔ (1, \frac{x_{1}}{r} , \frac{x_{2}}{r}) ∈ 𝕃^{2} ,  (2.2)

which basically means we can fix r = 1 and look for a polyhedral approximation of the resulting set

{x ∈ ℝ^{2} | (1, x) ∈ 𝕃^{2}} = {x ∈ ℝ^{2} | x_{1}^{2} + x_{2}^{2} ≤ 1} = B_{2}(1) ,

which is exactly the disc of radius one in ℝ^{2}. Any approximating polyhedron for B_{2}(2) will be then later straightforwardly converted into a polyhedral cone approximating 𝕃^{2}, using the additional homogenizing variable r.

At this point, we have to introduce a measure of the quality of our approximations. A natural choice for this measure is to state that a polyhedron P ⊆ ℝ^{2} is a ϵ-approximation of B_{2}(1) if and only we have the double inclusion B_{2}(1) ⊆ P ⊆ B_{2}(1 + ϵ), i.e. the polyhedron contains the unit disc but lies entirely within the disc of radius 1 + ϵ.

The most obvious approximation of the unit disc is the regular m-polyhedron P_{m}, which is described by m linear inequalities. We have the following theorem:

Theorem 2.1. The regular polyhedron with m sides is an approximation of the unit disc B_{2}(1) with accuracy ϵ = cos(\frac{π}{m})^{−1} − 1.

Proof. The proof is quite straightforward: looking at Figure 1 (which represents the case m = 8), we see that angle ∠AOM is equal to \frac{π}{m} and thus that |OA| cos(\frac{π}{m}) = |OM| = 1. Our measure of quality is then equal to ϵ = |OA| − 1 = cos(\frac{π}{m})^{−1} − 1, as announced.

This result is not very satisfying: since cos(x)^{−1} ≈ 1 + \frac{x^{2}}{2} when x is small, we have that ϵ ≈ \frac{π^{2}}{2m^{2}} when m is large, which means that doubling the number of inequalities only divides the accuracy by four. For example, approximating B_{2}(1) with the relatively modest accuracy 10^{−4} would already take a 223-sided polyhedron, i.e. more than 200 linear inequalities.

Figure 1: Approximating B_{2}(1) with a regular octagon.

2.4 A better approximation of 𝕃^{2}

As outlined in Section 2.1, the key idea introduced by Ben-Tal and Nemirovski to obtain a better polyhedral approximation is to consider the projection of a polyhedron belonging to a higher dimensional space. The construction we are going to present here is a variation of the one described in [BTN98], featuring slightly better parameters and a more transparent proof.

Let us introduce an integer parameter k ≥ 2 and consider the set D_{k} ⊆ ℝ^{2k+2} defined as

\{(α_{0}, …, α_{k}, β_{0}, …, β_{k}) ∈ ℝ^{2k+2} | \{ \begin{matrix}α_{i+1} & = & α_{i} cos \frac{π}{2^{i}} & + & β_{i} sin \frac{π}{2^{i}} \\ β_{i+1} & ≥ & β_{i} cos \frac{π}{2^{i}} & − & α_{i} sin \frac{π}{2^{i}} \\ −β_{i+1} & ≤ & β_{i} cos \frac{π}{2^{i}} & − & α_{i} sin \frac{π}{2^{i}} \\ 1 & = & α_{k} cos \frac{π}{2^{k}} & + & β_{k} sin \frac{π}{2^{k}} \end{matrix} ∀ 0 ≤ i < k\} .

This set is obviously a polyhedron^{1}, since its defining constraints consist in k + 1 linear equalities and 2k inequalities. The following theorem gives some insight about the structure of this set.

Theorem 2.2. The projection of the set D_{k} on the subspace of its two variables (α_{0}, β_{0}) is equal to the regular 2^{k}-sided polyhedron, i.e. we have

(α_{0}, β_{0}) ∈ P_{2^{k}} ⇔ ∃(α_{1}, …, α_{k}, β_{1}, …, β_{k}) ∈ ℝ^{2k} | (α_{0}, …, α_{k}, β_{0}, …, β_{k}) ∈ D_{k} .

Proof. To fix ideas, we are going to present some figures corresponding to the case k = 3, but our reasoning will of course be valid for all k ≥ 2. Looking at Figure 1, which depicts P_{2^{3}}, we see that the last equality in the definition of D_{k} describes the line AM. Indeed, we have A = (cos(\frac{π}{2^{k}})^{−1}, 0) and M = (cos(\frac{π}{2^{k}}), sin(\frac{π}{2^{k}})) and it is straightforward that both of these points satisfy the last equality in the definition of D_{k}.

Recall now that the application

R_{θ} : ℝ^{2} ↦ ℝ^{2} : (\begin{matrix}x \\ y \end{matrix}) ↦ R_{θ} (\begin{matrix}x \\ y \end{matrix}) = (\begin{matrix}x cos θ + y sin θ \\ −x sin θ + y cos θ \end{matrix})

1Strictly speaking, this set is not a full-dimensional polyhedron in ℝ^{2k+2} because of the additional linear constraints but this has no incidence on our purpose.

is a clockwise rotation around the origin with angle θ. Calling P_{i} the point of ℝ^{2} whose coordinates are (α_{i}, β_{i}) and P̂_{i} = (α̂_{i}, β̂_{i}) the image of P_{i} by rotation R_{π/2^{i}}, we have that the first three constraints in the definition of D_{k} are equivalent to α_{i+1} = α̂_{i}, β_{i+1} ≥ β̂_{i} and −β_{i+1} ≤ β̂_{i}. These last two inequalities rewritten as −β_{i+1} ≤ β̂_{i} ≤ β_{i+1} immediately imply that β_{i+1} has to be nonnegative. Under this assumption, we call P̄_{i} the points whose coordinates are (α_{i}, −β_{i}) and find that these three constraints are equivalent to saying that P̂_{i} ∈ [P_{i+1} P̄_{i+1}]. In other words, the point P̂_{i} has to belong to a vertical segment [P_{i+1} P̄_{i+1}] such that P_{i+1} has its second coordinate nonnegative. Since P̂_{i} is the image of P_{i} by a rotation of angle π/2^{i}, saying that P̂_{i} belongs to some set is equivalent to saying that P_{i} belongs to the image of this set by the inverse rotation. In our case, this means in fine that P_{i} has to belong to the image by a rotation of angle −π/2^{i} of a segment [P_{i+1} P̄_{i+1}] such that P_{i+1} has its second coordinate nonnegative.

We can now specialize this result to i = k − 1. Recall that P_{k} is known to belong to the line AB. According to the above discussion, we have first to restrict this set to its points with a nonnegative β_{k}, which gives the half line [AB. Taking the union of all segments [P_{k} P̄_{k}] for all possible P_{k}’s gives the region bounded by half lines [AB and [AB′. Taking finally the image of this set by a rotation of angle −π/2^{k−1}, we find that P_{k−1} has to belong to the region bounded by half lines [BA and [BC.

We can now iterate this procedure and describe the set of points P_{k−2}, P_{k−3}, etc. Indeed, using exactly the same reasoning, we find that the set of points P_{i−1} can be deduced from the set of points P_{i} with a three-step procedure:

a. Restrict the set of points P_{i} to those with a nonnegative β_{i} coordinate.

b. Consider the union of segments [P_{i} P̄_{i}] where P_{i} belongs to the above restricted set, i.e. add for each point (α_{i}, β_{i}) the set of points (α_{i}, x) for all x ranging from −β_{i} to β_{i}.

c. Rotate this union counterclockwise around the origin with an angle equal to π/2^{i} to find the set of points P_{i−1}.

In the case of our example with k = 3, we have already shown that the set {P_{k}} = {P_{3}} = [AB and {P_{k−1}} = {P_{2}} is the region bounded by [BA ∪ [BC. Going on with the procedure described above, we readily find that {P_{k−2}} = {P_{1}} is the region bounded by the polygonal line [ABCDE] while {P_{k−3}} = {P_{0}} is the complete octagon ABCDED′C′B′A, which is the expected result (see Figure 2 for the corresponding pictures). It is not difficult to see that in the general case {P_{k−i}} is a set bounded by 2^{i} consecutive sides of P_{2^{k}}, which means we always end up with {P_{0}} equal to the whole regular 2^{k}-sided polyhedron P_{2^{k}}. This completes the proof since the set of points P_{0} is the projection of D_{k} on the subspace of the two variables (α_{0}, β_{0}).

This theorem allows us to derive quite easily a polyhedral approximation for B_{2}(1).

Corollary 2.1. The projection of D_{k} on the subspace of its two variables (α_{0}, β_{0}) is a polyhedral approximation of B_{2}(1) with accuracy ϵ = cos(\frac{π}{2^{k}})^{−1} − 1.

Proof. Straightforward application of Theorems 2.1 and 2.2.

This approximation is much better than the previous one: we have here that ϵ ≈ \frac{π^{2}}{2^{2k+1}} , which means that dividing the accuracy by four can be achieved by increasing k by 1, which corresponds to adding 2 variables, 1 equality and 2 inequality constraints (compare to the

Figure 2: The sets of points P_{3}, P_{2}, P_{1} and P_{0} when k = 3.

previous situation which needed to double the number of inequalities to reach to same goal). For example, an accuracy of ϵ = 10^{−4} can be obtained with k = 8, i.e. with 16 inequalities, 9 equalities and 18 variables (as opposed to 223 inequalities with the previous approach).

We are now in position to convert this polyhedral approximation of B_{2}(1) into an approximation of 𝕃^{2}. We define the set L_{k} ∈ ℝ^{2k+3} as

\{(r, α_{0}, …, α_{k}, β_{0}, …, β_{k}) ∈ ℝ^{2k+3} | \{ \begin{matrix}α_{i+1} & = & α_{i} cos \frac{π}{2^{i}} & + & β_{i} sin \frac{π}{2^{i}} \\ β_{i+1} & ≥ & β_{i} cos \frac{π}{2^{i}} & − & α_{i} sin \frac{π}{2^{i}} \\ −β_{i+1} & ≤ & β_{i} cos \frac{π}{2^{i}} & − & α_{i} sin \frac{π}{2^{i}} \\ r & = & α_{k} cos \frac{π}{2^{k}} & + & β_{k} sin \frac{π}{2^{k}} \end{matrix} ∀ 0 ≤ i < k\} .

Note the close resemblance between this set and D_{k}, the only difference being the introduction of an additional variable r in the last equality constraint. This set L_{k} is our final polyhedral approximation of 𝕃^{2}. Obviously, before we give proof of this fact, we need a measure of the quality of an approximation in the case of a second-order cone. This is the purpose of the next definition.

Definition 2.1. A set S ⊆ ℝ^{n+1} is said to be an ϵ-approximation of the second-order cone 𝕃^{n} if and only if we have

𝕃^{n} ⊆ S ⊆ 𝕃_{ϵ}^{n} = {(r, x) ∈ ℝ × ℝ^{n} | ∥x∥ ≤ (1 + ϵ)r}

where 𝕃_{ϵ}^{n} is an ϵ-relaxed second-order cone.

This definition extends our definition of ϵ-approximation for the unit disc B_{2}(1). The next theorem demonstrates how Corollary 2.1 on the accuracy of the polyhedral approximation D_{k} for B_{2}(1) can be converted into a result on the accuracy of L_{k} for 𝕃^{2}.

Theorem 2.3. The projection of L_{k} on the subspace of its three variables (r, α_{0}, β_{0}) is a polyhedral approximation of 𝕃^{2} with accuracy ϵ = cos(\frac{π}{2^{k}})^{−1} − 1.

Proof. Assuming r > 0 for the moment, we first establish a link between D_{k} and L_{k}. It is indeed straightforward to check using the corresponding definitions that the following equivalence holds

(r, α_{0}, …, α_{k}, β_{0}, …, β_{k}) ∈ L_{k} ⇔ (\frac{α_{0}}{r} , …, \frac{α_{k}}{r} , \frac{β_{0}}{r} , …, \frac{β_{k}}{r}) ∈ D_{k} ,  (2.3)

since

\{ \begin{matrix}α_{i+1} & = & α_{i} cos \frac{π}{2^{i}} & + & β_{i} sin \frac{π}{2^{i}} \\ β_{i+1} & ≥ & β_{i} cos \frac{π}{2^{i}} & − & α_{i} sin \frac{π}{2^{i}} \\ −β_{i+1} & ≤ & β_{i} cos \frac{π}{2^{i}} & − & α_{i} sin \frac{π}{2^{i}} \\ r & = & α_{k} cos \frac{π}{2^{k}} & + & β_{k} sin \frac{π}{2^{k}} \end{matrix} ⇔ \{ \begin{matrix}\frac{α_{i+1}}{r} & = & \frac{α_{i}}{r} cos \frac{π}{2^{i}} & + & \frac{β_{i}}{r} sin \frac{π}{2^{i}} \\ \frac{β_{i+1}}{r} & ≥ & \frac{β_{i}}{r} cos \frac{π}{2^{i}} & − & \frac{α_{i}}{r} sin \frac{π}{2^{i}} \\ −\frac{β_{i+1}}{r} & ≤ & \frac{β_{i}}{r} cos \frac{π}{2^{i}} & − & \frac{α_{i}}{r} sin \frac{π}{2^{i}} \\ 1 & = & \frac{α_{k}}{r} cos \frac{π}{2^{k}} & + & \frac{β_{k}}{r} sin \frac{π}{2^{k}} \end{matrix}

which means that L_{k} is nothing more than the homogenized polyhedral cone corresponding to D_{k}.

Let us now suppose (r, x_{1}, x_{2}) ∈ 𝕃^{2}. Equivalence (2.2) implies (\frac{x_{1}}{r} , \frac{x_{2}}{r}) ∈ B_{2}(1), which in turn implies by Corollary 2.1 that there exists a vector (α, β) ∈ ℝ^{2k} such that (\frac{x_{1}}{r} , α, \frac{x_{2}}{r} , β) belongs to D_{k}. Using the link (2.3), this last inclusion is equivalent to (r, x_{1}, rα, x_{2}, rβ) ∈ L_{k}, which means that (r, α_{0}, β_{0}) belongs to the projection of L_{k} on the subspace (r, α_{0}, β_{0}). We have thus shown that this projection is a relaxation of 𝕃^{2}, the first condition for it to be an ϵ-approximation of 𝕃^{2}.

Supposing now (r, x_{1}, x_{2}) belongs to the projection of L_{k}, there exists a vector (α, β) ∈ ℝ^{2k} such that (r, x_{1}, α, x_{2}, β) ∈ L_{k}. The equivalence (2.3) implies then that (\frac{x_{1}}{r} , \frac{α}{r}, \frac{x_{2}}{r} , \frac{β}{r}) ∈ D_{k}, which means that (\frac{x_{1}}{r} , \frac{x_{2}}{r}) belongs to the projection of D_{k} on its subspace (α_{0}, β_{0}). Using now Corollary 2.1, which states that this projection is an ϵ-approximation of B_{2}(1) with ϵ = cos(\frac{π}{2^{k}})^{−1} −1, we can write that ∥(\frac{x_{1}}{r} , \frac{x_{2}}{r})∥ ≤ 1+ϵ, which can be rewritten as ∥(x_{1}, x_{2})∥ ≤ (1+ϵ)r, which is the exactly the second condition for this projection to be an ϵ-approximation of 𝕃^{2}.

The last task we have to accomplish is to check what happens in the case where r ≤ 0. Suppose (r, x_{1}, α, x_{2}, β) ∈ L_{k}. Looking at the definition of L_{k}, and using the same reasoning as in the proof of Theorem 2.2, it is straightforward to show that the variables α_{0} and β_{0} can only be equal to 0 when r = 0, and that they cannot satisfy the constraints if r < 0 (i.e. in the first case the set {P_{0}} is equal to {(0, 0)} while in the second case {P_{0}} = ∅). Since this is also the situation of the second-order cone 𝕃^{2}, our approximation is exact when r ≤ 0, and we can conclude that the projection of L_{k} on the subspace of its three variables (r, α_{0}, β_{0}) is an ϵ-approximation of the three dimensional second-order cone 𝕃^{2} with ϵ = cos(\frac{π}{2^{k}})^{−1} −1.

2.5 Reducing the approximation

Our polyhedral approximation L_{k} features 2k + 3 variables, 2k linear inequalities and k + 1 linear equalities. It is possible to reduce these numbers by pivoting out a certain number of variables. Namely, using the set of constraints α_{i+1} = α_{i} cos \frac{π}{2^{i}} + β_{i} sin \frac{π}{2^{i}} for 0 ≤ i < k, we can replace α_{k} by a linear combination of α_{k−1} and β_{k−1}, then replace α_{k−1} by a linear combination of α_{k−2} and β_{k−2}, and so on until all variables α_{i} have been replaced except α_{0} (which cannot and should not be pivoted out since it belongs to the projected approximation). The last equality r = α_{k} cos \frac{π}{2^{k}} + β_{k} sin \frac{π}{2^{k}} can also be used to pivot out β_{k}.

The resulting polyhedron has then k+2 variables (r, α_{0}, β_{0}, …, β_{k−1}), 2k linear inequalities and no linear equality. However, it should be noted that the constraint matrix describing the reduced polyhedron is denser than in the original approximation, i.e. it contains many more nonzero elements, as depicted on Figure 3 in the case k = 15 (which also mentions the number of nonzero elements in each case).

Figure 3: Constraint matrices for L_{15} and its reduced variant.

This denser constraint matrix has of course a negative impact on the efficiency of the algorithm used to solve the approximation problems, so that computational experiments are needed to decide whether this is enough to counterbalance the advantage of a reduced number of equalities and variables Indeed, preliminary testing on a few problems representative of the ones we are going to consider in Section 3 led us to the conclusion that pivoting out the variables is beneficial, leading roughly to a 20% reduction of computing times.

Another interesting remark can be made when we have to approximate a second-order cone whose components are restricted to be nonnegative. Namely, if we know beforehand that x_{1} and x_{2} cannot be negative, the polyhedral approximation of 𝕃^{2} can be reduced. Indeed, looking back at the proof of Theorem 2.2, we see that the set of points P_{2} is bounded by 2^{k−2} consecutive sides of the regular 2^{k}-sided polyhedron (see for example the set {P_{2}} depicted in Figure 2). Combining this with the restriction that α_{2} and β_{2} are nonnegative, we have that the set {P_{2}} is exactly equal to the restriction of P_{2^{k}} to the positive orthant, and is thus a valid ϵ-approximation of 𝕃^{2} on this positive orthant with ϵ = cos(\frac{π}{2^{k}})^{−1} − 1.

This observation leads to the formulation of a reduced polyhedral approximation L′_{k} defined by

\{(r, α_{2}, …, α_{k}, β_{2}, …, β_{k}) ∈ ℝ^{2k−1} | \{ \begin{matrix}α_{i+1} & = & α_{i} cos \frac{π}{2^{i}} & + & β_{i} sin \frac{π}{2^{i}} \\ β_{i+1} & ≥ & β_{i} cos \frac{π}{2^{i}} & − & α_{i} sin \frac{π}{2^{i}} \\ −β_{i+1} & ≤ & β_{i} cos \frac{π}{2^{i}} & − & α_{i} sin \frac{π}{2^{i}} \\ r & = & α_{k} cos \frac{π}{2^{k}} & + & β_{k} sin \frac{π}{2^{k}} \end{matrix} ∀ 2 ≤ i < k\} ,

whose projection of the subspace of (r, α_{2}, β_{2}) approximates the nonnegative part 𝕃^{2}. This approximation features 2k − 1 variables, 2k − 4 linear inequalities and k − 1 linear equalities and can be reduced to k variables, 2k − 4 linear inequalities and no linear equality if we perform the pivoting described above.

At this stage, we would like to compare our approximation with the one presented in [BTN98]. Both of these feature the same accuracy ϵ = cos(\frac{π}{2^{k}})^{−1} − 1 (with parameter ν in [BTN98] equal to k − 1 in our setting). However, Ben-Tal and Nemirovski do not make explicit that the projection of their polyhedral approximation is equal to the regular 2^{k}-sided polyhedron in ℝ^{2} and only prove the corresponding accuracy result.

Table 1 compares the sizes of the polyhedral approximations in three cases: the original approximation, the reduced approximation where variables α_{i} are pivoted out and the nonnegative approximation L′_{k} (also with variables α_{i} pivoted out).

Table 1: Comparison of our approximation L_{k} with [BTN98].

Original  Reduced  Nonnegative
[BTN98] L_{k} [BTN98] L_{k} [BTN98] L′_{k}
Variables  2k + 3 2k + 3 k + 4 k + 2 k + 2  k
Inequalities 2k + 4  2k  2k + 4  2k  2k  2k − 4
Equalities  k − 1  k + 1  0  0  0  0

Our version uses 4 less inequality constraints in all three cases. It also features 2 more equality constraints in the original approximation, which turns out to be an advantage since it allows us to pivot out more variables in the reduced versions. Both the reduced and the nonnegative versions of L_{k} use 2 less variables than their counterparts in the original article of Ben-Tal and Nemirovski.

2.6 An approximation of 𝕃^{n}

We are now going to use the decomposition presented in Section 2.2 and our polyhedral approximation L_{k} for 𝕃^{2} to build an approximation for 𝕃^{n}. Recall that expression (2.1) decomposed 𝕃^{n} into ⌊n/2⌋ three-dimensional second-order cones 𝕃^{2} and a single larger cone 𝕃^{⌈ \frac{n}{2} ⌉}. Applying this decomposition recursively, we can decompose 𝕃^{⌈ \frac{n}{2} ⌉} into ⌊⌈n/2⌉/2⌋ secondorder cones 𝕃^{2} with a remaining larger cone 𝕃^{⌈⌈n/2⌉/2⌉}, which can be again decomposed into ⌊⌈⌈n/2⌉/2⌉/2⌋ cones 𝕃^{2}, etc.

Calling q_{k} the number of three-dimensional second-order cones appearing in the decomposition at each stage of this procedure and r_{k} the corresponding size of the remaining cone, we have initially q_{0} = 0, r_{0} = n and

q_{k} = ⌊\frac{r_{k−1}}{2} ⌋ and r_{k} = ⌈\frac{r_{k−1}}{2} ⌉ ∀ k > 0 .

Obviously, r_{k} is strictly decreasing and we must eventually end up with r_{k} equal to 2. Indeed, it is easy to see that 2^{i−1} < r_{k−1} ≤ 2^{i} implies 2^{i−2} < r_{k} ≤ 2^{i−1} and a simple recursive argument shows then that if 2^{m−1} < n ≤ 2^{m} we have 2^{m−k−1} < r_{k} ≤ 2^{m−k} and thus that r_{m−1} = 2. At this stage, the remaining second-order cone is 𝕃^{2}, which we can add to the decomposition in the last stage with q_{m} = 1 to have r_{m} = 0. Our decomposition has thus in total m stages. We also showed in Section 2.2 that the total of 𝕃^{2} cones in the final decomposition is equal to ∑_{i=1}^{m} q_{i} = n − 1, and we also note for later use that 2^{m−k} < r_{k−1} ≤ 2^{m−k+1} implies 2^{m−k−1} ≤ q_{k} ≤ 2^{m−k}.

We ask ourselves now what happens if each of the second-order cones appearing in this decomposition is replaced by an ϵ-approximation. Namely, suppose each of the ⌊n/2⌋ secondorder cones 𝕃^{2} in expression (2.1) is replaced by an ϵ^{(i)}-approximation (0 ≤ i ≤ ⌊n/2⌋), while the remaining larger cone is replaced by an ϵ′-approximation. We end up with the set

\{(r, x, y) ∈ ℝ_{+}×ℝ^{n+⌊ \frac{n}{2} ⌋} | (y_{i}, x_{2i−1}, x_{2i}) ∈ 𝕃_{ϵ^{(i)}}^{2}, 1 ≤ i ≤ ⌊\frac{n}{2}⌋, \begin{cases}(r, y) ∈ 𝕃_{ϵ′}^{⌈n/2⌉} & (n even) \\ (r, y, x_{n}) ∈ 𝕃_{ϵ′}^{⌈n/2⌉} & (n odd) \end{cases}\} ,

whose constraints are equivalent to

x_{2i−1}^{2} + x_{2i}^{2} ≤ (1 + ϵ^{(i)})^{2}y_{i}^{2}, 1 ≤ i ≤ ⌊\frac{n}{2}⌋, \begin{cases}∑i=1^{\frac{n}{2}} y_{i}^{2} & ≤ (1 + ϵ′)^{2}r^{2} (n even) \\ ∑i=1^{\frac{n}{2}} y_{i}^{2} + x_{n}^{2} & ≤ (1 + ϵ′)^{2}r^{2} (n odd) \end{cases} \} .

Ideally, we would like this decomposition to be an ϵ-approximation of 𝕃^{n}. We already know that it is a relaxation of 𝕃^{n}, since each approximation of 𝕃^{2} is itself a relaxation. We have thus to concentrate on the second condition defining an ϵ-approximation, ∥x∥ ≤ (1 + ϵ)r. Writing

∑_{i=1}^{2⌊n/2⌋} x_{i}^{2} ≤ ∑_{i=1}^{⌊n/2⌋}(1 + ϵ^{(i)})^{2}y_{i}^{2} ,

we would like to bound the quantity on the right hand-side. Unfortunately, we only know a bound on the sum of y_{i}^{2}’s, which forces us to write

∑_{i=1}^{2⌊n/2⌋} x_{i}^{2} ≤ (1 + max_{i} ϵ^{(i)})^{2} ∑_{i=1}^{⌊n/2⌋} y_{i}^{2} ⇒ ∑_{i=1}^{n} x_{i}^{2} ≤ (1 + max_{i} ϵ^{(i)})^{2}(1 + ϵ′)^{2}r^{2} .

This shows that our decomposition is an approximation of 𝕃^{n} with accuracy ϵ = (1 + max_{i} ϵ_{i})(1 + ϵ′). This immediately implies that there is no point in approximating with different accuracies the ⌊n/2⌋ small second-order cones 𝕃^{2} appearing in the decomposition, since only the largest of these accuracies has an influence on the resulting approximation for 𝕃^{n}. Applying now our decomposition recursively to the remaining cone, and choosing at each stage k a unique accuracy ϵ_{k} for all the 𝕃^{2} cones, we find that 1 + ϵ = ∏_{k=1}^{m}(1 + ϵ_{k}) ,, i.e. that the final accuracy of our polyhedral approximation is the product of the accuracies chosen at each stage of the decomposition (note that, unlike the situation for a single stage, there is no reason here to choose all ϵ_{k} accuracies to be equal to each other).

2.7 Optimizing the approximation

The previous section has shown how to build a polyhedral approximation of 𝕃^{n} and how its quality depends on the accuracy of the approximations used at each stage of the decomposition. Our goal is here to optimize these quantities, i.e. given a target accuracy ϵ for 𝕃^{n}, find the values of ϵ_{k} (1 ≤ k ≤ m) that lead to the smallest polyhedral approximation, i.e. the one with the smallest number of variables and constraints.

Let us suppose we use at stage k the approximation L_{u_{k}} with u_{k} + 2 variables and 2u_{k} linear inequalities (i.e. with variables α pivoted out of the formulation), which has an accuracy ϵ_{k} = cos(\frac{π}{2^{u_{k}}})^{−1} − 1. Recalling notation q_{k} for the number of cones 𝕃^{2} introduced at stage k of the decomposition, the final polyhedral approximation has thus an accuracy equal to ∏_{k=1}^{m} cos(\frac{π}{2^{u_{k}}})^{−1} with 2 ∑_{k=1}^{m} q_{k}u_{k} inequalities and n + ∑_{k=1}^{m} q_{k}u_{k} variables. Indeed, we have n original x_{i} variables and u_{k} additional variables for each of the q_{k} approximations at stage k, since the first two variables in these approximations are coming from the previous stage. We observe that the main quantity to be minimized is ∑_{k=1}^{m} q_{k}u_{k} for both the number of variables and inequalities, which leads to the following optimization problem:

σ_{n,ϵ} = min_{u∈ℕ^{m}} ∑_{k=1}^{m} q_{k}u_{k} s.t. ∏_{k=1}^{m} cos(\frac{π}{2^{u_{k}}})^{−1} ≤ 1 + ϵ .  (2.4)

A possible choice for variables u_{k} is to take them all equal. Plugging this unique value into the accuracy constraint, we readily find that u_{k} has to be equal to

u_{k} = ⌈log_{2}(π/ arccos(1 + ϵ)^{−1/m})⌉

and that when the dimension of the cone 𝕃^{n} (and thus m) tends to +∞, we have u_{k} = O (log \frac{m}{ϵ}) and σ_{n,ϵ} = O (n log \frac{m}{ϵ}).

This obviously does not lead to an optimal solution of (2.4). Indeed, since the number of approximations is decreasing as we move from one stage of the decomposition to the next, it is intuitively clear that trading a lower accuracy for first stages against a higher accuracy for the last stages will be beneficial, since the lowering of the number of variables and inequalities in the first stages will affect many more constraints than the increase of size for the last stages. This implies that the components u_{k} of any optimal solution of (2.4) will have to be in increasing order.

Finding a closed form optimal solution of (2.4) does not appear to be possible, but we can find a good suboptimal solution using some approximations. We first introduce variables v_{k} such that v_{k} = 4^{−u_{k}} ⇔ u_{k} = − log_{4} v_{k} and rewrite problem (2.4) as

σ_{n,ϵ} = min_{v∈ℝ^{m}} − log_{4} ∏_{k=1}^{m} v_{k}^{q_{k}} s.t. ∑_{k=1}^{m} log(cos(π\sqrt{v_{k}})^{−1}) ≤ log(1 + ϵ) .

Since u_{k} ≥ 2, we have π\sqrt{v_{k}} ≤ \frac{π}{4} and we can use the easily proven^{2} inequality

log(cos(x)^{−1}) ≤ (\frac{3x}{4})^{2}

valid for all 0 ≤ x ≤ \frac{π}{4} to write

σ_{n,ϵ} = − max_{v∈ℝ^{m}} log_{4} ∏_{k=1}^{m} v_{k}^{q_{k}} s.t. ∑_{k=1}^{m} v_{k} ≤ \frac{16}{9}π^{−2} log(1 + ϵ) = K(ϵ) ,

2This inequality can be easily checked by plotting the graphs of its two sides on the interval [0 \frac{π}{4} ].

which is thus a restriction of our original problem. It amounts to maximizing a product of variables whose sum is bounded, a problem whose optimality conditions are well-known. In our case, they can written as

\frac{v_{1}}{q_{1}} = \frac{v_{2}}{q_{2}} = … \frac{v_{m}}{q_{m}} = \frac{∑_{k=1}^{m} v_{k}}{∑_{k=1}^{m} q_{k}} = \frac{K(ϵ)}{n − 1} ⇒ v_{k} = \frac{q_{k}}{n − 1}K(ϵ) ⇔ u_{k} = log_{4} \frac{n − 1}{q_{k}} − log_{4} K(ϵ) .

However u_{k} must be integer, so that we have to degrade this solution further and round it towards a larger integer. Using that fact that n − 1 ≤ 2^{m} and q_{k} ≥ 2^{m−k−1}, we have

\frac{n − 1}{q_{k}} ≤ \frac{2^{m}}{2^{m−k−1}} = 2^{k+1} ⇒ log_{4} \frac{n − 1}{q_{k}} ≤ log_{4} 2^{k+1} = \frac{k + 1}{2} ,

so that we can take u_{k} = ⌈\frac{k+1}{2}⌉ − ⌊log_{4} K(ϵ)⌋ as our suboptimal integer solution for (2.4). Let us plug now these values into the objective function σ_{n,ϵ}: we find

∑_{k=1}^{m} q_{k}u_{k} = ∑_{k=1}^{m} q_{k}⌈\frac{k + 1}{2}⌉ − ∑_{k=1}^{m} q_{k}⌊log_{4} K(ϵ)⌋

≤ ∑_{k=1}^{m} q_{k}(\frac{k}{2} + 1) − (n − 1)⌊log_{4} K(ϵ)⌋ (using ∑_{k=1}^{m} q_{k} = n − 1)

≤ ∑_{k=1}^{m} 2^{m−k}(\frac{k}{2} + 1) − (n − 1)⌊log_{4} K(ϵ)⌋ (using q_{k} ≤ 2^{m−k})
≤ 2^{m+1} − \frac{m}{2} − 2 − (n − 1)⌊log_{4} K(ϵ)⌋
≤ 4(n − 1) − (n − 1)⌊log_{4} K(ϵ)⌋ (using 2^{m−1} ≤ n − 1)
≤ (n − 1)⌈4 − log_{4} K(ϵ)⌉ = (n − 1)⌈4 − log_{4} \frac{16}{9} + log_{4} π^{2} − log_{4} log(1 + ϵ)⌉
≤ (n − 1)⌈5.3 − log_{4} log(1 + ϵ)⌉

(where we have used at the fourth line the fact that ∑_{k=1}^{m} 2^{m−k}(\frac{k}{2} +1) = 2^{m+1} − \frac{m}{2} −2, which is easily proved recursively). We can wrap this result into the following theorem:

Theorem 2.4. For every ϵ < \frac{1}{2}, there exists a polyhedron with no more than

2 + (n − 1)⌈5.3 − log_{4} log(1 + ϵ)⌉ = O(n log \frac{1}{ϵ}) variables

and

2 + 2(n − 1)⌈4.3 − log_{4} log(1 + ϵ)⌉ = O(n log \frac{1}{ϵ}) inequalities

whose projection on a certain subspace of n + 1 variables is an ϵ-approximation of the secondorder cone 𝕃^{n} ⊆ ℝ^{n+1}.

Proof. This a consequence of the previous derivation, which showed that choosing u_{k} = ⌈\frac{k+1}{2}⌉ − ⌊log_{4} K(ϵ)⌋ lead to an ϵ-approximation of 𝕃^{n} with n + σ_{n,ϵ} variables and 2σ_{n,ϵ} linear inequalities, with σ_{n,ϵ} = (n − 1)⌈5.3 − log_{4} log(1 + ϵ)⌉. However, the size of this polyhedron can be further reduced using L′_{k}, the polyhedral approximation of the nonnegative part of 𝕃^{2}. Indeed, looking at the decomposition (2.1), we see that all the y variables used in the second stage of the decomposition are guaranteed to be nonnegative, since we have in our approximation (y_{i}, x_{2i−1}, x_{2i}) ∈ L_{u_{1}} which implies y_{i} ≥ 0. This means that we can use for the second stage and the following our reduced approximation L′_{u_{k}}, known to be valid when its first two variables are restricted to the nonnegative orthant 𝕃^{2}, which uses 2 less variables and 4 less inequalities per cone. Since there are \frac{n}{2} cones in the first stage of the decomposition and n − 1 cone in total, we can use \frac{n}{2} − 1 reduced approximations L′, which give us a total saving of n − 2 variables and 2n − 4 constraints^{3}. Combining this with the value of σ_{n,ϵ}, we find that our approximation has 2+(n−1)⌈5.3−log_{4} log(1+ϵ)⌉ variables and 2+⌈4.3−log_{4} log(1+ϵ)⌉ inequalities.

We also have to prove the asymptotic behaviour of σ_{n,ϵ} when n tends to infinity. Indeed, we have log(1 + ϵ) ≥ \frac{ϵ}{2} when ϵ < \frac{1}{2}, which implies − log_{4} log(1 + ϵ) ≤ log_{4} \frac{1}{ϵ}. This leads to ⌈5.3 − log_{4} log(1 + ϵ)⌉ = O(log \frac{1}{ϵ}), which is enough to prove the theorem.

This result is better than the one we previously obtained choosing all u_{k}’s equal to each other: indeed, we had in that case σ_{n,ϵ} = O(n log \frac{m}{ϵ}) = O(n log \frac{log n}{ϵ}) while we have here O(n log \frac{1}{ϵ}). For a fixed accuracy ϵ, our first choice translates into σ_{n,ϵ} = O (n log log n) when n tends to infinity while our optimized choice of u_{k}’s leads to σ_{n,ϵ} = O (n), which is better. We note however that if we fix n and let ϵ tends to 0, the asymptotic behaviour is the same in both cases , namely we have σ_{n,ϵ} = O(log \frac{1}{ϵ}).

Ben-Tal and Nemirovski achieve essentially the same result in [BTN98], albeit in the special case when n is a power of two. Our proof has the additional advantage of providing a closed form for the parameters u_{k} as well as for the total size of the polyhedral approximation, for all values of n. They also prove that the number of inequalities of a ϵ-approximation of 𝕃^{n} must be greater that O(n log \frac{1}{ϵ}), i.e. that the order of the result of Theorem 2.4 is not improvable.

2.8 An approximation of second-order cones optimization

The previous sections have proven the existence of a polyhedral approximation of the secondorder cone with a moderate size (growing linearly with the dimension of the cone and the logarithm of the accuracy). However, we have to point out that these polyhedrons are not strictly speaking approximations of the second-order cone: more precisely, it is their projection on a certain subspace than is an ϵ-approximation of 𝕃^{n}.

This does not pose any problem when trying to approximate a second-order cone optimization problem with linear optimization. Let us suppose we want to approximate problem (1.3), which we recall here for convenience,

inf c^{T} x s.t. Ax = b and x^{k} ∈ 𝕃^{n_{k}} ∀k = 1, 2, …, r  (1.3)

with ϵ-approximations of the second order cones 𝕃^{n_{k}}. Theorem 2.4 implies the existence of a polyhedron

Q_{k} = \{(x^{k}, y^{k}) ∈ ℝ^{n_{k}+1} × ℝ^{O(n_{k} log \frac{1}{ϵ})} | A_{k}(x^{k}, y^{k})^{T} ≥ 0\}

with

A_{k} ∈ ℝ^{O(n_{k} log \frac{1}{ϵ})×O(n_{k} log \frac{1}{ϵ})} ,

3This reasoning was made for an even n. In the case of an odd n, we have one cone in the decomposition for which only the first variable is known to be nonnegative. It is possible to show that there exists a polyhedral approximation adapted to this situation that uses 1 less variable and 2 less inequalities than the regular approximation, which allows us to write exactly the same results as for an even n.

whose projection on the subspace (r, x) is an ϵ-approximation of 𝕃^{n_{k}}, which allows us to write the following linear optimization problem^{4}

min c^{T} x s.t. Ax = b and A_{k}(x^{k}, y^{k})^{T} ≥ 0 ∀k = 1, 2, …, r .  (2.5)

We note that the fact that our approximations are projections is handled in a seamless way by this formulation: the only difference with the use of a direct approximation of the cones 𝕃^{n_{k}} is the addition of the auxiliary variables y^{k} to the formulation. This problem features ∑_{k=1}^{r} O(n_{k} log \frac{1}{ϵ}) = O(n log \frac{1}{ϵ}) variables, m equality constraints and ∑_{k=1}^{r} O(n_{k} log \frac{1}{ϵ}) = m + O(n log \frac{1}{ϵ}) homogeneous inequality constraints. We also point out as a minor drawback of this formulation the fact that it involves irrational coefficients, namely the quantities sin \frac{π}{2^{i}} and cos \frac{π}{2^{i}} occurring in the definition of L_{k}. However, it rational coefficients are really needed (for example if one wants to work with a complexity model based on exact arithmetic), it is possible to replace those quantities with rational approximations while keeping an essentially equivalent accuracy for the resulting polyhedral approximation, i.e. featuring the same asymptotic behaviour.

To conclude this section, we are going to compare the algorithmic complexity of solving problem (1.3) either directly or using our polyhedral approximation. The best complexity obtained so far^{5} for solving a linear program with v variables up to accuracy ϵ is O(v^{3.5} log(\frac{1}{ϵ})) arithmetic operations (using for example a short-step path-following method, see e.g. [NN94]). In our case, assuming we solve the approximate problem (2.5) up to the same accuracy that the one used to approximate the second-order cones, this leads to a complexity equal to O(n^{3.5} log(\frac{1}{ϵ})^{4.5}).

On the other hand, solving problem (1.3) can be done using O(\sqrt{r} n^{3} log(\frac{1}{ϵ})) arithmetic operations, using for example a potential reduction approach, see e.g. [LVBL98]. If r = O (1), i.e. if the number of cones used in the formulation is bounded, the second complexity is better, both if n → +∞ or ϵ → 0. However, if r = O (n), which means that the dimension of the cones used in the formulation is bounded, both complexity become equivalent from the point of view of the dimension n, but the second one is still better when letting the accuracy tend to 0. We conclude that the direct solving of (1.3) as a second-order cone problem is superior from the point of view of algorithmic complexity. The purpose of the second part of this article will be to test whether this claim is also valid for computational experiments.

2.9 Accuracy of the approximation

The linearizing scheme for second-order cone optimization presented in the previous section is based on a polyhedral approximation whose accuracy is guaranteed in the sense of Definition 2.1. It is important to realize that this bound on the accuracy of the approximation does not imply a bound on the accuracy of the solutions (or the objective value) of the approximated problem.

Indeed, let us consider the following set:

\{(r, x_{1}, x_{2}) ∈ ℝ^{3} | r − x_{2} = \frac{1}{2} and (r, x_{1}, x_{2}) ∈ 𝕃^{2}\} .

4We could replace the inf of problem (1.3) by a min since it is well-known that linear optimization problems always attain their optimal objectives.

5Using standard linear algebra and without partial updating.

This set can be seen as the feasible region of a second-order cone problem. Using the fact that

x_{1}^{2} + x_{2}^{2} ≤ r^{2} ⇔ x_{1}^{2} + x_{2}^{2} ≤ (x_{2} + \frac{1}{2})^{2} ⇔ x_{1}^{2} − \frac{1}{4} ≤ x_{2} ,

we find that the projection of this set on the subspace (x_{1}, x_{2}) is the epigraph of the parabola x ↦ x^{2} − \frac{1}{4}. Let us now replace 𝕃^{2} by the polyhedral approximation L_{k}. Since the resulting set will be polyhedral, its projection on the subspace (x_{1}, x_{2}) will also be polyhedral, and we can deduce without difficulties that it is the epigraph of a piecewise linear function, as shown by Figure 4 (depicting the cases k = 1, 2, 3 and 4).

Figure 4: Linear approximation of a parabola using L_{k} for k = 1, 2, 3, 4.

Because a polyhedron has a finite number of vertices, this piecewise linear function must have a finite number of segments. Considering the rightmost piece, i.e. whose x_{1} values span an interval of the type [α + ∞[, it is obvious that it cannot approximate the parabola with a guaranteed accuracy. Indeed, the difference between the approximation and the parabola grows quadratically on this segment, which means that even the ratio of variable x_{2} between the parabola and its linear approximation is not bounded.

Let us now consider the following parameterized family of second-order cone optimization problems

min x_{2} s.t. x_{1} = λ, r − x_{2} = \frac{1}{2} and (r, x_{1}, x_{2}) ∈ 𝕃^{2} ,  (PB_{λ})

which is using the same feasible set as above with an additional constraint fixing the variable x_{1} to λ. Denoting the optimal objective value of (PB_{λ}) by p^{∗}(λ), we have in light of the previous discussion that p^{∗}(λ) = λ^{2} − \frac{1}{4}. However, we also showed that the optimal objective value p_{k}^{∗}(λ) of the approximated problem

min x_{2} s.t. x_{1} = λ, r − x_{2} = \frac{1}{2} and (r, x_{1}, x_{2}, y) ∈ L_{k}

must be a piecewise linear function of λ with a finite number of segments. Indeed, simple computations^{6} show that the endpoints of these segments occur for

λ = \frac{sin(iθ)}{2 cos \frac{θ}{2} − 2 cos(iθ)} for i = 1, 2, …, 2^{k} − 1 with θ = \frac{π}{2^{k−1}} ,

which shows that p^{∗}(λ) is linear as soon as λ ≥ _{\frac{sin θ}{2 cos \frac{θ}{2} −2 cos θ}}.

The discrepancy between the real optimum p^{∗}(λ) and the approximated optimum p_{k}^{∗}(λ) is thus unbounded when λ goes to infinity. Moreover, we have that the relative accuracy of p_{k}^{∗}(λ) tends to 1, the worst possible value, i.e.

\frac{p^{∗}(λ) − p_{k}^{∗}(λ)}{p^{∗}(λ)} → 1 .

Another interesting feature of this small example is that performing a complete parametric analysis for parameter λ ranging from −∞ to +∞ would lead to 2^{k} −1 different break points.

We conclude that we cannot give an a priori bound on the accuracy of the optimal objective value of the linear approximation of a second-order cone optimization problem (this remark is also valid for accuracy of the optimal solution itself, since we have in our example (PB_{λ}) that the optimal value of x_{2} is equal to p^{∗}).

3 Computational experiments

In this section, we present computational experiments with an implementation of the linearizing scheme for second-order cone optimization we have just described.

3.1 Implementation

The computer used to conduct those experiments is an Intel 500 MHz Pentium III with 128 megabytes of memory. We chose to use the MATLAB programming environment, developed by The MathWorks, for the following reasons:

⋄ MATLAB is a flexible and modular environment for technical computing, two very important characteristics when developing research code. Although MATLAB may be somehow slower than a pure C or FORTRAN approach, we think that this loss of performance is more than compensated by the ease of development (especially from the point of view of graphic capabilities and debugging). Moreover, the critical (i.e. time consuming) parts of the algorithms can be coded separately in C or FORTRAN and used in MATLAB via MEX files (this is the approach taken by the solvers we mention below), which allows a well designed MATLAB program to be nearly as efficient as an equivalent pure C or FORTRAN program.

6Simply observe that the extremal rays of L_{k} obey to the relation x_{2} = x_{1} tan iθ with i = 1, 2, …, 2^{k} and θ = π/2^{k−1}.

⋄ Efficient interior-point solvers are available on the MATLAB platform. Indeed, we used in our experiments

– The MOSEK optmization toolbox for MATLAB by EKA Consulting ApS [AA99, ART00], a full-featured optimization package including a simplex solver and primaldual interior-point solvers for linear optimization, convex linearly and quadratically constrained optimization, second-order cone optimization, linear least square problems, linear l_{1} and l_{∞}-norm optimization and geometric and entropy optimization. When compared with the standard optimization toolbox from MATLAB, MOSEK is particularly efficient on large-scale and sparse problems. MOSEK can be downloaded for research and evaluation purposes at http://www.mosek.com.

– SeDuMi by Jos Sturm [Stu99b], another primal-dual interior-point solver which is able to handle linear, second-order cone and semidefinite optimization problems. SeDuMi is designed to take into account sparsity and complex values, and has the advantage of dealing with the very important class of semidefinite optimization, but is a little more restrictive than MOSEK concerning the input format, since problems must be entered in the standard conic form (1.3). SeDuMi can be

    downloaded at http://www.unimaas.nl/~sturm/.

The main routines we implemented are the following:

⋄ PolySOC2(k) generates the polyhedron L_{k} with accuracy ϵ_{k} = cos(\frac{π}{2^{u_{k}}})^{−1} −1. Variables α_{i} are pivoted out, so that this routine returns a polyhedron with k + 2 variables and 2k inequalities. An optional parameter is available to use the reduced approximation L′_{k}, valid on the nonnegative restriction of 𝕃^{2}.

⋄ Steps(q, e) computes the optimal choice for the size of the cones at each stage of the decomposition of 𝕃^{n}. Indeed, q contains our vector q (i.e. the number of cones at each stage) and e is the target accuracy ϵ.

⋄ PolySOCN(n, e) generates a e-approximation of 𝕃^{n}. It uses the output of PolySOC2 and the optimal sizes for the cones computed by Steps.

⋄ PolySOCLP(p, e) linearizes the second-order cone optimization problem p, replacing each second-order cone constraint with a polyhedral e-approximation using PolySOCN and outputting a linear optimization problem.

The procedure Steps we implemented features some improvements when compared with the theory we presented in the previous section. Indeed, Theorem 2.4 shows that the choice u_{k} = ⌈\frac{k+1}{2}⌉ − ⌊log_{4} K(ϵ)⌋ leads to a polyhedron of size O(n log \frac{1}{ϵ}), but is not optimal for two reasons:

⋄ We approximated the formula giving the accuracy of the approximation to derive u_{k} (namely, we used log(cos(x)^{−1}) ≤ (\frac{3x}{4})^{2}).

⋄ The optimal solution for this approximated accuracy was not guaranteed to be integer and had to be rounded to the smallest greater integer.

However, one can easily improve this choice in practice as follows. Let us suppose theory predicts some optimal values v_{k} for the sizes of the cones at stage k, which have to be rounded to ⌈v_{k}⌉. Because of this rounding, the actual accuracy of the approximation will be much better than our target ϵ. Recalling now that this accuracy is equal to (1 + ϵ_{1})(1 + ϵ′), where ϵ_{1} is the accuracy of the cones in the first stage, and is thus equal to cos(\frac{π}{2^{⌈v_{1}⌉}})^{−1} − 1, and ϵ′ is the accuracy from the cone modelled by all the remaining stages 𝕃^{⌈n/2⌉}, we can compute an upper bound for ϵ′, according to

(1 + ϵ_{1})(1 + ϵ′) ≤ ϵ ⇔ ϵ′ ≤ \frac{ϵ}{1 + ϵ_{1}} − 1

which will be better (i.e. higher) than in the theoretical derivation since it takes into account the exact accuracy ϵ_{1} of the first stage, rounding included. We can now apply this procedure to the second stage, i.e. computing a theoretical value for ϵ_{2} and an upper bound on the accuracy of 𝕃^{⌈⌈n/2⌉/2⌉}, and so on, obtaining in the end a smaller polyhedral approximation, since the required accuracies for the cones at every stage (except the first one) are higher and hence need less constraints and variables.

Still, this improved rounding does not address the first reason why our u_{k}’s are not optimal, the fact that we do not optimize the actual formula for the accuracy. Since it seems impossible to deal with it in closed form, we implemented a dynamic programming approach to optimize it. This algorithm uses the theoretical suboptimal solution described above (including the improved rounding procedure) to provide bounds on the optimal solution and therefore reduce the computing time.

Figure 5 presents the graphs of the size (measured by σ_{n,ϵ}) of our approximations in two situations: fixed dimension (n = 50, 200, 800) with accuracy ranging from 10^{−1} to 10^{−8} and fixed accuracy (ϵ = 10^{−2}, 10^{−5}, 10^{−8}) with dimension ranging from 10 to 1000. The asymptotic behaviour of σ_{n,ϵ} is very clear on these graphs: we have a linear increase when n tends to +∞ for a fixed accuracy and a logarithmic increase when ϵ tends to 0 for a fixed dimension (since the first graph has a logarithmic scale of abscissas).

Figure 5: Size of the optimal approximation versus accuracy (left) and dimension (right).

Finally, in order to give an idea of the efficiency of our improved rounding procedure and dynamic programming resolution, we provide in Table 2 the value of σ_{n,ϵ} for different strategies, using a target accuracy ϵ = 10^{−8}.

Table 2: Different approaches to optimize the size of a 10^{−8}-approximation of 𝕃^{n}

n Theory All equal (rounded) Theory (rounded) Dynamic programming
10  171  139  141  139
100  1881  1584  1552  1537
1000 18981  15987  15688  15522
10000 189981  160140  157063  155392

The first column represents the theoretical value σ_{n,ϵ} = (n − 1)⌈5.3 − log_{4} log(1 + ϵ)⌉, the second column describes the choice of all u_{k}’s, equal to each other, albeit using the improved rounding procedure presented above, the third column reports the choice of the theoretical value for u_{k}, this time with the improved rounding procedure, and the last column gives the true optimal value via our dynamic programming approach. We observe that our iterative rounding procedure improves the theoretical value of σ_{n,ϵ} in a noticeable way, lowering it by approximately 15%. The differences between the last three columns are less important, the dynamic programming approach giving a few additional percents of decrease in the size of the approximation.

3.2 Truss-topology design

We first tested our linearizing scheme for second-order cone optimization on a series of truss topology design problems. A truss is a structure composed of elastic bars connecting a set of nodes, like a railroad bridge or the Eiffel tower. The task consists in determining the size (i.e. the cross sectional areas) of the bars that lead to the stiffest truss when submitted to a set of forces, subject to a total weight limit. The problem we want to solve here is a multi-load truss topology design problem, which means we are simultaneously considering a set of k loading scenarios. This problem can be formulated as follows (see [BTN94]):

min ∑_{i=1}^{n} σ_{i} s.t. ∥(q_{i1}, …, q_{ik})∥ ≤ σ_{i} ∀1 ≤ i ≤ n and B (q_{1j} ; q_{2j} ; ⋮ ; q_{nj}) = f_{j} ∀1 ≤ j ≤ k , (TTD)

where n is the number of bars, k is the number of loading, vector σ ∈ ℝ^{n} and matrix Q ∈ ℝ^{n×k} are the design variables, B ∈ ℝ^{m×n} is a matrix describing the physical configuration of the truss and f_{j} ∈ ℝ^{m}, 1 ≤ j ≤ k are vectors of forces describing the loadings scenarios.

It is easily cast as a second-order cone problem in the form (1.3), since the norm constraints can be modelled as (σ_{i}, q_{i1}, …, q_{ik}) ∈ 𝕃^{k} for all 1 ≤ i ≤ n. Indeed, we have that the variables x ∈ ℝ^{k(n+1)}, the objective c ∈ ℝ^{k(n+1)} and the equality constraints Ax = b with A ∈ ℝ^{km×k(n+1)} and b ∈ ℝ^{km} are given by

x = ( σ ; q_{11} ; q_{21} ; ⋮ ; q_{nk}) , c = (1^{n×1} ; 0 ; 0 ; ⋮ ; 0) , A = (\begin{matrix}0^{m×n} & B \\ 0^{m×n} & & B \\ ⋮ & & & ... \\ 0^{m×n} & & & & B \end{matrix}) and b = (f_{1} ; f_{2} ; ⋮ ; f_{k})

to give the following second-order cone optimization problem equivalent to (TTD):

min (1^{n×1} ; 0 ; 0 ; ⋮ ; 0)^{T} ( σ ; q_{11} ; q_{21} ; ⋮ ; q_{nk}) s.t. \{ \begin{matrix}(\begin{matrix}0^{m×n} & B \\ 0^{m×n} & & B \\ ⋮ & & & ... \\ 0^{m×n} & & & & B \end{matrix}) ( σ ; q_{11} ; q_{21} ; ⋮ ; q_{nk}) = (f_{1} ; f_{2} ; ⋮ ; f_{k}) \\ (σ_{i}, q_{i1}, …, q_{ik}) ∈ 𝕃^{k} for all 1 ≤ i ≤ n \end{matrix} . (TP)

This allows us to write a dual problem in the form (1.4) in a straightforward manner:

max ∑_{j=1}^{k} f_{j}^{T} y_{j} s.t. \{ \begin{matrix}(\begin{matrix}0^{m×n} & 0^{m×n} & ⋯ & 0^{m×n} \\ B^{T} \\ & B^{T} \\ . & . & . \\ & & & B^{T} \end{matrix}) (y_{σ} ; y_{1} ; ⋮ ; y_{k}) + ( σ^{∗} ; q_{11}^{∗} ; q_{21}^{∗} ; ⋮ ; q_{nk}^{∗}) = (1^{n×1} ; 0 ; 0 ; ⋮ ; 0) \\ (σ_{i}^{∗}, q_{i1}^{∗}, …, q_{ik}^{∗}) ∈ 𝕃^{k} for all 1 ≤ i ≤ n \end{matrix} (TD)

where σ^{∗} ∈ ℝ^{n}, Q^{∗} ∈ ℝ^{n×k}, y_{σ} ∈ ℝ^{m} and y_{j} ∈ ℝ^{m} for all 1 ≤ j ≤ k are the dual variables. Variables σ^{∗} and Q^{∗} can be pivoted out of the formulation using the linear constraints

σ^{∗} = 1^{n×1} and q_{ij}^{∗} = −b_{i}^{T} y_{j} ∀1 ≤ i ≤ n, 1 ≤ j ≤ k

(where b_{i} ∈ ℝ^{m} is the i^{th} column of B), which gives then

max ∑_{j=1}^{k} f_{j}^{T} y_{j} s.t. (1, −b_{i}^{T} y_{1}, …, −b_{i}^{T} y_{k}) ∈ 𝕃^{k} for all 1 ≤ i ≤ n

and finally

max ∑_{j=1}^{k} f_{j}^{T} y_{j} s.t. ∑_{j=1}^{k}(b_{i}^{T} y_{j})^{2} ≤ 1 for all 1 ≤ i ≤ n ,  (TQC)

which is a convex quadratically constrained problem with a linear objective. We can thus solve a truss-topology design problem in at least three different manners: either solving the secondorder cone optimization problems (TP) or (TD) or solving the quadratically constrained problem (TQC).

The problems we used for our computational experiments we randomly created using a generator developed by A. Nemirovski. Given three integers p, q and k, it produced the matrix B and vectors f_{j} corresponding to k loading scenarios for a truss using a 2-dimensional p × q

Table 3: Dimensions of the truss-topology design problems.

Problem description  Formulation  Primal  Dual
2 × 2 grid with 2 loads 5 cones 𝕃^{2}  15 × 8  23 × 15
2 × 2 grid with 4 loads 5 cones 𝕃^{4}  25 × 16  41 × 25
2 × 2 grid with 8 loads 5 cones 𝕃^{8}  45 × 32  77 × 45
2 × 2 grid with 16 loads 5 cones 𝕃^{16}  85 × 64  149 × 85
2 × 2 grid with 32 loads 5 cones 𝕃^{32}  165 × 128  293 × 165
4 × 4 grid with 2 loads 114 cones 𝕃^{2}  342 × 48  390 × 342
4 × 4 grid with 4 loads 114 cones 𝕃^{4}  570 × 96  666 × 570
4 × 4 grid with 6 loads 114 cones 𝕃^{6}  798 × 144  942 × 798
6 × 6 grid with 2 loads 615 cones 𝕃^{2} 1845 × 120 1965 × 1845
6 × 6 grid with 4 loads 615 cones 𝕃^{4} 3075 × 240 3315 × 3075
8 × 8 grid with 2 loads 1988 cones 𝕃^{2} 5964 × 224 6188 × 5964

nodal grid, with n ≈ \frac{1}{2}p^{2}q^{2} and m ≈ 2pq. We tested 11 combinations of parameters p, q and k. The dimensions of the corresponding problems are reported in Table 3 (the last two columns report the number of variables × the number of constraints). We see that the last problems involve a fairly large number of small second-order cones.

Polyhedral approximations of these problems were computed for three different accuracies, namely ϵ = 10^{−2}, 10^{−5} and 10^{−8}. The dimensions of the resulting linear optimization problems are reported in Table 4. It is interesting to note that problems with accuracy 10^{−8} are only approximately three times larger than problems with accuracy 10^{−2} and 50% larger than problems with accuracy 10^{−5} (with several dozens of thousands of variables for the largest among them).

Before we present computing times, we have to mention a special feature of this class of problems. Contrary to the general assertion that is stated in Section 2.9, it is possible to give here an estimation of the quality of the optimum objective value of the approximated problem. Indeed, let us call t^{∗} the optimal objective value of problem (TTD) and t_{ϵ}^{∗} the optimal objective value of the approximated problem with accuracy ϵ. Since our approximation is a relaxation, we obviously have t_{ϵ}^{∗} ≤ t^{∗}, and the optimal solution of the approximated problem (Q_{ϵ}^{∗}, σ_{ϵ}^{∗}) is not necessarily feasible for the original problem. However, Definition 2.1 of a ϵ-approximation of a second-order cone implies in our case that

∥(q_{ϵ,i1}^{∗}, …, q_{ϵ,ik}^{∗})∥ ≤ (1 + ϵ)σ_{ϵ,i}^{∗} ∀1 ≤ i ≤ n ,

which means that (Q_{ϵ}^{∗}, (1 + ϵ)σ_{ϵ}^{∗}) is feasible for the original problem, with an objective value equal to (1 + ϵ)t_{ϵ}^{∗}. Since we must have then t^{∗} ≤ (1 + ϵ)t_{ϵ}^{∗}, we conclude that

\frac{t^{∗}}{1 + ϵ} ≤ t_{ϵ}^{∗} ≤ t^{∗} ⇔ 0 ≤ \frac{t^{∗} − t_{ϵ}^{∗}}{t^{∗}} ≤ \frac{ϵ}{1 + ϵ} ,

i.e. we have a bound on the relative accuracy of our approximated optimum objective value.

We generated three random problems for each of the 11 combinations of parameters (p, q, k) presented in Table 3, and report in Table 5 the average computing time in seconds. Each column in this table corresponds to a different way to solve the truss-topology design problem:

Table 4: Dimensions of the approximated problems (primal above, dual below).

\begin{matrix}p × q × k & ϵ = 10^{−2} & ϵ = 10^{−5} & ϵ = 10^{−8} \\ 2 × 2 × 2 & 35 × 58 & 60 × 108 & 85 × 158 \\ 2 × 2 × 4 & 85 × 146 & 160 × 296 & 235 × 446 \\ 2 × 2 × 8 & 200 × 352 & 370 × 692 & 545 × 1042 \\ 2 × 2 × 16 & 420 × 744 & 795 × 1494 & 1165 × 2234 \\ 2 × 2 × 32 & 860 × 1528 & 1635 × 3078 & 2410 × 4628 \\ 4 × 4 × 2 & 798 × 1188 & 1368 × 2328 & 1938 × 3468 \\ 4 × 4 × 4 & 1938 × 3060 & 3648 × 6480 & 5358 × 9900 \\ 4 × 4 × 6 & 3306 × 5388 & 6156 × 11088 & 9006 × 16788 \\ 6 × 6 × 2 & 4305 × 6270 & 7380 × 12420 & 10455 × 18570 \\ 6 × 6 × 4 & 10455 × 16230 & 19680 × 34680 & 28905 × 53130 \\ 8 × 8 × 2 & 13916 × 20104 & 23856 × 39984 & 33796 × 59864 \end{matrix}

\begin{matrix}p × q × k & ϵ = 10^{−2} & ϵ = 10^{−5} & ϵ = 10^{−8} \\ 2 × 2 × 2 & 43 × 65 & 68 × 115 & 93 × 165 \\ 2 × 2 × 4 & 101 × 155 & 176 × 305 & 251 × 455 \\ 2 × 2 × 8 & 232 × 365 & 402 × 705 & 577 × 1055 \\ 2 × 2 × 16 & 484 × 765 & 859 × 1515 & 1229 × 2255 \\ 2 × 2 × 32 & 988 × 1565 & 1763 × 3115 & 2538 × 4665 \\ 4 × 4 × 2 & 846 × 1482 & 1416 × 2622 & 1986 × 3762 \\ 4 × 4 × 4 & 2034 × 3534 & 3744 × 6954 & 5454 × 10374 \\ 4 × 4 × 6 & 3450 × 6042 & 6300 × 11742 & 9150 × 17442 \\ 6 × 6 × 2 & 4425 × 7995 & 7500 × 14145 & 10575 × 20295 \\ 6 × 6 × 4 & 10695 × 19065 & 19920 × 37515 & 29145 × 55965 \\ 8 × 8 × 2 & 14140 × 25844 & 24080 × 45724 & 34020 × 65604 \end{matrix}

Table 5: Computing times to solve truss-topology problems using different approaches.

QCO  SOCO  10^{−2}  10^{−5}  10^{−8}
p × q × k  (D) (P) (D)  (D’) (P) (D) (P) (D)  (P) (D)
2 × 2 × 2  0.00 0.00 0.00  0.58 0.01 0.00 0.01 0.02  0.04 0.03
2 × 2 × 4  0.01 0.01 0.00  0.10 0.02 0.02 0.05 0.05  0.11 0.12
2 × 2 × 8  0.01 0.01 0.01  0.12 0.05 0.05 0.14 0.14  0.33 0.34
2 × 2 × 16 0.01 0.02 0.01  0.19 0.10 0.11 0.37 0.37  0.85 0.83
2 × 2 × 32 0.02 0.03 0.03  0.35 0.26 0.27 0.87 0.90  1.96 1.86
4 × 4 × 2  0.09 0.06 0.29  0.60 0.30 0.21 0.73 0.69  1.54 1.52
4 × 4 × 4  0.11 0.13 0.74  1.74 1.52 0.92 3.19 2.86  9.30 5.61
4 × 4 × 6  0.20 0.26 2.22  5.03 3.86 2.09 7.98 5.50 13.50 10.84
6 × 6 × 2  0.47 0.61 1.96 30.13 2.54 1.88 21.59 5.16 34.95 10.45
6 × 6 × 4  1.28 1.32 6.48 475.73 19.08 7.82 40.80 23.89 396.04 43.43
8 × 8 × 2  4.30 2.76 11.81 339.66 12.08 8.89 53.29 24.55 127.68 48.42

a. the first column reports computing times using MOSEK on the quadratically constrained formulation (TQC),

b. the following three columns report computing times using MOSEK on the primal and the dual second-order cone formulations (TP)–(TD) in columns (P) and (D), as well as the results of SeDuMi on the dual formulation in column (D’).

c. the last six columns report computing times using the interior-point code in MOSEK to solve the polyhedral approximations of the primal and dual second-order cone problems (TP)–(TD) with three different accuracies.

Our first constatation is that solving the quadratically constrained formulation (TQC) and the primal second-order cone formulation (TP) directly are the two fastest methods (with similar computing times). The quadratically constrained formulation has less variables and constraints, but this advantage seems to be counterbalanced by a more efficient second-order cone solver.

Solving the dual second-order cone formulation (TP) directly is also very fast with a 2 × 2 nodal grid but noticeably slower on the larger problems (3 to 8 times slower). This is most probably due to the greater dimensions of the problem. The SeDuMi solver is much less efficient to solve these dual problems, and is really slow on the three largest problems.

Let us now look at the approximated problems. First of all, we checked whether the accuracy of the optimum approximated objective was below the theoretical bound \frac{ϵ}{1+ϵ}, since rounding errors in the computations and handling of irrational coefficients could affect this result. Unsurprisingly, the accuracy was below the theoretical threshold for all experiments.

Computing times are worse than for the direct approaches, even with the low accuracy 10^{−2}. The difference grows up to one or two orders of magnitude for the larger problems. We also observe that despite slightly greater dimensions, solving the approximated dual problem is more efficient than solving the primal problem, especially with the largest problems. The reasons for this behaviour, which is opposite to the situation for direct resolutions, are unclear to us, but could be related to sparsity issues.

Finally, let us mention that we also tried to solve the linear approximations using the simplex algorithm instead of an interior-point method. This lead to surprisingly bad computing times: for example, solving problem 4 × 4 × 4 with accuracy 10^{−2} using the MOSEK^{7} simplex code took 21.57 seconds, instead of 1.52 seconds with the interior-point algorithm. We believe this disastrous behaviour of the simplex algorithm is due to the presence of an exponential number of vertices in the approximation, which leads to very slow progress.

3.3 Quadratic optimization

Second-order cone formulations of truss-topology design problems feature a relatively large number of small cones. Since our approximation procedure has not proven to be more efficient than direct methods on these problems, we would like to turn our attention to the opposite configuration, i.e. a small number of large cones. We are going to show that convex quadratic optimization can be formulated such as to meet this requirement.

More specifically, we are going to consider linearly constrained convex quadric optimization problems. Such problems can be formulated as

min \frac{1}{2}x^{T} Qx + c^{T} x + c_{0} s.t. l_{c} ≤ Ax ≤ u_{c} and l_{x} ≤ x ≤ u_{x} ,  (QO)

where x ∈ ℝ^{n} denotes the vector of design variable. The objective is defined by a matrix Q ∈ ℝ^{n×n}, required to be positive semidefinite to ensure convexity of the problem, a vector c ∈ ℝ^{n} and a scalar c_{0} ∈ ℝ. Variables are bounded by two vectors l_{x} ∈ ℝ^{n} and u_{x} ∈ ℝ^{n} (note that some components of l_{x} or u_{x} can be equal to −∞ or +∞ if a variable has no lower or upper bound). Finally, the linear constraints are described by a matrix A ∈ ℝ^{m×n} and two vectors l_{c} ∈ ℝ^{m} and u_{c} ∈ ℝ^{m} (with the same remark holding about possible infinite values for some components of l_{c} and u_{c}).

In order to model problem (QO) with a second-order cone formulation, we first write the Cholevsky factorization of matrix Q. Indeed, we have Q = L^{T} L with L ∈ ℝ^{k×n}, where k ≤ n is the rank of Q. Introducing a vector of auxiliary variables z ∈ ℝ^{k} such that z = Lx, we have that x^{T} Qx = x^{T} L^{T} Lx = (Lx)^{T} Lx = z^{T} z, which allows us to write the following problem:

min \frac{r + v}{2} +c^{T} x+c_{0} s.t. (r, v, z) ∈ 𝕃^{k+1}, r−v = 1, l_{c} ≤ Ax ≤ u_{c} and l_{x} ≤ x ≤ u_{x} . (QO’)

It is indeed equivalent to (QO), since the conic constraint (r, v, z) ∈ 𝕃^{k+1} combined with the equality r − v = 1 leads to

v^{2} + ∑_{i=1}^{k} z_{i}^{2} ≤ r^{2} ⇔ ∑_{i=1}^{k} z_{i}^{2} ≤ r^{2} − v^{2} ⇔ z^{T} z ≤ (r − v)(r + v) ⇔ z^{T} z ≤ r + v ,

which is why the quadratic term \frac{1}{2}x^{T} Qx in the objective of (QO) could be replaced by \frac{r+v}{2} in (QO’).

The problems we tested come from the convex quadratic optimization library QPDATA, collected by Maros and Mészáros [MM99]. As for our tests with truss-topology design problems, we decided to formulate approximations with three different accuracies 10^{−2}, 10^{−5} and 10^{−8}. Table 6 lists for each problem its original size (variables × constraints), the number of nonzero elements in the constraint matrix A, the number of nonzero elements in the upper triangular part of Q and the size (variables × constraints) of each of the three polyhedral approximations.

7In order to make sure that this behaviour was not caused by a flaw in the MOSEK simplex solver, we performed a similar comparison with the CPLEX solver, which lead to the same conclusion.

Table 7 reports computing times (in seconds) needed to solve these convex quadratic optimization problems in three different ways:

a. the first column reports computing times using MOSEK directly on the original quadratic formulation (QO),

b. the following two columns report computing times needed to solve the second-order cone formulation (QO’) of these problems with MOSEK and SeDuMi (in columns labelled (SOCO) and (SOCO’) respectively),

c. the last three columns report computing times using MOSEK to solve the polyhedral approximations of the second-order cone problem (QO’) with three different accuracies.

Once again, the direct approach, i.e. solving (QO), is the most efficient method. Solving these problems with a second-order cone formulation is slower, especially on larger problems. Using SeDuMi instead of MOSEK degrades further the computing times.

It is also manifest that the linear approximations take more time than the direct approach to provide a solution, even with the lowest accuracy 10^{−2} (however, we note that this low accuracy approximation is faster than the SeDuMi resolution on a few problems).

Although we only tested small-scale and medium-scale problems, it is pretty clear from the trend present for the last problems that large-scale problems would also be most efficiently solved directly as convex quadratic optimization problems.

We pointed out in Section 2.9 that our bound on the accuracy of the polyhedral approximation did not imply anything on the quality of the optimum of the approximated problems. Indeed, Table 8 reports the relative accuracy for a few representative approximated problems. Some problems (GENHS28, DUALC5) behave very well, with a relative accuracy well below the target accuracy. Other problems (GOULDQP2, MOSARQP1) have higher relative accuracies, but still decreasing when the target accuracy is decreased. Problem CVXQP3S shows a worse behaviour, with virtually no improvement between the second and the third approximation. Finally, HS21 is a toy problem with the surprising property that its approximation is exact for any accuracy.

We were able to compute these relative accuracies because the true optimal objective values were known by other means. In a real-world situation where such a piece information would not be available, it would be still possible to estimate roughly this accuracy. Indeed, since our approximation is a relaxation, we have that the approximated optimal objective p_{ϵ}^{∗} is lower than the true optimum p^{∗}. On the other hand, the optimal solution x_{ϵ}^{∗} of the approximation must be feasible, since it satisfies the linear constraints. Computing the objective function corresponding to this solution, i.e. letting p′_{ϵ}^{∗} = \frac{1}{2}x_{ϵ}^{∗}^{T} Qx_{ϵ}^{∗} + c^{T} x_{ϵ}^{∗} + c_{0}, we have finally that p_{ϵ}^{∗} ≤ p^{∗} ≤ p′_{ϵ}^{∗}, which allows us to estimate a posteriori the true optimum objective value.

However, in the special case where the objective is purely quadratic, i.e. when c = 0 and c_{0} = 0, it is possible to slightly modify the formulation so that we have a bound on the accuracy of the objective^{8}. Indeed, letting again z = Lx, we add this time the conic constraint

8 A similar improvement can be made in the case when c belongs to the column space of L^{T}, the Cholevsky factor of Q. However, we have been unable to generalize this construction to the case where c does not belong to this column space, e.g. for an objective equal to x_{1}^{2} + x_{2}.

Table 6: Statistics for the convex quadratic optimization problems.

Name  Size ANZ QNZ  Size 10^{−2}  Size 10^{−5}  Size 10^{−8}
TAME  2 × 1  2  3  9 × 13  14 × 23  19 × 33
HS21  2 × 1  2  2  14 × 22  24 × 42  34 × 62
ZECEVIC2  2 × 2  4  1  9 × 14  14 × 24  19 × 34
HS35  3 × 1  3  5  20 × 31  35 × 61  50 × 91
HS35MOD  3 × 1  3  5  20 × 31  35 × 61  50 × 91
HS52  5 × 3  7  7  29 × 46  49 × 86  69 × 126
HS76  4 × 3  10  6  28 × 46  48 × 86  68 × 126
HS51  5 × 3  7  7  29 × 46  49 × 86  69 × 126
HS53  5 × 3  7  7  29 × 46  49 × 86  69 × 126
S268  5 × 5  25  15  34 × 57  59 × 107  84 × 157
HS268  5 × 5  25  15  34 × 57  59 × 107  84 × 157
GENHS28  10 × 8  24  19  61 × 100  106 × 190  151 × 280
LOTSCHD  12 × 7  54  6  47 × 70  76 × 128  106 × 188
HS118  15 × 17  39  15  99 × 169  174 × 319  248 × 467
QPCBLEND  83 × 74 491  83  545 × 914  959 × 1742  1373 × 2570
CVXQP2S  100 × 25  74 386 647 × 1020  1135 × 1996  1624 × 2974
CVXQP1S  100 × 50 148 386 647 × 1045  1135 × 2021  1624 × 2999
CVXQP3S  100 × 75 222 386 647 × 1070  1135 × 2046  1624 × 3024
QPCBOEI2 143 × 166 1196 143 939 × 1614  1652 × 3040  2366 × 4468
DUALC5  8 × 278 2224  36  54 × 61  94 × 441  134 × 521
PRIMALC1  230 × 9 2070 229 1505 × 2329  2646 × 4611  3789 × 6897
PRIMALC5  287 × 8 2296 286 1878 × 2903  3304 × 5755  4732 × 8611
DUAL4  75 × 1  75 2799  493 × 761  867 × 1509  1241 × 2257
GOULDQP2 699 × 349 1047 697 2635 × 3872  4370 × 7342 6107 × 10816
DUAL1 85 × 1 85 3558 558 × 861 982 × 1709 1406 × 2557
PRIMALC8  520 × 8 4160 519 3410 × 5268 5996 × 10440 8586 × 15620
GOULDQP3 699 × 349 1047 1395 4584 × 7420 8062 × 14376 11546 × 21344
DUAL2 96 × 1 96 4508 632 × 976 1110 × 1932 1589 × 2890
MOSARQP2 900 × 600 2390 945 5911 × 9721 10395 × 18689 14887 × 27673

Table 7: Computing times to solve convex quadratic optimization problems

Name  QO SOCO SOCO’ 10^{−2} 10^{−5} 10^{−8}
TAME  0.06  0.11  0.46 0.01 0.01 0.02
HS21  0.00  0.01  0.17 0.01 0.01 0.02
ZECEVIC2 0.01  0.00 0.101 0.00 0.00 0.03
HS35 0.00 0.00 0.1 0.01 0.01 0.02
HS35MOD 0.00  0.03 0.361 0.01 0.01 0.02
HS52  0.00  0.01  0.09 0.01 0.02 0.03
HS76  0.00  0.00  0.11 0.00 0.02 0.04
HS51  0.00  0.01  0.08 0.00 0.02 0.04
HS53  0.01  0.00  0.09 0.00 0.02 0.03
S268  0.02  0.01  0.20 0.01 0.04 0.07
HS268  0.01  0.00  0.20 0.01 0.04 0.07
GENHS28 0.00  0.01  0.08 0.01 0.03 0.07
LOTSCHD 0.02  0.01  0.16 0.02 0.03 0.08
HS118  0.00  0.01  0.25 0.05 0.11 0.17
QPCBLEND 0.03  0.12  0.83 0.41 1.43 2.75
CVXQP2S 0.02  0.23  1.34 0.58 1.91 3.58
CVXQP1S 0.03  0.14  1.55 0.57 1.88 4.47
CVXQP3S 0.05  0.21  1.86 0.60 2.32 3.93
QPCBOEI2 0.08  0.59  4.01 1.26 5.53 11.55
DUALC5  0.01  0.01  1.56 0.04 0.04 0.09
PRIMALC1 0.04  0.91  2.17 1.27 7.93 17.85
PRIMALC5 0.03  0.83  0.69 3.15 9.49 20.83
DUAL4  0.04  0.13  0.57 0.39 0.83 1.42
GOULDQP2 1.21  0.41  1.87 3.54 10.14 20.65
DUAL1  0.06  0.18  0.72 0.72 1.66 2.53
PRIMALC8 0.08  4.72  4.32 3.97 30.32 50.99
GOULDQP3 0.98 12.16  6.32 9.16 38.44 74.15
DUAL2 0.07 0.16 0.93 0.88 2.03 3.03
MOSARQP2 0.38  1.64 8.512 7.69 39.03 76.59

Table 8: Relative accuracy of the optimum of some approximated problems.

Name  10^{−2}  10^{−5}  10^{−8}
GENHS28  5.2e-4 8.8e-7 7.6e-10
HS21  3.3e-10 3.3e-10 3.3e-10
CVXQP3S  6.7e-3 7.2e-4 4.7e-4
DUALC5  3.9e-3 7.2e-7 3e-10
GOULDQP2 5.2e-1 5.6e-3 2.1e-5
MOSARQP2 8.6e-1 2.5e-4 8.5e-7

(r, z) ∈ 𝕃^{k}, which implies to z^{T} z ≤ r^{2} ⇔ x^{T} Qx ≤ r^{2}. We can now choose r as our objective, which is equivalent to minimizing \sqrt{x^{T} Qx} and is obviously the same thing as minimizing the true quadratic objective \frac{1}{2}x^{T} Qx. This leads to a situation that is very similar to the case of truss-topology design problems, and one can show without difficulties that this approximated problem provides an estimation of the true optimum with a relative accuracy equal to (\frac{ϵ}{1+ϵ})^{2}.

4 Concluding remarks

In this paper, we presented a polyhedral approximation of the second-order cone originally developed by Ben-Tal and Nemirovski [BTN98]. Our presentation features several improvements, including smaller dimensions for the approximation, a more transparent proof, complete developments valid for any size of the second-order cone (i.e. not limited to powers of two) and explicit constants in the derivation of a theoretical bound on the size of the approximation (Theorem 2.4).

This scheme was implemented in MATLAB and optimized as much as possible. Indeed, we developed several approaches to reduce the size of the resulting linear problems (including pivoting out some variables and using dynamic programming to choose the best accuracies for each stage of the decomposition). Our experiments mainly showed that solving the original second-order cone problems or alternative equivalent formulations is more efficient than solving the linear approximations, even at low accuracies. On a side note, we noticed that these approximate problems are particularly difficult to solve with the simplex algorithm.

However, we would like to point out this approximating scheme can still prove very useful in certain well-defined circumstances, such as a situation where a user is equipped with a solver that is only able to solve linear optimization problems. In this case, this procedure provides him with an inexpensive and relatively straightforward way to test improved versions of his linear models that make use of second-order cones.

Moreover, we have to admit that we tested two very specific classes of second-order cone optimization problems for which either a simplified formulation or a well-understood dedicated algorithm was available. It might well be possible that this linearizing scheme becomes competitive for other types of difficult (i.e. that cannot be simplified and for which no dedicated solver is available) second-order cones optimization problems.

We would also like to insist on the fact that it is not possible to guarantee a priori the accuracy of a linear approximation of a general second-order optimization problem (see the example in Section 2.9). It is nevertheless possible to provide such a bound in some special cases (e.g. truss-topology design problems or convex quadratic optimization problems with a pure quadratic objective).

It is worth to point out that a straightforward modification of our polyhedral approximation of 𝕃^{2} can lead to a restriction instead of a relaxation of second-order cone optimization problems. This would then provide an upper bound instead of a lower bound on the true optimum objective value, and optimal solutions of the approximate problems would always be feasible for the original problem. However, this approach can be problematic in some cases since it might happen that the approximated problem is unfeasible, even if the original problem admits some feasible solutions.

An interesting topic for further research is the generalization of the polyhedral approximation of 𝕃^{2} or, more precisely, of the unit ball B_{2}(1), to other convex sets. Indeed, finding a similar polyhedral approximation for a set like {(x_{1}, x_{2}) ∈ ℝ^{2} | |x_{1}|^{p} + |x_{2}|^{p} ≤ 1} with p > 1, i.e. the unit ball for the p-norm, would lead to linearizing scheme for other classes of convex problems, such as l_{p}-norm optimization. However, it is unclear to us at this stage whether this goal is achievable or not, since the symmetry of the standard unit ball, which is not present for other norms, seems to play a great role in the construction of the approximation.

Acknowledgements. The author wishes to thank Professor Kees Roos and the SSOR department at the Delft University of Technology for their kind hospitality, as well as Erling Andersen for making the MOSEK solver available to him.

References

[AA99] E. D. Andersen and K. D. Andersen, The MOSEK interior point optimizer for linear programming: an implementation of the homogeneous algorithm, High Performance Optimization (H. Frenk, C. Roos, T. Terlaky, and S. Zhang, eds.), Applied optimization, vol. 33, Kluwer Academic Publishers, 1999.

[ART00] E. D. Andersen, C. Roos, and T. Terlaky, On implementing a primal-dual interiorpoint method for conic quadratic optimization, in preparation, 2000.

[BTN94] A. Ben-Tal and A. Nemirovski, Potential reduction polynomial-time method for truss topology design, SIAM Journal of Optimization 4 (1994), 596–612.

[BTN98] A. Ben-Tal and A. Nemirovski, On polyhedral approximations of the second-order cone, Tech. report, Minerva Optimization Center, Faculty of Industrial Engineering and Management, Technion - Israel Institute of Technology, Haifa, Israel, 1998, to appear in Mathematics of Operations Research.

[LVBL98] M. S. Lobo, L. Vandenberghe, S. Boyd, and H. Lebret, Applications of second-order cone programming, Linear Algebra and its Applications 284 (1998), 193–228.

[MM99] I. Maros and Cs. Mészáros, A repository of convex quadratic programming problems, Optimization Methods and Software 11-12 (1999), 671–681, special issue on interior-point methods (CD supplement with software), guest editors: Florian Potra, Cornelis Roos and Tamás Terlaky.

[NN94] Y. E. Nesterov and A. S. Nemirovski, Interior-point polynomial methods in convex programming, SIAM Studies in Applied Mathematics, SIAM Publications, Philadelphia, 1994.

[Stu99a] J. F. Sturm, Duality results, High Performance Optimization (H. Frenk, C. Roos, T. Terlaky, and S. Zhang, eds.), Applied optimization, vol. 33, Kluwer Academic Publishers, 1999, pp. 21–60.

[Stu99b] J. F. Sturm, Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones, Optimization Methods and Software 11-12 (1999), 625–653, special issue on interior-point methods (CD supplement with software), guest editors: Florian Potra, Cornelis Roos and Tamás Terlaky.

[SW70] J. Stoer and Ch. Witzgall, Convexity and optimization in finite dimensions I, Springer Verlag, Berlin, 1970.