Polynomial Optimization:
Enhancing RLT relaxations with Conic Constraints

The large volume of theoretical and computational research on conic optimization has led to important advances over the last few years in the efficiency and robustness of the associated algorithmic procedures to solve them. Leading state-of-the-art mixed integer linear programming (MILP) solvers such as Gurobi (Gurobi Optimization, 2022), CPLEX (IBM Corp., 2022), and Xpress (FICO, 2022) have recently added functionalities that allow to efficiently solve second-order cone programming (SOCP) problems and, further, Mosek (MOSEK ApS, 2022, Andersen and Andersen, 2000) has positioned itself as a reliable solver for general semidefinite programming (SDP) problems.

Importantly, despite the convex nature of conic optimization problems, they are also proving to be a powerful tool in the design of branch-and-bound algorithms for general non-convex mixed integer nonlinear programming (MINLP) problems and, specially, polynomial optimization problems. The latter are very general and encompass many problems arising in operations research: integer programming, linear programming, mixed integer programming, and quadratic optimization, to name a few. The developments of the last several decades have shown that conic optimization is a central tool in addressing non-convexities. These advances have been more prominent in the case of polynomial optimization problems, and even more so in the particular case of (quadratically-constrained) quadratic optimization problems, where a variety of convex relaxations have been thoroughly studied (Shor, 1987, Ghaddar et al., 2011b , Burer and Ye, 2020, Bonami et al., 2019, Elloumi and Lambert, 2019). These conic-based relaxations include semidefinite programming and second-order cone programming as their powerhouses. They often provide fast and guaranteed approaches for computing bounds on the global value of the non-convex optimization problem at hand. For instance, Lasserre, (2001) introduced semidefinite relaxations corresponding to liftings of the polynomial programs into higher dimensions. The construction is motivated by results related to representations of non-negative polynomials as sum-of-squares (Parrilo, 2003) and the dual theory of moments. This led to the development of systematic approaches for solving polynomial optimization problems to global optimality, the main limitation of these approaches currently being that they are computationally prohibitive in general.

Nowadays, state-of-the-art solvers tackle MINLP problems and, in particular, polynomial optimization ones through different relaxations to be solved at the nodes of the branch-and-bound tree. A crucial direction of current research focuses on integrating tighter and tighter relaxations while preserving reasonable computational properties, with relaxations that build upon different families of conic constraints becoming increasingly important. One relaxation for polynomial optimization problems is based on the Reformulation-Linearization Technique (RLT) and the most common approach to strengthen this relaxation comes from the identification of linear SDP-based cuts that tighten the RLT relaxation as the algorithm progresses (Sherali et al., 2012, Baltean-Lugojan et al., 2019, González-Rodríguez et al., 2020). In this paper we are interested in the alternative approach of directly adding SOCP and SDP constraints to the relaxations as in Burer and Vandenbussche, (2008) and Buchheim and Wiegele, (2013). Probably the main reason why this approach has received less attention so far is that the resulting algorithms need to rely on SDP solvers which, until recently, were probably not reliable enough and too expensive computationally in some instances. One of the primary goals of this paper is to show that the situation has changed, and that general branch-and-bound schemes based on the solution of SOCPs and SDPs at each and every node of the branch-and-bound tree can be very competitive and even superior to previous approaches.

In order to achieve the above goal, in this paper we introduce and compare linear SDP-based constraints as well as SOCP and SDP constraints. Importantly, their definition ensures that they are efficient in the sense of preserving the size and sparsity of the original RLT-based relaxation for polynomial optimization problems introduced in Sherali and Tuncbilek, (1992) and further refined in Dalkiran and Sherali, (2013). We consider a total of nine different versions of constraints to be added: one based on linear SDP-based cuts, four based on SOCP constraints, and four based on SDP constraints. These conic constraints are then integrated into the polynomial optimization solver RAPOSa (González-Rodríguez et al., 2020), whose core is an RLT-based branch-and-bound algorithm. As a second step, we then develop thorough computational studies on well established benchmarks including randomly generated instances from Dalkiran and Sherali, (2016), MINLPLib instances (Bussieck et al., 2003), and QPLIB instances (Furini et al., 2018), which provide a wide variety of classes of polynomial optimization problems. One of the main findings is that, in around $50\%$ of the instances, the best performance is achieved by one of the versions explicitly incorporating SOCP and SDP conic constraints. The remaining $50\%$ is split quite evenly between the baseline RLT and the version that incorporates SDP-based linear cuts. Importantly, the analysis also allowed to identify particular classes of problems where one specific family of SOCP/SDP conic constraints is consistently superior to the linear versions.

