I want to solve the (convex) optimisation task:

$max_{r,z}\quad r$
subject to the following two constraints
$r\|x_i\| - x_i^Tz \leq 0 \qquad \forall i=1,\dots, N $
$\|z\| \leq 1$

$r$ is a scalar, $z$ is a vector, the $x_i$'s are vectors of the same dimension and $\|\cdot\|$ is the simple eucl. norm. One can assume that the feasible region is non-empty.

Is there an easy way to solve this? I think this should be easy because without the $\|z\|\leq 1$ constraint this is just a linear program. Before refering my to software packages can you give an hint about the general approach which is usefull for this kind of task?
Thanks DG


You have a couple options, depending on how crucial it is that you take the Euclidean norm of $z$.

  1. Use your formulation as is, with a small tweak:

    $max_{r,z}\quad r$
    subject to the following two constraints
    $r\|x_i\| - x_i^Tz \leq 0 \qquad \forall i=1,\dots, N $
    $z^{T}z \leq 1$

    This problem is a quadratically-constrained program, for which there are many fast solvers out there, such as CPLEX and Gurobi. This particular program is also a second-order cone program, semidefinite program, and convex nonlinear program, so you could also use any of those solvers as well. The reason I replaced the Euclidean norm constraint with a dot product is that the two constraints are equivalent, but the latter is differentiable, whereas the former is not. Nondifferentiable functions require more expensive algorithms, and this problem doesn't require that type of machinery, so it's best to avoid it.

  2. Change the norm of $z$. The 1-norm and infinity-norm are both linear functions of the elements of $z$, and replacing the Euclidean norm in your formulation with either of those norms results in a linear program, for which the best solvers tend to be commercial (Gurobi, CPLEX), but slower free solvers exist (GLPK, solvers in the COIN-OR suite). Using the 1-norm means that any solution to that formulation would also be a feasible solution of your current formulation (that is, using the 1-norm would yield a restriction of your current formulation). Using the infinity-norm means that any solution to your current formulation would also be feasible in the infinity-norm formulation (that is, using the infinity-norm would yield a relaxation of your current formulation).

Although it is true that linear programming solvers are very efficient, I would select option 1, because quadratically-constrained programming solvers are also very efficient (relative to convex programming solvers, and other types of nonlinear programming solvers) and can solve large formulations (at least hundreds of thousands of decision variables, last time I looked at the literature). Unless your formulation were astonishingly large, you should be fine using a quadratically-constrained programming solver in serial, and needn't change the norm in your formulation unless you absolutely have to.

A final remark: I would scale the vectors $x_{i}$ prior to constructing your formulation so that they all have unit norm, which will likely help with the conditioning of your formulation when you solve it numerically. It's another trick that costs you practically nothing, but guards you against numerical difficulties.

| cite | improve this answer | |
  • $\begingroup$ You should not call this a semidefinite program. It is a quadratically constrained program, or slightly more general a second-order cone program. Calling it a semidefinite program would be akin to calling a linear program a second order cone program (following the classes linear program - quadratically constrained program - second order cone program - semidefinite program) $\endgroup$ – Johan Löfberg Aug 9 '13 at 12:15
  • $\begingroup$ You're right that calling it a QP would be better. For some reason I saw a Gramian (even a small one) and jumped at the pattern. $\endgroup$ – Geoff Oxberry Aug 9 '13 at 16:13
  • $\begingroup$ I think quadratic program sends the wrong signals, it is typically limited to problems with quadratic objective and linear constraints. quadratically constrained problem, or, perhaps most common but slightly more general than the problem here, a second order cone program. $\endgroup$ – Johan Löfberg Aug 9 '13 at 16:31

As we see in Geoffs answer, this is a very simple quadratically constrained problem, or more generally a second-order cone program. If you don't have extreme performance requirements or enormous dimensions, solving it using a standard nonlinear solver in the quadratic form $z^Tz \leq 1$ or using an SOCP solver in the norm formulation $\|z\|\leq 1$ will work perfectly fine.

If you have to improve performance, there are methods to exploit the single cone feature. Here is one example

SIAM J. Optim., 17(2), 459–484. (26 pages) An Active Set Method for Single-Cone Second-Order Cone Programs E. Erdougan and G. Iyengar

I would like to point out that replacing the norm with a 1-norm probably won't work well. The quadratic norm has its origin in the geometric background of this problem (which I interpret as finding a vector which has the smallest angle to a given set of vectors).

Interestingly, a QP approximation of the problem seems to work extremely well. Remove the quadratic constraint, and add a penalty $\alpha z^Tz$ to the objective instead. I would not be surprised if it is possible to prove something about this.

In the code below, implemented using YALMIP (Disclaimer, developed by me) in MATLAB, using CPLEX as the solver, the average distance from the true $z$ and the $z$ computed using the QP heuristics, is on the order of $10^{-6}$, while the distance to the solution from the LP (1-norm) formulation is on the order $10^{-1}$.

z = sdpvar(5,1);
r = sdpvar(1);

err1 = [];
err2 = [];
for i = 1:1000
    X = randn(5,10);
    Con = [r*sqrt(sum(X.^2,1)) <= z'*X,norm(z,2) <= 1]
    sol = solvesdp(Con,-r)
    if sol.problem == 0 & double(r)>1e-3
        zSOCP = double(z);
        Con = [r*sqrt(sum(X.^2,1)) <= z'*X];
        sol = solvesdp(Con,-r+0.001*z'*z);
        zQP = double(z/norm(double(z)));
        err1 = [err1 norm(zQP-zSOCP)];
        Con = [r*sqrt(sum(X.^2,1)) <= z'*X, norm(z,1)<=1];
        sol = solvesdp(Con,-r);
        zLP = double(z/norm(double(z)));
        err2 = [err2 norm(zLP-zSOCP)];

Finally, using geometric insight might lead to a much better approach to solve this problem. You are essentially looking for a particularly defined center of a set of points on the unit-sphere.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.