Linear Optimization with Quadratic Constraints: Efficient, arbitrary precision solver

Specific Domains Optimization (Mathematical)
, ,
1

I want to solve a linear optimization problem with quadratic constraints of the form
\max_{\boldsymbol x} y \text{ such that } f_i(y \boldsymbol x) \leq 0 \: \forall \: i = 1, \dots, C
and the f_i are quadratic in the argument, i.e., it is of form f_i(y \boldsymbol x) = \alpha_i + \sum_{j=1}^N \sum_{k=1}^N y \beta_{kl} x_j x_k.

Since the number of constraints C is in the order of thousands, I want to solve this problem rather efficiently. I tried using Convex.jl and a certain optimizer, but the equivalent cone optimization problem that Convex.jl generates is usually significantly larger in the sense that a fistful optimization variables get blown up to multiple thousand ones.

I think now about using JuMP.jl but I am not sure if it supports arbitrary precision (through the use of BigFloat, preferably).

2

JuMP doesn’t support BigFloat, only Float64.

But note that your problem isn’t linear, or quadratic. It seems to have cubic terms? That means it’s nonconvex and nonlinear. You could use Ipopt though and see what sort of solution you get.

1 Like
3

Yeah you are right, as written now it’s cubic. I could reformulate the problem maybe as
\max_{\boldsymbol x} \Vert \boldsymbol x \Vert \text{ such that } f_i(y \boldsymbol x) \leq 0 \: \forall \: 1 , \dots, C \text{ for } y \text{ given } . The important thing is that the \boldsymbol x satisfies the constraints, the actual objective is not that important.

Does Ipopt.jl work with BigFloat?

4

No, sorry I meant Ipopt and JuMP together.

Why do you need BigFloat?

5

The constraints are poorly conditioned and solving the problem on Matlab with SDPT3 has shown that with default precision, not more than ~14 variables may be optimized at once, which is much less then I want to optimize.

This is why considered julia (besides Mathematica) since it offers (in theory) arbitrary precision arithmetics.

I will now give Optim.jl a try, based on the suggestion made in this post.

6

How did you solve it with SDPT3? The problem isn’t convex? Did you fix y, and then solve a SDP relaxation? That would cause numerical issues.

It might help if you could provide a small set of data to test with. What is N (in practice)? What is α? What is β? Why doesn’t the summation in f depend on i? What is l?

7

The way it is done in Matlab with CVX/SDPT3 is to formulate the problem as
\min_{\boldsymbol x} \max_i \vert g_i(\boldsymbol x) \vert where the g_i are linear but complex in (\boldsymbol x) . Since my desired outcome is that \max_i \vert g_i(\boldsymbol x) \vert \leq 1, I thought about re-formulating this as \max_i \vert g_i(\boldsymbol x) \vert^2 \leq 1 \Leftrightarrow \vert g_i(\boldsymbol x) \vert^2 \leq 1 \: \forall i with f_i basically being the \vert g_i(\boldsymbol x) \vert^2.

8

Are you able to share the problem written in Matlab?

1 Like
9

What is the formulation you actually want to solve, before you reformulate it in different ways?

10

If it turns out you don’t need arbitrary precision after all then you might want to try Gurobi. The latest version can handle large nonconvex quadratic problems. It may be helpful even if your problem is too large to solve to the global optimum before the end of the universe because it provides bounds on the global optimum. This can help you evaluate local optima from IPOPT or Gurobi’s own solve sequence.

Here’s a presentation that explains how the algorithm works:

2 Likes
11

The base problem is:

\max_{\boldsymbol x} y \text{ such that } \max_i \vert g_i(y \boldsymbol x) \vert \leq 1
where y \in \mathbb R^+, \boldsymbol x \in \mathbb R, g_i: \mathbb R \mapsto \mathbb C and g_i are linear in \boldsymbol x.

For fixed y, this is a convex problem in \boldsymbol x which I solve with CVX.

To optimize y, I apply bisection in an outer loop since the g_i are essentially complex-valued polynomials with coefficients \boldsymbol x, i.e., \max_i \vert g_i(y \boldsymbol x) \vert is expexted to decrease if y decreases. Thus, a bisection method is applied to determine y. Details can be found in this paper.

12

If the g_i are linear in x, then g_i(yx) = y g_i(x) and the constraint can be rewritten as

-1 \leq y g_i(x) \leq 1

which is itself equivalent to

-z \leq g_i(x) \leq z

where z = 1/y (under the assumption that y > 0; note that you gave us y \geq 0 and any x is feasible for y=0).

Now you can solve the LP

\min z subject to -z \leq g_i(x) \leq z

using any LP solver.

13

Recall that the g_i(\boldsymbol x) map to the complex numbers \mathbb C and thus \vert g_i(\boldsymbol x)\vert \leq 1 is not equivalent to -1 \leq g_i(\boldsymbol x) \leq 1.

I could try to formulate this condition using polar numbers, maybe I will look into that.

14

Maybe something like this?

using JuMP, Ipopt
model = Model(Ipopt.Optimizer)
@variable(model, 1 <= x[1:3] <= 2)
@variable(model, y >= 0, start = 1)
@objective(model, Max, y)
@expression(model, g_re[i=1:3], x[i]) # Replace these with g_i
@expression(model, g_im[i=1:3], x[i]) #
@NLconstraint(model, [i=1:3], g_re[i]^2 + g_im[i]^2 <= 1 / y^2)
optimize!(model)

If Ipopt can’t find a solution, you could get rid of y and do the bisection search.

1 Like
15

But can’t you simplify the NLConstraint to a (convex, quadratic) Constraint by replacing 1/y^2 with z and changing the objective to \min z?

1 Like