To the best of our knowledge, our contribution represents one of very few implementations of branch-and-bound schemes with conic relaxations for broad classes of problems, with the added value of the generality of the resulting scheme, since it can be applied to any given polynomial optimization problem. The most related approaches are Burer and Vandenbussche, (2008) for non-convex quadratic problems with linear constraints and Buchheim and Wiegele, (2013) for unconstrained mixed-integer quadratic problems. In Burer and Vandenbussche, (2008), the authors present a computational analysis in which they compare the relative performance of different SDP relaxations, with the main highlight being that “only a small number of nodes are required” to fully solve the problems at hand. Buchheim and Wiegele, (2013) introduce a new branch-and-bound algorithm, Q-MIST (Quadratic Mixed-Integer Semidefinite programming Technique), and develop a thorough computational analysis in which they show that, for a wide variety of families of instances within the scope of their algorithm, Q-MIST notably outperforms Couenne (Belotti et al., 2009). It is worth noting that, if looking at specific classes of optimization problems, then one can find additional successful implementations of branch-and-bound algorithms with the inclusion of tailor-made conic constraints, such as Ghaddar and Jabr, (2019) for optimal power flow problems and combinatorial optimization problems, Krislock et al., (2017) and Rendl et al., (2010) for maximum cut problems, Piccialli et al., (2022) for minimum sum-of-squares clustering, and Ghaddar et al., 2011a for maximum $k$-cut problems.

Last, but not least, in our computational experiments we also observe that there is a lot of variability in the best performing version of conic constraints for the different instances, with each version beating the rest for a non negligible number of instances. This observation motivates the last contribution of this paper, in which we exploit this variability by learning to choose the best version among our portfolio of different constraints. Building upon the framework in Ghaddar et al., (2022), we show that the resulting machine learning version significantly outperforms each and every one of the underlying versions. This last contribution naturally fits into the rapidly emerging strand of research on “learning to optimize”, whose advances are nicely presented in the survey papers Lodi and Zarpellon, (2017) and Bengio et al., (2021). The closest approach to this last part of our contribution is Baltean-Lugojan et al., (2019), in which deep neural networks are used to rank SDP-based cuts for quadratic problems. Then, only the top scoring cuts are added, aiming to obtain a good balance between the tightness and the complexity of the relaxations. A common feature of the SDP-based cuts used in Baltean-Lugojan et al., (2019) and those in this paper is that they are designed with a big emphasis in sparsity considerations (number of nonzeros in the constraints), with the goal of obtaining computationally efficient relaxations. From the point of view of the learning process the approaches are quite different, since they focus on one type of constraints (the SDP-based cuts) and they want to learn the best SDP-based cuts to add at a given node, whereas we want learn to choose for any given instance which conic constraints to include in the relaxations.

The contribution of this paper can be summarized as follows. First, we define and show the potential of different SOCP and SDP strengthenings of the classic RLT relaxations for general polynomial optimization problems in a branch-and-bound scheme. Second, we design a machine learning approach to learn the best strenghtening to use on a given instance and obtain promising results for the resulting algorithm.

The remainder of this paper is organized as follows. In Section 2 we present a brief overview of the classic RLT scheme. In Section 3 we describe the different families of conic constraints that will be integrated within the baseline RLT implementation. In Section 4 we present a first series of computational results. Then, in Section 5 we show how the conic constraints can be further exploited within a machine learning framework. Finally, we conclude in Section 6 and discuss future research directions.

The Reformulation-Linearization Technique was originally developed in Sherali and Tuncbilek, (1992). It was designed to find global optima in polynomial optimization problems of the following form:

