# Linearize or approximate a square root constraint

I encounter a nonlinear constraint that contains the square root of a sum of integer variables. Of course one could use nonlinear solvers and techniques; but I like linear programming. Are there any standard results on linearizing or approximating a square root of the sum of integer variables?

For example, the constraints look like this:

$$\sqrt{\sum_{i \in \mathcal{I}} a_{ij}x_{ij} } \leq \theta_j, \quad \quad \forall j \in \mathcal{J}$$

where $$x_{ij} \in \{0,1\}$$ are binary variables, $$\theta_j \in \mathbb{R}$$ are continuous variables, and $$a_{ij} \geq 0$$ are parameters. $$\mathcal{I}$$ and $$\mathcal{J}$$ are any given sets of polynomial size.

Of course, this constraint is part of a larger MIP, but as I am curious to general methods and results regarding this constraint I believe it not to be of interest to post it here.

• Could you add some additional information about the constraint (maybe a simplified version)? Is the constraint convex? – Kevin Dalmeijer Jul 21 at 9:25

This can be handled as an MISOCP, Mixed-Integer Second Order Cone problem. The leading commercial MILP solvers can also handle MISOCP.

Specifically, due to $$x_{ij}$$ being binary, $$x_{ij}^2 = x_{ij}$$. Therefore, the left-hand side is the two-norm of the vector over $$i \in I$$ having elements $$\sqrt{a_{ij}} x_{ij}$$.

I don't know whether this is the best way to handle this constraint, but it is a way, and it is "exact".

• I agree with this answer. To add to this: if you really want to solve an LP rather than an SOCP, you can approximate the SOCP constraint with a polynomial number of LP constraints, as discussed here. AFAIK it's generally understood that solving the problem as one SOCP is faster. – Ryan Cory-Wright Jul 21 at 21:22
• I like this answer a lot. To add: for a good reference on modeling SOCP, see this paper. – Kevin Dalmeijer Jul 21 at 22:06
• @Daniel Duque If $x \in \{0,1\}$ then you can simply replace $x$ by $x^2$. This is due to the fact that $0^2 = 0$ and $1^2 = 1$ and thus $x^2 = x$ for all permitted values. – Kevin Dalmeijer Jul 22 at 9:26
• I see now, thanks. – Daniel Duque Jul 22 at 12:04
• @Alexandre Frias "Convex Optimization" by Boyd and Vandenberghe web.stanford.edu/~boyd/cvxbook . Also, the MOSEK Modeling Cookbook (chapter 3) docs.mosek.com/MOSEKModelingCookbook-letter.pdf – Mark L. Stone Oct 1 at 20:14

Please also have a look at the very similar question in math.stackexchange. As @Mark L. Stone mentioned in his answer, all you need is a second-order cone model to solve your problem.

To linearize that constraint as it is can be hard since it is non-convex. Assuming you still want to do that, you would need to introduce binary variables that allow you characterize the function.

Focusing on a single $$j$$, let first define $$w_j=\sum\limits_{I\in I}a_{i,j} x_{i,j}$$, with $$w_j\geq 0$$ and assume you have a bound on such that $$w_j\leq UB_j$$. Now let $$n$$ be the number of pieces (linear inequalities) you want to use to describe $$\sqrt{w_j}$$, and for each piece, let $$m_{k,j}$$ and $$b_{k,j}$$ be the slope and intercept of the $$k$$th piece of the $$j$$th constraint for $$k=1,\ldots,n$$, which are tangent lines of $$\theta_j=\sqrt{w_j}$$ at (finite) points $$w_{k,j}\in[0,UB_j]$$ (these are the breakpoints in the $$w_j$$ space), $$k=1,\ldots,n+1$$. Since the constraint are not convex, only one piece can be "on" in an optimal solution, hence, let $$\lambda_{k,j}\in\{0,1\}$$ be a binary variable that is one if the piece is "on" for constraint $$j\in J$$, zero otherwise. Putting all together,

