# Linear Optimization with Quadratic Constraints: Efficient, arbitrary precision solver

**URL:** <https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105>\
**Category:** Optimization (Mathematical)\
**Tags:** question, optimization, extended-precision\
**Created:** [June 2, 2022, 6:51am UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105 "2022-06-02T06:51:19Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![DanDoe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dandoe/32/52717_2.png) [@DanDoe](https://discourse.julialang.org/u/DanDoe)\
**Post date:** [June 2, 2022, 6:51am UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/1 "2022-06-02T06:51:20Z")

</div>

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`](https://jump.dev/Convex.jl/stable/) 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`](https://jump.dev/JuMP.jl/stable/) but I am not sure if it supports arbitrary precision (through the use of `BigFloat`, preferably).

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 2, 2022, 7:01am UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/2 "2022-06-02T07:01:36Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![DanDoe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dandoe/32/52717_2.png) [@DanDoe](https://discourse.julialang.org/u/DanDoe)\
**Post date:** [June 2, 2022, 7:13am UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/3 "2022-06-02T07:13:45Z")

</div>

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`?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 2, 2022, 7:14am UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/4 "2022-06-02T07:14:31Z")

</div>

No, sorry I meant Ipopt and JuMP together.

Why do you need BigFloat?

---

<div class="post-metadata">

**Author:** ![DanDoe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dandoe/32/52717_2.png) [@DanDoe](https://discourse.julialang.org/u/DanDoe)\
**Post date:** [June 2, 2022, 7:28am UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/5 "2022-06-02T07:28:48Z")

</div>

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`](https://julianlsolvers.github.io/Optim.jl/stable/#) a try, based on the [suggestion made in this post](https://scicomp.stackexchange.com/questions/34495/arbitrary-precision-optimization-libraries/34496#34496).

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 2, 2022, 9:13am UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/6 "2022-06-02T09:13:41Z")

</div>

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`?

---

<div class="post-metadata">

**Author:** ![DanDoe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dandoe/32/52717_2.png) [@DanDoe](https://discourse.julialang.org/u/DanDoe)\
**Post date:** [June 2, 2022, 9:52am UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/7 "2022-06-02T09:52:46Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [June 2, 2022, 2:05pm UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/8 "2022-06-02T14:05:17Z")

</div>

Are you able to share the problem written in Matlab?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 2, 2022, 7:24pm UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/9 "2022-06-02T19:24:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![NiclasMattsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/niclasmattsson/32/21988_2.png) [@NiclasMattsson](https://discourse.julialang.org/u/NiclasMattsson)\
**Post date:** [June 2, 2022, 8:49pm UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/10 "2022-06-02T20:49:48Z")

</div>

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:

[![](https://global.discourse-cdn.com/julialang/original/3X/6/1/6107a213ebf212cee6508f421cc33b05d008acff.jpeg "Non-Convex Quadratic Optimization Webinar") ](https://www.youtube.com/watch?v=0g5cMvOV7KY)

---

<div class="post-metadata">

**Author:** ![DanDoe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dandoe/32/52717_2.png) [@DanDoe](https://discourse.julialang.org/u/DanDoe)\
**Post date:** [June 5, 2022, 11:47am UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/11 "2022-06-05T11:47:29Z")

</div>

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.](https://msp.org/camcos/2012/7-2/p04.xhtml)

---

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [June 5, 2022, 1:36pm UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/12 "2022-06-05T13:36:52Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![DanDoe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dandoe/32/52717_2.png) [@DanDoe](https://discourse.julialang.org/u/DanDoe)\
**Post date:** [June 5, 2022, 2:21pm UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/13 "2022-06-05T14:21:08Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 5, 2022, 11:19pm UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/14 "2022-06-05T23:19:42Z")

</div>

Maybe something like this?

```julia
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.

---

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [June 6, 2022, 8:39am UTC](https://discourse.julialang.org/t/linear-optimization-with-quadratic-constraints-efficient-arbitrary-precision-solver/82105/15 "2022-06-06T08:39:17Z")

</div>

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?