where $N=\{1,\dots,n\}$ denotes the set of variables, each $\phi_{r}(\mathbf{x})$ is a polynomial of degree $\delta_{r}\in\mathbb{N}$ and $\Omega=\{\mathbf{x}\in\mathbb{R}^{n}:0\leq l_{j}\leq x_{j}\leq u_{j}<\infty,\,\forall j\in N\}\subset\mathbb{R}^{n}$ is a hyperrectangle containing the feasible region. Then, $\delta=\max_{r\in\{0,\ldots,R\}}\delta_{r}$ is the degree of the problem and $(N,\delta)$ represents all possible monomials of degree $\delta$.

The Reformulation-Linearization Technique consists of a branch and bound algorithm based on solving linear relaxations of the polynomial problem (PO). These linear relaxations are built by working on a lifted space, where each monomial of the original problem is replaced with a corresponding RLT variable. For example, associated to monomials of the form $x_{1}x_{2}x_{4}$ and $x_{1}^{2}x_{3}^{2}$ one would define the RLT variables $X_{124}$ and $X_{1133}$, respectively. More generally, RLT variables are defined as

where $J$ is a multiset containing the information about the multiplicity of each variable in the underlying monomial. Then, at each node of the branch-and-bound tree, one would solve the corresponding linear relaxation. Whenever we get a solution of a linear relaxation in which the identities in (1) hold, we get a feasible solution of (PO). Otherwise, the violations of these identities are used to choose the branching variable.

In order to get tighter relaxations and ensure convergence, new constraints, called bound-factor constraints, must be added. They are of the following form:

Thus, for each pair of multisets $J_{1}$ and $J_{2}$ such that $J_{1}\cup J_{2}\subset(N,\delta)$ and $|J_{1}\cup J_{2}|=\delta$, the corresponding bound-factor constraint is added to the linear relaxation.

Dalkiran and Sherali, (2013) show that it is not necessary to add all bound-factor constraints to the linear relaxations, since certain subsets of them are enough to ensure the convergence of the algorithm. More precisely, they proved that convergence to a global optimum only requires the inclusion in the linear relaxation of those bound-factor constraints where $J_{1}\cup J_{2}$ is a monomial that appears in (PO), regardless of its degree. Further, they also showed that convergence is also preserved if, whenever the bound-factor constraints associated to a monomial $J$ are present, all bound-factor constraints associated to monomials $J^{\prime}\subset J$ are removed. Motivated by these results, $J$-sets are defined as those monomials of degree greater than one present in (PO) which, moreover, are not included in any other monomial (multiset inclusion). Consider the polynomial programming problem

The monomials with degree greater than one are $\{1,1\}$ ,$\{2,2\}$, $\{1,2,3\}$, $\{1,2\}$, and $\{1,4\}$. Since $\{1,2\}$ is included in $\{1,2,3\}$, it is removed. Therefore, the $J$-sets are $\{1,1\}$, $\{2,2\}$, $\{1,2,3\}$, and $\{1,4\}$.

The analysis developed in this paper builds upon the above theoretical results from Dalkiran and Sherali, (2013) and, therefore, only the bound-factor constraints associated to $J$-sets are incorporated to the linear relaxations of (PO). Such relaxations are less tight but, on the other hand, they are smaller in size and, hence, faster to solve. Note that the use of $J$-sets reduces not only the number of constraints but, more importantly, also the number of RLT variables. As it can be seen in Dalkiran and Sherali, (2013) and González-Rodríguez et al., (2020), the performance of the RLT algorithm is clearly superior when the $J$-set approach is followed.

As already discussed earlier, the use of semidefinite programming to improve the performance of branch-and-bound schemes is not new, and Baltean-Lugojan et al., (2019) provides a thorough and up-to-date review of the field. Typically, the goal is to rely on semidefinite programming to tighten the relaxations of the original non-convex optimization problems targeted by a branch-and-bound algorithm. We start this section by reviewing the main ingredient of such methods and then present various families of SDP-driven constraints that can be incorporated into the RLT relaxations in an efficient way. In particular, they should preserve the underlying dimensionality and sparsity of the problem, which is crucial for these approaches to be competitive.