• Choose only one piece for crt $$j$$: $$\sum\limits_{k=1}^n{\lambda_{k,j}}=1 \quad\forall j\in J$$

• $$w_j$$ need to be in the right interval if you choose piece $$k$$ $$-M(1-\lambda_{k,j}) + w_{k,j}\le w_j \le w_{k+1,j} + M(1-\lambda_{k,j}) \quad \forall j \in J,\,k=1,\ldots,n$$

• Definition of $$w_j$$: $$w_j = \sum\limits_{I\in I}a_{i,j} x_{i,j} \quad\forall j \in J$$

• This is the linearized constraint, where $$\theta_j$$ is greater or equal to the piece that is selected: $$\theta_j\ge m_{k,j} w_j + b_{k,j} - M(1-\lambda_{k,j}) \quad\forall j\in J,\, k=1,\ldots,n$$

As a side note, you have to choose the breakpoints upfront. A plot of $$\theta_j\ge \sqrt{w_j}$$ (for a single $$j$$, this a 2D-plot) can help to clarify the linearization.

If your constraints are convex (e.g., the inequality is $$\ge$$ or you treat it as an SOCP as described in the answer above), then you could implement Kelley's cutting-plane1 method which is an outer approximation method. Those cuts are not cuts in the integer programming sense, so don't add them as cuts. Rather, in B&B add them as lazy constraints. Alternatively, if the MIP is easy to solve, generate a single (Kelley's) cut at a time an re-optimize.

Reference

[1] Kelley, J. E., Jr. (1960). The Cutting-Plane Method For Solving Convex Programs. Journal of the Society for Industrial and Applied Mathematics. 8(4):703-712.

• Is all this really needed? The square root function is concave, so you could just use a series of cutting planes: linearisation of square root at $\sum\limits_i a_i x_i = \text{sth}\le\theta_j$; with enough cutting planes, you are drawing an approximation that is good enough. (Or you can use an iterative process to add new cutting planes, to reduce the total number of constraints.) I guess that's what pubsonline.informs.org/doi/abs/10.1287/moor.26.2.193.10561 does. – dourouc05 Jul 22 at 1:59
• It is necessary due to the sense of the inequality. In the way it’s written in the question, that constraint defines a non-convex feasible region. If the inequality is the other way around, then it would be much easier as it would be a convex feasible region (ignoring the integrality of x). – Daniel Duque Jul 22 at 2:09

You can manipulate this inequality as follows

$$\sqrt{\sum_{i \in \mathcal{I}} a_{ij}x_{ij} } \leq \theta_j, \quad \quad \forall j \in \mathcal{J}$$

$$\sum_{i \in \mathcal{I}} a_{ij}x_{ij} \leq \theta_j^2, \quad \quad \forall j \in \mathcal{J}$$

Now, you need to linearize $$\theta_j^2$$ using McCormick Envelopes. To do this, assume $$-M_j\leq \theta_j \leq M_j$$ and consider $$w_j=\theta_j^2$$

\begin{align} 0\leq (\theta_j + M_j)(\theta_j + M_j) & \implies & -w_j - 2M_j\theta \leq M_j^2\\ 0\leq (M_j - \theta_j)(M_j - \theta_j) & \implies & -w_j + 2M_j\theta \leq M_j^2\\ 0\leq (\theta_j + M_j)(M_j - \theta_j) & \implies & w_j \leq M_j^2\\ \end{align}

The final set of constraints is

\begin{align} \sum_{i \in \mathcal{I}} a_{ij}x_{ij} \leq w_j, \quad \quad \forall j \in \mathcal{J}\\ -w_j - 2M_j\theta_j \leq M_j^2, \quad \quad \forall j \in \mathcal{J}\\ -w_j + 2M_j\theta_j \leq M_j^2, \quad \quad \forall j \in \mathcal{J}\\ 0 \leq w_j \leq M_j^2, \quad \quad \forall j \in \mathcal{J}\\ -M_j \leq \theta_j \leq M_j, \quad \quad \forall j \in \mathcal{J}\\ \end{align}