The main ingredient that has to be specified is the matrix or matrices on which positive semidefiniteness is to be imposed. To each (multi-)set of variables $\{x_{j}\}_{j\in J}$, with $J\subset(N,\delta)$, one can associate a vector $\omega=(x_{j})_{j\in J}$ (resulting from the concatenation of all the variables in $J$, including repetitions). To any such vector one can associate matrix $M=\omega^{T}\omega$, which is trivially positive semidefinite. Now, let $M_{L}=[\omega^{T}\omega]_{L}$ be the matrix obtained when each monomial in $M$ is replaced by the corresponding RLT variable in the lifted space. The constraint $M_{L}\succcurlyeq 0$ is a valid cut because it never removes feasible solutions of (PO) and, hence, it does not compromise convergence of the RLT algorithm to a global optimum. In practice, vectors of the form $(1,(x_{j})_{j\in J})$ are often preferred, since they result in $M_{L}$ matrices containing also variables in the original space and not only RLT variables, leading to tighter relaxations. In Sherali et al., (2012), for instance, the authors discuss different ways of defining $w$ for a given $J$. Specifically, they mainly work with $\omega^{1}=(x_{j})_{j\in J}$, $\omega^{2}=(1,(x_{j})_{j\in J})$, and $\omega^{3}=(1,x_{1},x_{2},\ldots)$, where $\omega^{3}$ is defined by concatenating also monomials of degree greater than one, while ensuring that no monomial in the resulting matrix $M$ has degree larger than $\delta$, the degree of (PO).

We now move to the definition of the specific SDP-driven constraints for the RLT algorithm that constitute the subject of study in this paper.

In this section we want to exploit the wide variability in the performance of the different RLT versions shown above to try to learn in advance which one should be chosen for a given instance. The goal is to design a machine learning procedure that can be trained using the different features of the instances and then choose the most promising RLT version when confronted with a new instance. The performance of the hypothetical “Optimal version” in Table 1 represents an upper bound on the improvement that can be attained by such a machine learning version.

We follow the framework in Ghaddar et al., (2022), where the authors use learning techniques to improve the performance of the RLT-based solver RAPOSa by learning to choose between different branching rules. The improvements reported there are substantial, with the machine learning version delivering improvements of up to $25\%$ with respect to the best original branching rule. Table 2 presents the full list of the input variables (features).²²2VIG and CMIG stand for two graphs that can be associated to any given polynomial optimization problem: variables intersection graph and constraints-monomials intersection graph, and whose precise definitions is given in Ghaddar et al., (2022). They capture diverse characteristics of each instance and are a key ingredient of the machine learning framework, which consists of predicting the performance ($\operatorname{npace}^{LB}$) of each branching rule on a new instance based on a regression analysis of its performance on the training instances. Then, the rule with a highest predicted performance for the given instance is chosen.

For the learning we rely on quantile regression models since as argued in Ghaddar et al., (2022), the presence of outliers and of an asymmetric behavior of $\operatorname{npace}^{LB}$ (negative skewness) makes quantile regression more suitable than conventional regression models (based on the conditional mean). More precisely, we use quantile regression forests as the core tool for our analysis. Random forests were introduced in Breiman, (2001) as ensemble methods aggregate the information on several individual decision trees to make a single prediction. Quantile regression forests, introduced in Meinshausen, (2006), are a generalization of random forests that compute an estimation of the conditional distribution of the response variable by taking into account all the observations in every leaf of every tree and not just their average. As discussed in Ghaddar et al., (2022), one advantage of random forests is that the reporting of the results can be done for the complete data set in terms of Out-Of-Bag predictions; there is no need to do the reporting with respect to the underlying training and test sets. The statistical analysis is conducted in programming language R (R Core Team, 2021), using library ranger (Wright and Ziegler, 2017). Finally, the learning is conducted jointly on all sets of instances.

Table 3 shows the remarkable improvement obtained with the machine learning (ML) version of RLT that chooses, for each given instance, the most promising version of the 10-version portfolio. Considering all instances, the ML version improves $32.9\%$ with respect to RLT, being the upper bound for learning $50.5\%$. In instances with high densities from DS test set, it improves a remarkable $69.5\%$ out of the optimal $74.4\%$, and in “water”-related instances in MINLPLib it improves $56.7\%$ out of $92.5\%$. It is worth noting that the size of these improvements is comparable, and even superior, to those obtained in Ghaddar et al., (2022) when this learning scheme was introduced to learn to choose between branching rules, which shows the robustness of the proposed approach.

Figure 3 represents, side by side, how often each of the ten RLT versions is selected by the ML version and by the optimal one. We can see that the behavior of the former mimics quite well that of the latter in the three sets of instances. Given that the learning is conducted jointly on the whole set of instances, the fact that the ML version adapts to the instances in DS, MINLPLib, and QPLib, is reassuring about the quality of the learning process. In particular, the SOCP versions with Mosek are primarily chosen for DS test set and the SDP versions are mainly chosen for QPLIB test set, the one in which they are optimal more often. In Figure 4 we further explore the behavior for the high density problems in DS and for “water”-related instances in MINLPLib. We see that, again, the ML version mimics the patterns of the optimal version. The dominant version in DS, SOCP^M, is chosen in almost $75\%$ of the instances. Similarly, the three SDP dominant versions in “water”-related instances are chosen almost $50\%$ of the time. The dominant version in DS, SOCP^M, is chosen in almost $75\%$ of the instances.

Figure 5 presents boxplots summarizing the performance according to $\operatorname{npace}^{LB}$ for all instances and for each individual test set. Recall that, by definition, values close to 1 mean that the corresponding version is almost the best one, whereas values close to 0 mean that its pace is much worse than the best one. We can see that, although RLT and SDP cuts versions are, on aggregate, the best ones in all three test sets, they are significantly outperformed by the ML version in the three of them. The improvement is particularly noticeable in DS, where the learning criterion is, by far, better than all underlying versions. This observation is further reinforced by the performance profiles (Dolan and Moré, (2002)) represented in Figure 6. Again, the ML version clearly outperforms all others, specially in DS instances.

The main contribution of this paper is to show that the solution to global optimality via branch-and-bound schemes of non-convex optimization problems and, in particular, polynomial optimization ones, can benefit from tightening the underlying relaxations with conic constraints. We explore different families of such constraints, building upon either second-order cones or positive semidefiniteness. We also show that the potential of these conic constraints can be successfully exploited by embedding it into a learning framework. The main goal is to predict which is the most promising type of constraints to add to the RLT relaxation at each node of the underlying branch-and-bound algorithm when confronted with a new instance. The results in Section 4 show that the versions with SOCP/SDP conic constraints deliver consistently good results for instances in specific subclasses of problems: high density problems in DS and of those based on positive semidefinite constraints for “water”-related instances in MINLPLib.

As a future step, one may wonder to what extent one might get an even superior performance if the learning analysis was further specialized for the current setting: for example including features capturing some “conic” characteristics of the polynomial optimization problems and fine-tuning the regression techniques. Furthermore, an important direction for future research is to improve the understanding on the structure of these problems and the specificities that lead to the superior performance of SOCP and SDP constraints, respectively. Another direction is to investigate the number of SOCP and SDP constraints added at different nodes of the branch-and-bound tree. Additionally, we aim to extend this framework for various relaxations of polynomial optimization problems besides RLT.

This research has been funded by FEDER and the Spanish Ministry of Science and Technology through projects MTM2014-60191-JIN, MTM2017-87197-C3, and PID2021-124030NB-C32. Brais González-Rodríguez acknowledges support from the Spanish Ministry of Education through FPU grant 17/02643. Raúl Alvite-Pazó and Samuel Alvite-Pazó acknowledge support from CITMAga through proyect ITMATI-R-7-JGD. Bissan Ghaddar’s research is supported by Natural Sciences and Engineering Research Council of Canada Discovery Grant 2017-04185 and by the David G. Burgoyne Faculty Fellowship.