Solving unconstrained 0-1 polynomial programs through quadratic convex reformulation

TL;DR

This paper introduces PQCR, a three-phase method for solving unconstrained binary polynomial optimization problems by reformulating them into quadratic convex programs, leveraging semidefinite relaxations for optimal parameter selection.

Contribution

The paper presents a novel three-phase approach that reformulates binary polynomial problems into quadratic convex programs using semidefinite relaxations for parameter optimization.

Findings

PQCR effectively solves instances of image restoration and autocorrelation problems.

The method outperforms existing convexification techniques and the Baron solver in computational tests.

Abstract

We propose a solution approach for the problem (P) of minimizing an unconstrained binary polynomial optimization problem. We call this method PQCR (Polynomial Quadratic Convex Reformulation). The resolution is based on a 3-phase method. The first phase consists in reformulating (P) into a quadratic program (QP). For this, we recursively reduce the degree of (P) to two, by use of the standard substitution of the product of two variables by a new one. We then obtain a linearly constrained binary program. In the second phase, we rewrite the quadratic objective function into an equivalent and parametrized quadratic function using the equality x 2 i = x i and new valid quadratic equalities. Then, we focus on finding the best parameters to get a quadratic convex program which continuous relaxation's optimal value is maximized. For this, we build a semidefinite relaxation (SDP) of (QP). Then, we prove that the standard linearization inequalities, used for the quadratization step, are redundant in (SDP) in presence of the new quadratic equalities. Next, we deduce our optimal parameters from the dual optimal solution of (SDP). The third phase consists in solving (QP *), the optimal reformulated problem, with a standard solver. In particular, at each node of the branch-and-bound, the solver computes the optimal value of a continuous quadratic convex program. We present computational results on instances of the image restoration problem and of the low autocorrelation binary sequence problem. We compare PQCR with other convexification methods, and with the general solver Baron 17.4.1 [39]. We observe that most of the considered instances can be solved with our approach combined with the use of Cplex [24].

Tables

Table 1. Table 1: Comparison of the maximum time on 4 solution methods for the vision instances - time limit one hour

Instance				PQCR			Q+QCR			Q+Cplex		Baron
Name	n	m	N	Gap	tSdp	tTotal	Gap	tSdp	tTotal	Gap	tTotal	Gap	tTotal
v.10.10 1	100	668	352	0,59	66	68	396	7	(250 %)	1113	2	1098	15
v.10.10 2	100	668	352	0,28	64	66	536	8	(343 %)	1549	2	1529	10
v.10.10 3	100	668	352	0,05	65	67	973	8	(573 %)	3375	1	3332	6
v.10.10 4	100	668	352	0,12	63	65	957	8	(561 %)	3377	1	3334	6
v.10.10 5	100	668	352	0,13	65	66	1006	8	(585 %)	3568	1	3523	5
v.10.10 6	100	668	352	0,11	73	74	359	9	(229 %)	984	2	972	11
v.10.10 7	100	668	352	0,02	64	65	305	8	(194 %)	829	2	817	14
v.10.10 8	100	668	352	1,37	64	66	1376	8	(804 %)	4765	1	4705	7
v.10.10 9	100	668	352	3,02	65	67	1749	8	(1026 %)	6187	1	6110	4
v.10.10 10	100	668	352	3,64	66	68	1879	8	(1075 %)	6843	1	6757	4
v.10.10 11	100	668	352	0,36	70	72	489	8	(316 %)	1388	2	1370	35
v.10.10 12	100	668	352	0,20	70	72	361	9	(232 %)	997	2	984	23
v.10.10 13	100	668	352	0,00	60	61	709	8	(392 %)	2654	1	2620	2
v.10.10 14	100	668	352	0,00	60	61	546	8	(297 %)	2027	1	2001	2
v.10.10 15	100	668	352	0,00	118	119	541	8	(285 %)	2048	1	2022	1
v.10.15 1	150	1033	542	0,31	290	294	447	24	(351 %)	1245	5	1234	80
v.10.15 2	150	1033	542	0,00	285	287	367	24	(287 %)	999	5	990	36
v.10.15 3	150	1033	542	0,05	280	283	1027	27	(772 %)	3549	3	3520	6
v.10.15 4	150	1033	542	0,33	276	280	845	27	(640 %)	2840	3	2817	7
v.10.15 5	150	1033	542	0,07	269	271	799	27	(595 %)	2808	3	2785	4
v.10.15 6	150	1033	542	0,55	297	302	462	25	(366 %)	1277	5	1266	47
v.10.15 7	150	1033	542	0,08	288	291	360	26	(283 %)	981	5	972	38
v.10.15 8	150	1033	542	0,79	281	284	1792	26	(1356 %)	6202	3	6152	19
v.10.15 9	150	1033	542	1,80	283	286	1525	26	(1160 %)	5209	3	5167	10
v.10.15 10	150	1033	542	1,38	275	279	1510	25	(1124 %)	5500	3	5456	7
v.10.15 11	150	1033	542	0,10	283	286	391	25	(305 %)	1102	5	1092	41
v.10.15 12	150	1033	542	0,60	275	279	453	25	(355 %)	1269	5	1258	125
v.10.15 13	150	1033	542	0,00	254	256	634	27	(469 %)	2254	3	2236	4
v.10.15 14	150	1033	542	0,04	269	273	731	27	(547 %)	2590	3	2569	2
v.10.15 15	150	1033	542	0,00	258	259	576	28	(423 %)	2183	2	2165	2
v.15.15 1	225	1598	827	0,12	1234	1244	365	70	(320 %)	998	9	993	(95 %)
v.15.15 2	225	1598	827	0,43	1251	1265	482	83	(421 %)	1350	9	1343	(138 %)
v.15.15 3	225	1598	827	0,10	1167	1175	678	65	(582 %)	2326	5	2313	(83 %)
v.15.15 4	225	1598	827	0,04	1251	1256	877	65	(753 %)	2996	5	2980	(127 %)
v.15.15 5	225	1598	827	0,03	1167	1174	641	67	(546 %)	2252	5	2240	(76 %)
v.15.15 6	225	1598	827	0,28	1238	1249	403	64	(353 %)	1104	10	1098	(107 %)
v.15.15 7	225	1598	827	0,36	1237	1246	525	67	(463 %)	1455	10	1447	(144 %)
v.15.15 8	225	1598	827	0,29	1197	1205	1148	73	(979 %)	4104	5	4082	(137 %)
v.15.15 9	225	1598	827	0,27	1170	1176	1542	66	(1315 %)	5570	5	5541	(171 %)
v.15.15 10	225	1598	827	0,31	1173	1179	1194	67	(1020 %)	4380	5	4357	1154
v.15.15 11	225	1598	827	0,27	1224	1230	529	69	(462 %)	1528	8	1520	(144 %)
v.15.15 12	225	1598	827	0,25	1225	1235	461	68	(414 %)	1273	12	1266	(133 %)
v.15.15 13	225	1598	827	0,00	1124	1128	651	63	(551 %)	2398	4	2385	1239
v.15.15 14	225	1598	827	0,02	1171	1177	651	63	(553 %)	2359	5	2346	(79 %)
v.15.15 15	225	1598	827	0,00	1100	1103	609	65	(513 %)	2320	4	2308	263

Table 2. Table 2: Results of PQCR  and Baron  for the 45 instances of the LABS problem. Time limit 5 hours. - means that the time limit of 3h on the SDP phase is reached.

Instance			PQCR					Baron 17.4.1
Name	n	m	N	Gap	tSdp	tTotal	Nodes	Gap	tTotal	Nodes
b.20.03	20	38	20	0	1	2	0	100	1	1
b.20.05	20	207	65	23	22	23	5886	1838	2	1
b.20.10	20	833	124	8	837	846	24183	2918	125	7
b.20.15	20	1494	164	5	1228	1242	9130	3202	728	9
b.25.03	25	48	25	0	1	2	0	100	0	1
b.25.06	25	407	105	17	461	469	163903	2307	65	27
b.25.13	25	1782	206	4	1552	1603	76828	3109	3750	75
b.25.19	25	3040	265	4	-	13433	224550	3356	14399	129
b.25.25	25	3677	289	5	-	13395	167423	3405	(12 %)	100
b.30.04	30	223	82	23	58	78	134635	1347	7	7
b.30.08	30	926	174	10	1940	2040	752765	2696	2778	237
b.30.15	30	2944	296	5	-	13525	438278	3221	(21 %)	103
b.30.23	30	5376	390	11	5953	6865	9337391	3450	(135 %)	8
b.30.30	30	6412	422	4	8500	15352	452460	3470	(161 %)	5
b.35.04	35	263	97	19	135	167	156085	1350	32	13
b.35.09	35	1381	234	10	2245	4630	8163651	2826	(29 %)	354
b.35.18	35	5002	419	644	-	(12 %)	4899872	3356	(133 %)	10
b.35.26	35	8347	530	30	-	(5 %)	5006407	3508	(229 %)	3
b.35.35	35	10252	579	12	-	(11 %)	134426	3499	(214 %)	3
b.40.05	40	447	145	25	430	1630	23459121	1856	3674	1021
b.40.10	40	2053	304	9	-	(4 %)	25480163	2953	(54 %)	147
b.40.20	40	7243	544	9	-	(4 %)	9783350	3405	(203 %)	3
b.40.30	40	12690	702	360	-	(25 %)	281134	3561	(274 %)	1
b.40.40	40	15384	762	62	-	(44 %)	57534	3536	(464 %)	1
b.45.05	45	507	165	24	1384	(4 %)	84159279	1854	16609	4727
b.45.11	45	2813	382	9	-	(2 %)	25114985	3018	(132 %)	33
b.45.23	45	10776	706	21	-	(16 %)	1225234	3470	(242 %)	2
b.45.34	45	18348	898	137	-	(105 %)	38513	3604	(375 %)	1
b.45.45	45	21993	969	187	-	(153 %)	25964	3559	(624 %)	1
b.50.06	50	882	230	19	1230	(9 %)	49490829	2321	(35 %)	1225
b.50.13	50	4457	506	8	-	(5 %)	12039566	3131	(192 %)	7
b.50.25	50	14412	866	676	-	(247 %)	684010	3511	(280 %)	1
b.50.38	50	25446	1118	242	-	(163 %)	309289	3646	(505 %)	1
b.50.50	50	30271	1202	360	-	(305 %)	49507	3541	(729 %)	1
b.55.06	55	977	255	21	-	(11 %)	23603952	2323	(54 %)	6
b.55.14	55	5790	607	11	-	(7 %)	7829649	3186	(373 %)	6
b.55.28	55	19897	1069	174	-	(106 %)	580827	3553	(646 %)	2
b.55.41	55	33318	1347	330	-	(244 %)	117912	3654	(639 %)	1
b.55.55	55	40402	1459	547	-	(493 %)	117027	3575	(705 %)	1
b.60.08	60	2036	384	12	-	(9 %)	24800852	2712	(175 %)	1
b.60.15	60	7294	716	16	-	(14 %)	4044387	3236	(404 %)	1
b.60.30	60	25230	1264	256	-	(165 %)	295197	3578	(471 %)	1
b.60.45	60	43689	1614	547	-	(439 %)	26955	704	(671 %)	1
b.60.60	60	52575	1742	784	-	(704 %)	23716	3604	(762 %)	1

Table 3. Table 3: Comparison of the best known solution and best lower bound values of PQCR  and of the minlplib  for the unsolved LABS instances. ∗ ⁣ ∗ ** : solved for the first time, # # \# : best known solution improved, and ∗ * : best known lower bound improved

Instance	PQCR (5h)		minlplib [35]
Name	Solution	$L{B}_f$	Solution	$L{B}_f$
b.25.19∗∗	-14644	-14644	-14644	-16108
b.25.25∗∗	-10664	-10664	-10664	-12494
b.30.15∗∗	-15744	-15744	-15744	-19780
b.30.23∗∗	-30460	-30460	-30420	-72030
b.30.30∗∗	-22888	-22888	-22888	-54014
b.35.09∗∗	-5108	-5108	-5108	-6312
b.35.18∗	-31144	-34964	-31160	-74586
b.35.26#∗	-55288	-57789	-55184	-191466
b.35.35∗	-41052	-45787	-41068	-290424
b.40.10#∗	-8248	-8551	-8240	-14618
b.40.20#∗	-50576	-52465	-50516	-162365
b.40.30#∗	-94872	-118324	-94768	-398617
b.40.40∗	-67528	-98031	-67964	-302028
b.45.05∗	-1068	-1112	-1068	-1145
b.45.11#∗	-12748	-13035	-12740	-30771
b.45.23#∗	-85423	-98984	-85248	-320397
b.45.34∗	-151352	-311627	-152368	-752427
b.45.45∗	-111292	-285811	-112764	-685911
b.50.06∗	-2160	-2363	-2160	-2921
b.50.13#∗	-23791	-24975	-23772	-74768
b.50.25∗	-124572	-433247	-124748	-562446
b.50.38∗	-232344	-611906	-232496	-1318325
b.50.50∗	-162640	-681105	-168216	-1173058
b.55.06∗	-2400	-2659	-2400	-3439
b.55.14#∗	-33272	-35698	-33168	-116748
b.55.28∗	-189896	-392929	-190472	-989145
b.55.41∗	-335388	-1160180	-337388	-2494477
b.55.55∗	-233648	-1434663	-241912	-1947633
b.60.08∗	-6792	-7388	-6792	-13915
b.60.15#∗	-45232	-51467	-44896	-169767
b.60.30∗	-259271	-692721	-261048	-1491016
b.60.45∗	-475504	-2579935	-478528	-3687344
b.60.60∗	-343400	-2816441	-350312	-3021077

Table 4. Table 4: Comparison of bounds and number of variables of PQCR after different quadratizations

		$PC1$		$PC2$		$PC3$		$Quad$
	Opt	N	$L{B}_i$	N	$L{B}_i$	N	$L{B}_i$	N	$L{B}_i$
b.20.05	-416	64	-435	56	-439	40	-436	65	-435
b.20.10	-2936	123	-3052	135	-3115	93	-3068	124	-3051
b.40.10	-8248	303	-8590	315	-8659	262	-8745	304	-8589

Equations

$$(P)$$

$$(QCQP)$$

$$(C_{i_1,i_2}^i)$$

$${\mathcal{F}}_{\mathcal{E}}=\{x\in \{0,1{\}}^N:C_{i_1,i_2}^i\text{ is satisfied }\forall (i,i_1,i_2)\in J\times (I\cup J{)}^2:\text{leavevmode}{\mathcal{E}}_i={\mathcal{E}}_{i_1}\cup {\mathcal{E}}_{i_2}\}.$$

$$(QP)$$

$$(Ex)\Big{\{}\min_{x\in\{0,1\}^{4}}2x_{1}+3x_{2}x_{3}-2x_{2}x_{3}x_{4}-3x_{1}x_{2}x_{3}x_{4}$$

$$(QE{x}_1)$$

$$(QE{x}_2)$$

$$(QE{x}_3)$$

$$(S_{\mathcal{E}})$$

$$x_i{x}_j$$

$$g_{\alpha ,\beta ,\delta ,\lambda }(x)$$

$$(Q{P}_{\alpha ,\beta ,\delta ,\lambda })$$

$$(CP):\max\limits_{\begin{subarray}{c}\alpha\in\mathbb{R}^{N},\beta\in\mathbb{R}^{T_{1}},\delta\in\mathbb{R}^{T_{2}},\lambda\in\mathbb{R}^{T_{3}}\\ Q_{\alpha,\beta,\delta,\lambda}\succeq 0\end{subarray}}\leavevmode\nobreak\ \leavevmode\nobreak\ \Big{\{}\min\limits_{\begin{subarray}{c}x\end{subarray}\in\overline{\mathcal{F}}_{\mathcal{E}}}\leavevmode\nobreak\ \leavevmode\nobreak\ g_{\alpha,\beta,\delta,\lambda}(x)\leavevmode\nobreak\ \Big{\}}$$

$$(SDP)$$

Full text

Abstract

We propose an exact solution approach for the problem $(P)$ of minimizing an unconstrained binary polynomial optimization problem. We call PQCR

(Polynomial Quadratic Convex Reformulation) this three-phase method. The first phase consists in reformulating $(P)$ into a quadratic program $(QP)$. To that end, we recursively reduce the degree of $(P)$ to two, by use of the standard substitution of the product of two variables by a new one. We then obtain a linearly constrained binary quadratic program. In the second phase, we rewrite the objective function of $(QP)$ into an equivalent and parameterized quadratic function using the identity $x_{i}^{2}=x_{i}$ and other valid quadratic equalities that we introduce from the reformulation of phase 1. Then, we focus on finding the best parameters to get a quadratic convex program which continuous relaxation’s optimal value is maximized. For this, we build a new semi-definite relaxation $(SDP)$ of $(QP)$. Then, we prove that the standard linearization inequalities, used for the quadratization step, are redundant in presence of the new quadratic equalities. Next, we deduce our optimal parameters from the dual optimal solution of $(SDP)$. The third phase consists in solving $(QP^{*})$, the optimally reformulated problem, with a standard solver. In particular, at each node of the branch-and-bound, the solver computes the optimal value of a continuous quadratic convex program. We present computational results where we compare PQCR  with other convexification methods, and with the solver Baron [41]. We evaluate our method on instances of the image restoration problem [17] and the low auto-correlation binary sequence problem [7] from minlplib [35]. For this last problem, 33 instances among the 45 were unsolved in minlplib. We solve to optimality 6 of them, and for the 27 others we improve primal and/or dual bounds.

1 Introduction

In this paper, we are interested in solving the unconstrained binary polynomial optimization problem that can be stated as follows:

[TABLE]

where $I=\{1,..,n\}$, $f(x)$ is an $n-$variable polynomial of degree $d$ and $m$ is the number of monomials. For a monomial $p$, $\mathcal{M}_{p}$ is the subset of $I$ containing the indexes of the variables involved in $p$. It follows that $d=\max_{p}|\mathcal{M}_{p}|$.

Unconstrained binary polynomial optimization is a general model that allows to formulate many important problems in optimization. The special case where the polynomial objective function of $(P)$ is a quadratic function (i.e. $d=2$) has been widely studied. In this case, $(P)$ has many applications, including those from financial analysis [31], cluster analysis [39], computer aided design [27] or machine scheduling [40]. Moreover, many graph combinatorial optimization problems such as determining maximum cliques, maximum cuts, maximum vertex packing or maximum independent sets can be formulated as quadratic optimization problems [5, 13, 37]. In the cubic case (i.e. $d=3$), the important class of satisfiability problems known as 3-SAT, can be formulated as $(P)$ [26]. In the case where $d\geq 3$, there also exists many applications including, for example: the construction of binary sequences with low aperiodic correlation [7] that is one of the most challenging problems in signal design theory, or the image restoration problem in computer vision [17].

Problem $(P)$ is $\mathcal{NP}$-hard [20]. Practical difficulties come from the non-convexity of $f$ and the integrality of its variables. During the last decade, several algorithms that can handle $(P)$ were introduced. In particular, methods were designed to solve the more general class of mixed-integer nonlinear programs. These methods are branch-and-bound algorithms based on a convex relaxation of $(P)$. More precisely, in a first step a convex relaxation is designed and then a branch-and-bound is performed based on this relaxation. The most classical relaxation consists in the complete linearization of $(P)$, but quadratic convex relaxations can also be used. For instance, the well known $\alpha-$branch-and-bound [2] computes convex under-estimators of nonlinear functions by perturbing the diagonal of the Hessian matrix of the objective function. Several implementations of these algorithms are available, see for instance Baron [41], Antigone [36], SCIP [1] or Couenne [6].

In the case where the objective function is a polynomial, but the variables are continuous, Lasserre proposes in [29] an algorithm based on a hierarchy of semi-definite relaxations of $(P)$. The idea is, at each rank of the hierarchy, to successively tighten semi-definite relaxations of $(P)$ in order to reach its optimal solution value. It is also proven in [29] that this hierarchy converges in a finite number of iterations to the optimal solution of the considered problem. Further, this work has been extended to hierarchies of second order conic programs [3, 21, 28], and of sparse doubly non-negative relaxation [25]. Although these algorithms were not originally tailored for binary programming, they can handle $(P)$ by considering the quadratic constraint $x_{i}^{2}=x_{i}$. Methods devoted to the binary polynomial case were also proposed. In [14, 30], the authors use separable or convex under-estimators to approximate a given polynomial. Other methods based on linear reformulations can be found in [17, 18, 42], in which linear equivalent formulations to $(P)$ are proposed and then improved. In [15], the authors focus on a polyhedral description of the linearization of a binary polynomial program. Finally, the work in [4] considers quadratizations with a minimal number of additional variables.

In this paper, we focus on finding equivalent quadratic convex formulations of $(P)$. Quadratic convex reformulation methods [8, 9] were introduced for the specific case where $d=2$. The idea of these approaches is to build tight equivalent reformulations to $(P)$ that have a convex objective function. This equivalent problem can be built using the dual solution of a semi-definite relaxation of $(P)$, and further solved by a branch-and-bound algorithm based on quadratic convex relaxation. Here, we consider the more general case where $d\geq 3$, and we propose to compute an equivalent convex formulation to $(P)$. Hence, we present an exact solution method for problem $(P)$ that can be split in three phases. The first phase consists in building an equivalent formulation to $(P)$ where both objective function and constraints are at most quadratic. For this, we need to add some auxiliary variables. We then obtain problem $(QP)$ that has a quadratic objective function and linear inequalities.

Then in the second phase, we focus on the convexification of the obtained problem. As illustrated in the experiments of Section 4, the original QCR and MIQCR methods are not able to handle $(QP)$. Indeed, QCR leads to a reformulation with a weak bound, and in method MIQCR the semi-definite program that we need to solve is too large. This is why, in this paper we introduce a tailored convexification phase. The idea is to apply convex quadratic reformulation to any quadratization of $(P)$. For this, we need null quadratic functions on the domain of $(QP)$ so as to perturb the Hessian matrix of the new quadratic objective function. One of these null functions comes from the classic binary identity, $x_{i}^{2}=x_{i}$. One contribution of this paper is the introduction of new null quadratic functions on the domain of $(QP)$. This set of functions varies according to the quadratization used in phase 1. Adding these functions to the new objective function, we get a family of convex equivalent formulations to $(QP)$ that depend on some parameters. We then want to choose these parameters such that the continuous relaxation bound of the convexified problem is maximized. We show that they can be computed thanks to a semi-definite program. Finally, the last phase consists in solving the convexified problem using general-purpose optimization software.

Our experiments show that PQCR is able to solve to global optimality $6$ unsolved instances of the low auto-correlation binary sequence problem and improves lower and/or upper bounds of $27$ of the $45$ instances available at the minlplib website, in comparison to the other available solvers.

The outline of the paper is the following. In Section 2, we define and present our quadratizations of $(P)$. In Section 3, we introduce our family of convex reformulations and we prove how we compute the best parameters. Then, in Section 4, we present our computational results on polynomial instances of degree $4$ coming from the literature and we discuss different possible quadratizations of $(P)$. Section 5 draws a conclusion.

2 Phase 1: Quadratization of $(P)$

In this section, we present how we build equivalent quadratic formulations to $(P)$. The basic idea is to reduce the degree of $f$ to $2$. For this, in each monomial of degree 3 or greater, we simply recursively replace each product of two variables by an additional variable.

More formally, we define the set of indices of the additional variables $J=\{n+1,..,N\}$, where $N$ is the total number of initial and additional variables. We also define the subsets $\mathcal{E}_{i}$ for the initial or additional variable $i$ as follows:

Definition 1

For all $i\in I\cup J$, we define $\mathcal{E}_{i}$ as the set of indices of the variables whose product is equal to $x_{i}$:

• •

If $i\in I$, i.e. $x_{i}$ is an initial variable, then we set $\mathcal{E}_{i}=\{i\}$

• •

If $i\in J$, i.e. $x_{i}$ is an additional variable, then there exist $(i_{1},i_{2})\in(I\cup J)^{2}$ such that $x_{i}$ replaces $x_{i_{1}}x_{i_{2}}$ and we set $\mathcal{E}_{i}=\mathcal{E}_{i_{1}}\cup\mathcal{E}_{i_{2}}$

∎

Using these sets, we define a valid quadratization as a reformulation with $N$ variables where any monomial of degree at least 3 is replaced by the product of two variables.

Definition 2

The sets $J=\{n+1,..,N\}$ and $\{\mathcal{E}_{i},\leavevmode\nobreak\ i\in I\cup J\}$ define a valid quadratization with $N$ variables if, for any monomial $p$ of degree greater than or equal to $3$ (i.e. $|\mathcal{M}_{p}|\geq 3$), there exist $(j,k)\in(I\cup J)^{2}$ such that $\mathcal{M}_{p}=\mathcal{E}_{j}\cup\mathcal{E}_{k}$ and $\prod\limits_{\begin{subarray}{c}i\in\mathcal{M}_{p}\end{subarray}}x_{i}=x_{j}x_{k}$. Then the monomial $p$ is replaced by a quadratic term.

∎

With this definition of a quadratization, we reformulate $(P)$ as a non-convex quadratically constrained quadratic program $(QCQP)$ with $N$ variables.

[TABLE]

As the variables are binary, Constraints (1) are equivalent to the classical set of Fortet inequalities [18]:

[TABLE]

We now define set $\mathcal{F}_{\mathcal{E}}$:

[TABLE]

We denote by $M=4(N-n)$ the number of constraints of $\mathcal{F}_{\mathcal{E}}$. We thus obtain the following linearly constrained quadratic formulation that is equivalent to $(P)$ and has $N$ variables and $M$ constraints:

[TABLE]

where $Q\in S_{N}$ (the set of $N\times N$ real symmetric matrices), and $c\in\mathbb{R}^{N}$.

In the following, we will focus on the solution of problem $(QP)$ that is an equivalent formulation to $(P)$. Let us observe that $(QP)$, as well as $(QCQP)$, are parameterized by the quadratization defined by sets $\mathcal{E}$. Indeed, several valid quadratizations can be applied to $(P)$, each of them leading to different sets $\mathcal{E}_{i}$.

Different valid quadratizations were introduced and compared from the size point of view in [4]. In our case the comparison criterion is the continuous relaxation bound value from which we present our experimental comparison in Section 4.

**Example 1  ***[Different valid quadratizations]  

Let us consider the following problem:*

[TABLE]

For instance, we can build three different equivalent functions:

• •

$g_{1}(x)=2x_{1}+3x_{2}x_{3}-2\underbrace{x_{2}x_{4}}_{x_{5}}x_{3}-3\underbrace{x_{1}x_{4}}_{x_{6}}\underbrace{x_{2}x_{3}}_{x_{7}}$**

• •

$g_{2}(x)=2x_{1}+3x_{2}x_{3}-2\underbrace{x_{3}x_{4}}_{x_{5}}x_{2}-3\underbrace{x_{1}x_{2}}_{x_{6}}\underbrace{x_{3}x_{4}}_{x_{5}}$**

• •

$g_{3}(x)=2x_{1}+3x_{2}x_{3}-2\underbrace{x_{2}x_{4}}_{x_{5}}x_{3}-3\underbrace{x_{1}x_{2}}_{x_{6}}\underbrace{x_{3}x_{4}}_{x_{7}}$**

[TABLE]

Here we obtain 3 different quadratizations of $(Ex)$ with different sets $\mathcal{E}$. They have different sizes: $(QEx_{1})$ and $(QEx_{3})$ have 7 variables and 12 constraints, while $(QEx_{2})$ has 6 variables and 8 constraints.

$\Box$

We have reduced the degree of the polynomial program $(P)$ by building an equivalent quadratic program to $(P)$. However, the solution of $(QP)$ still has two difficulties, the non-convexity of the objective function $g$ and the integrality of the variables.

Some state-of-the-art solvers can solve $(QP)$ to global optimality (e.g. Cplex 12.7 [24]). Unfortunately, these solvers may not be enough efficient for solving dense instances of $(P)$. Here, we propose to compute an equivalent quadratic convex formulation to $(QP)$. There exist several convexification methods devoted to quadratic programming (see, for example [9, 11, 16, 22, 34]). These approaches can be directly applied to $(QP)$. For instance, one can use the QCR method, described in [11], that consists in computing an equivalent convex formulation to $(QP)$ using semi-definite programming. The convexification is obtained thanks to a non uniform perturbation of the diagonal of the Hessian matrix. The semi-definite relaxation used can be easily solved due to its reasonable size. However, the bound obtained by continuous relaxation of the reformulation is very weak. As a consequence, for the considered instances of Section 4, the branch-and-bound used to solve the reformulation failed as soon as $n\geq 20$. Another alternative is to apply the MIQCR method [9]. In this method, the perturbation is generalized to the whole Hessian matrix and hence is more refined than the previous one. This leads to a reformulation with a significantly sharper bound. Unfortunately, the semi-definite relaxation used in this approach is too large and its computation failed even with instances of $(P)$ containing only 10 variables. In the next section, we present a new convexification that leads to sharper bounds than QCR but with a better tractability than MIQCR.

3 Phase 2: A quadratic convex reformulation of $(QP)$

In this section, we consider the problem of reformulating $(QP)$ by an equivalent quadratic 0-1 program with a convex objective function. To do this, we define a new convex function which value is equal to the value of $g(x)$, but which Hessian matrix is positive semi-definite. More precisely, we first add to $g(x)$ a combination of four sets of functions that vanish on the feasible set $\mathcal{F}_{\mathcal{E}}$. For each function we introduce a scalar parameter. Then we focus on computing the best parameters that lead to a convex function and that maximize the optimal value of the continuous relaxation of the obtained problem.

3.1 Valid quadratic equalities for $(QP)$

For a quadratization characterized by $\mathcal{E}$, we introduce null quadratic functions over the set $\mathcal{F}_{\mathcal{E}}$.

Lemma 1

*The following quadratic equalities characterize null functions over $\mathcal{F}_{\mathcal{E}}$: *

[TABLE]

*Proof.  * Constraints (2) trivially hold since $x_{i}\in\{0,1\}$. Constraints (4) come from Definition 1. We then prove the validity of the Constraints (3) and (5).

• •

Constraints (3): we have $x_{i}=\displaystyle{\prod_{i^{\prime}\in\mathcal{E}_{i}}}x_{i^{\prime}}$ and $x_{j}=\displaystyle{\prod_{j^{\prime}\in\mathcal{E}_{j}}x_{j^{\prime}}}$, then:

[TABLE]

• •

Constraints (5): by definition we have:

[TABLE]

$\Box$

3.2 An equivalent quadratic convex reformulation to $(QP)$

We now compute a quadratic convex reformulation of $(QP)$ and thus of $(P)$. For this, we add to the objective function $g$ the null quadratic forms in (2)–(5). For each of them, we associate a real scalar parameter: $\alpha_{i}$ for Constraints (2), $\beta_{ij}$ for Constraints (3), $\delta_{ijk}$ for Constraints (4), and $\lambda_{ijkl}$ for Constraints (5). We get the following parameterized function:

[TABLE]

Obviously $g_{\alpha,\beta,\delta,\lambda}(x)$ has the same value as $g(x)$ for any $x\in\mathcal{F_{E}}$. Moreover, there exist vector parameters $\alpha$, $\beta$, $\delta$ and $\lambda$ such that $g_{\alpha,\beta,\delta,\lambda}$ is a convex function. Take for instance, $\alpha$ equals to the opposite of the smallest eigenvalue of $Q$, and $\beta=\delta=\lambda=0$.

By replacing $g$ by the new function, we obtain the following family of quadratic convex equivalent formulation to $(QP)$:

[TABLE]

where $Q_{\alpha,\beta,\delta,\lambda}\in S_{N}$ is the Hessian matrix of $g_{\alpha,\beta,\delta,\lambda}(x)$, and $c_{\alpha,\beta,\delta,\lambda}\in\mathbb{R}^{N}$ is the vector of linear coefficients of $g_{\alpha,\beta,\delta,\lambda}(x)$.

In order to use $(QP_{\alpha,\beta,\delta,\lambda})$ within a branch-and-bound procedure, we are interested by parameters ($\alpha$, $\beta$, $\delta$, $\lambda$) such that $g_{\alpha,\beta,\delta,\lambda}$ is a convex function. Moreover, in order to have a good behavior of the branch-and-bound algorithm, we want to find parameters that give the tightest continuous relaxation bound. More formally, we want to solve the following optimization problem:

[TABLE]

where $T_{1}$, $T_{2}$ and $T_{3}$ are the number of Constraints (3), (4), and (5), respectively, and $\overline{\mathcal{F}}_{\mathcal{E}}$ is the set $\mathcal{F}_{\mathcal{E}}$ where the integrality constraints are relaxed, i.e. $x\in[0,1]^{N}$.

In the rest of the paper we will focus on solving $(CP)$. For this, we build a compact semi-definite relaxation that uses our new valid equalities and prove that its optimal dual variables provide an optimal solution to $(CP)$.

3.3 Computing an optimal solution to $(CP)$

The following theorem shows that problem $(CP)$ is equivalent to the dual of a semi-definite relaxation of $(QP)$.

Theorem 3.1

The optimal value of $(CP)$ is equal to the optimal value of the following semi-definite program $(SDP)$:

[TABLE]

The optimal values $(\alpha^{*},\beta^{*},\delta^{*},\lambda^{*})$ of problem $(CP)$ are given by the optimal values of the dual variables associated with constraints (6)–(9) respectively.

*Proof.  * For simplicity, we rewrite $\mathcal{F}_{\mathcal{E}}$ as follows: $\mathcal{F}_{\mathcal{E}}=\{x\in\{0,1\}^{N}:Ax\leq b\}$ where $A$ is a $M\times N$-matrix, $b\in\mathbb{R}^{M}$, and we introduce $T=N+T_{1}+T_{2}+T_{3}$ the number of Constraints (2)–(5) respectively.

We start by observing that $x\in[0,1]^{N}$ is equivalent to $x^{2}\leq x$, thus, $(CP)$ is equivalent to $(Q1)$:

[TABLE]

$(Q1)$ is a convex optimization problem over a convex set. If we consider the solution $\tilde{x}_{i}=0.5$ $\forall i\in I$ and $\tilde{x}_{i}=\tilde{x}_{j}\tilde{x}_{k}$ $\forall(i,j,k)\in J\times(I\cup J)^{2},\leavevmode\nobreak\ \mathcal{E}_{i}=\mathcal{E}_{j}\cup\mathcal{E}_{k}$, the obtained $\tilde{x}$ is an interior point and the Slater’s conditions are satisfied for the minimization sub-problem. Then, by Lagrangian duality, we have $(Q1)$ equivalent to $(Q2)$:

[TABLE]

Due to Constraints (2), it holds that $(Q2)$ is equivalent to $(Q3)$:

[TABLE]

It is well known that a necessary condition for the quadratic function $g_{\alpha,\beta,\delta,\lambda,\gamma}(x)+\gamma^{T}(Ax-b)$ to have a minimum not equal to $-\infty$ is that matrix $Q_{\alpha,\beta,\delta,\lambda}$ is positive semi-definite. Therefore $(Q3)$ is equivalent to $(Q4)$:

[TABLE]

We know from [32] that $(Q4)$ is equivalent to problem $(D)$:

[TABLE]

By semi-definite duality of program $(D)$, and with $\alpha,\beta,\delta,\lambda$ the dual variables associated with Constraints (6)–(9) respectively, we get $(SDP^{\prime})$:

[TABLE]

We now prove that there is no duality gap between $(D)$ and $(SDP^{\prime})$, which holds since:

1. (i)

The feasible domain of $(SDP^{\prime})$ is nonempty, as $(QP_{\alpha,\beta,\delta,\lambda})$ contains 0 as a feasible solution and $(D)$ is bounded 2. (ii)

$(D)$ satisfies Slater’s condition. It is sufficient to take $\beta$, $\delta$ and $\lambda$ equal to [math], $\alpha$ large enough so that $Q_{\alpha,\beta,\delta,\lambda}\succeq 0$ holds, and $t$ a large negative number that ensures the diagonal dominance of the first row and the first column of matrix \left(\begin{array}[]{c c}-\gamma^{T}b-t&\frac{1}{2}(c^{T}_{\alpha,\beta,\delta,\lambda}+\gamma^{T}A)\\ \frac{1}{2}(c_{\alpha,\beta,\delta,\lambda}+A^{T}\gamma)&Q_{\alpha,\beta,\delta,\lambda}\end{array}\right).

From these equivalences, we know that we can build an optimal solution of $(CP)$ from the optimal dual variables of $(SDP^{\prime})$. However, constraints $Ax\leq b$ are redundant in $(SDP^{\prime})$ and we thus prove in Lemma 2 that $(SDP^{\prime})$ and $(SDP)$ are equivalent. As a consequence, an optimal solution to $(CP)$ can be deduced from the optimal dual variables of $(SDP)$.

Lemma 2

Due to Constraints (6)–(8) and (12), inequalities $Ax\leq b$ are redundant in $(SDP^{\prime})$.

*Proof.  * Recall that $Ax\leq b$ are the inequalities of $(C_{j,k}^{i}),\,\forall(i,j,k)\in J\times(I\cup J)^{2}:\leavevmode\nobreak\ \mathcal{E}_{i}=\mathcal{E}_{j}\cup\mathcal{E}_{k}$, i.e $x_{i}\geq 0$ $(a)$, $x_{i}\leq x_{j}$ $(b)$, $x_{i}\leq x_{k}$ $(c)$, and $x_{i}\geq x_{j}+x_{k}-1$ $(d)$.

The basic idea used here is that, since matrix \left(\begin{array}[]{c c}1&x^{T}\\ x&X\end{array}\right) is positive semi-definite, all its symmetric minors are non-negative.

• •

Constraint $(a)$: $x_{i}\geq 0$. We consider the determinant $\begin{vmatrix}1&x_{i}\\ x_{i}&X_{ii}\end{vmatrix}$, which implies $X_{ii}-x_{i}^{2}\geq 0$. By (6) we obtain $x_{i}-x_{i}^{2}\geq 0$ and thus $x_{i}\geq 0$.

• •

Constraint $(b)$: $x_{i}\leq x_{j}$. Considering the determinant of the symmetric minor $\begin{vmatrix}X_{jj}&X_{ji}\\ X_{ij}&X_{ii}\end{vmatrix}$ implies $X_{ii}X_{jj}-X_{ij}^{2}\geq 0$. By (6) we have $x_{j}x_{i}-X_{ij}^{2}\geq 0$ and by (7) we obtain $x_{i}x_{j}-x_{i}^{2}\geq 0$. Either $x_{i}>0$ and then we have $x_{j}-x_{i}\geq 0$, or $x_{i}=0$ and the inequality comes from $x_{j}\geq 0$.

• •

Constraint $(c)$: $x_{i}\leq x_{k}$. By symmetry, i.e. considering the determinant $\begin{vmatrix}X_{kk}&X_{ki}\\ X_{ik}&X_{ii}\end{vmatrix}$ , the inequality holds.

• •

Constraint $(d)$: $x_{i}\geq x_{j}+x_{k}-1$. By definition (12) implies z^{T}\left(\begin{array}[]{c c}1&x^{T}\\ x&X\end{array}\right)z\geq\leavevmode\nobreak\ 0,\\ \forall z\in\mathbb{R}^{N+1}. By taking $\bar{z}=(1,0,..,0,\underbrace{-1}_{j},0,..,0,\underbrace{-1}_{k},0,..,0,\underbrace{1}_{i},0,..,0)$, we have:

[TABLE]

$\Box$

Let us state Corollary 1 that shows that from an optimal dual solution to $(SDP)$ we can build an optimal solution to $(CP)$.

Corollary 1

We have $v(CP)=v(SDP)$ where $v(.)$ is the optimal value of problem $(.)$. Consequently, an optimal solution $(\alpha^{*},\beta^{*},\delta^{*},\lambda^{*})$ of $(CP)$ corresponds to the optimal values of the dual variables associated with constraints (6)–(9) of $(SDP)$ respectively.

*Proof.  * We have:

1. (i)

$v(CP)=v(D)$ 2. (ii)

since there is no duality gap between $(D)$ and $(SDP^{\prime})$, we have $v(D)=v(SDP^{\prime})$ 3. (iii)

by Lemma 2, we get $v(CP)=v(D)=v(SDP^{\prime})=v(SDP)$

$\Box$

$\Box$

To sum up, we obtain $(QP^{*})$, the best equivalent convex formulation to $(QP)$:

[TABLE]

From Theorem 1, we deduce the Algorithm 1 to solve $(P)$.

4 Numerical results

In this section, we evaluate PQCR on two applications. The first one is the image restoration (vision) problem [17], which results are presented in Section 4.1. The instances of this application are quite sparse with an average ratio $\frac{m}{n}$ of about $7$. We choose to use these instances in order to compare PQCR  with existing convexifications and in particular with methods QCR and MIQCR that are not able to handle larger and/or denser instances. Then, in Section 4.2, we present the results of the second application, the low auto-correlation binary sequence (LABS) problem [7] which instances are much denser (average ratio $\frac{m}{n}$ of about $212$). These instances are available on the minlplib website [35], and are very hard to solve. For most of them, the optimal solution value is not known. For these experiments, we have chosen the quadratization described in Algorithm 2 for Step 1 of PQCR. This choice impacts the number of constraints within the sets $\mathcal{F}_{\mathcal{E}}$ and $\mathcal{S}_{\mathcal{E}}$, and the associated continuous relaxation bound value can thus vary. We further illustrate this variation on toy instances in Section 4.3.

The quadratization used in our experiments is presented in Algorithm 2.

Example 1

Applying the quadratization of Algorithm 2 to the monomial $x_{1}x_{2}x_{3}x_{4}x_{5}$ we obtain the following monomial of degree $3$ at the first iteration:

[TABLE]

we then obtain a quadratic reformulation of the monomial at the second iteration using $3$ additional variables:

[TABLE]

∎

Our experiments were carried out on a server with $2$ CPU Intel Xeon each of them having $12$ cores and $2$ threads of $2.5$ GHz and $4*16$ GB of RAM using a Linux operating system. For all algorithms, we used the multi-threading version of Cplex 12.7 with up to 48 threads.

In our experiments, we use three classes of solution algorithms:

1. i)

The first class includes 3-phase algorithms that consist in a quadratization and a convexification followed by the solution of the equivalent convex problem with the solver Cplex 12.7: these methods are PQCR and Q+QCR. For both methods, the quadratization is implemented in C, and we used the solver csdp to solve the semi-definite programs. For denser instances (Section 4.2), we used the solver csdp [12] together with the Conic Bundle algorithm [23] to solve the semi-definite program of PQCR, as described in [10]. Then, we used the ampl [19] interface of the solver Cplex 12.7 [24] to solve the obtained quadratic convex problem with binary variables. 2. ii)

The second class includes a 2-phase algorithm, called Q+Cplex, that consists in a quadratization followed by the direct submission to Cplex 12.7. Here again, the quadratization is implemented in C, and we used the ampl interface of the solver Cplex 12.7. 3. iii)

The third class includes the direct submission to the general mixed-integer non-linear solver Baron 17.4.1 [41]. Here, we used the gams interface of the solver Baron 17.4.1.

Parameters of the solvers

• •

Cplex : we let the default parameters, except the parameter qptolin that is set to 0 for methods PQCR  and Q+QCR.

• •

csdp : Parameters axtol, aytol of Csdp are set to $10^{-3}$.

• •

Conic Bundle : the precision is set to $10^{-3}$. Parameter $p$ (see [10]) is set to $0.2*|\mathcal{F_{E}}|$.

Legends of Tables 1-3

• •

Name: Name of the considered instance.

• •

n: number of variables in the polynomial formulation.

• •

m: number of monomials.

• •

BKN: is the optimal solution value or the best known solution value of the instance.

• •

N: number of variables after quadratization.

• •

gap: is the initial gap, i.e. the gap at the root node of the branch-and-bound, $gap=\displaystyle{\left|\frac{BKN-LB_{i}}{BKN}\right|*100}$, where $LB_{i}$ is the initial lower bound.

• •

Solution: best solution value found within the time limit.

• •

tSdp: CPU time in seconds for solving semi-definite programs in PQCR and Q+QCR. The time limit is set to 2400 seconds for the vision problem and 3 hours for the LABS problem. If the solver reaches the time limit, tSdp is labeled as "-".

• •

tTotal: total CPU time in seconds of the associated method. The time limit is set to 1 hour for the vision problem and 5 hours for the LABS problem. If an instance remains unsolved within the time limit, we put the final gap$=\displaystyle\left|\frac{BKN-LB_{f}}{BKN}\right|*100$, where $LB_{f}$ is the final lower bound.

• •

Nodes: number of nodes visited by the branch-and-bound algorithm.

4.1 The image restoration problem

The vision instances are inspired from the image restoration problem, which arises in computer vision. The goal is to reconstruct an original sharp base image from a blurred image. An image is a rectangle containing $n=l\times h$ pixels. This rectangle is modeled as a binary matrix of the same dimension. A complete description of these instances can be found in [17]. The problem is modeled by the minimization of a degree $4$ polynomial of binary variables where each variable represents a pixel. The coefficients of the monomials are indicative of how likely a configuration is to appear on the sharp base image. The size of the considered instances are $l\times h=10\times 10$, $10\times 15$, and $15\times 15$, or in the polynomial formulation $n=100,150$ and $225$, with a number of monomials of $m=668$, $1033$, and $1598$ respectively. In our experiments, 15 instances of each size are considered obtaining a total of $45$ instances. Observe that the $15$ instances of the same size have identical monomials with different coefficients, because they represent different images with the same number of pixels.

We now focus on the comparison of several convexification methods after quadratization. Indeed, several ways are possible to solve the quadratic non-convex program $(QP)$. For instance, the standard solver Cplex can directly handle it, or one can apply the QCR [11] or MIQCR [9] methods. We compare PQCR with these three approaches. We do not report the results for method Q+MIQCR since it was not able to start the computation due to the size of the considered instances. We also give the computational results coming from the direct submission of $(P)$ to the solver Baron 17.4.1. Our observations for these instances are summed up in Table 1, where each line corresponds to one instance, where the $i^{th}$ instance of $l\times h$ pixels is labeled v.l.h i.

We start by comparing the convexification phase of our new algorithm with the original QCR and MIQCR methods. We observe that none of these convexifications are able to handle any considered instances: QCR because of the weakness of its initial gap, and MIQCR because of the size of the semidefinite problem considered for computing the best reformulation. These experiments confirm the interest of designing PQCR, an algorithm devoted to polynomial optimisation. Then, we can see that Q+Cplex dominates PQCR. However, one have to note that these instances are very sparse (average ratio $\frac{m}{n}$ of about $7$). It is well known that the standard linearization performs very well on sparse instances. Clearly, for these instances, the time spent on solving a large semidefinite program, even once, is not profitable in comparison to the efficiency of LP heuristic or cut methods implemented in cplex 12.7. Indeed, Q+Cplex solves all the considered instances at the root node of its branch-and-bound. Moreover, it is interesting to remark that $99\%$ of the CPU time of PQCR is spent for solving $(SDP)$, while the CPU time for solving $(QP^{*})$ is always smaller than $14$ seconds. Finally, we compare PQCR with the direct submission to the solver Baron. We observe that Baron is faster than PQCR on the medium size instances ($n=100$ or $150$), but is not able to solve all the larger instances within the time limit. Indeed, for $n=225$, it solves only $3$ instances out of $15$. On the contrary, PQCR seems quite stable to the increase of the size of the instances. Indeed, the initial gap remains stable ($0.42\%$ on average) while the total CPU time increases reasonably.

4.2 The Low Auto-correlation Binary Sequence problem

We consider the problem of binary sequences with low off-peak auto-correlations. More formally, let $S$ be a sequence $S=(s_{1},\ldots,s_{n})$ with $s\in\{-1,1\}^{n}$, and for a given $k=0,\ldots,n-1$, we define the auto-correlations $C_{k}(S)$ of $S$:

[TABLE]

The problem is to find a sequence $S$ of length $n$ that minimizes $E(S)$, a degree $4$ polynomial:

[TABLE]

This problem has numerous practical applications in communication engineering, or theoretical physics [7]. For our experiments, we consider truncated instances, i.e. sequences of length $n$ where we compute low off-peak auto-correlation up to a certain distance $n_{0}\leq n$, i.e. we consider the following function to minimize:

[TABLE]

In order to apply PQCR, which is initially tailored for $\{0,1\}$ polynomial programs, we convert the variables from $\{-1,1\}$ to $\{0,1\}$ using the standard transformation $x=\frac{s+1}{2}$.

This problem admits a lot of symmetries. In particular the correlations $C_{k}$ are identical for a sequence $S$ and its complement. We exploited this symmetry by fixing to 0 the variable that appears the most. Each instance is labeled b.n.n${}_{\texttt{0}}$. These instances were introduced by [33] and can be found on the minlplib [35] or the polip [38] websites. We do not report the results for methods Q+QCR and Q+MIQCR since they have failed to solve all the considered instances. Two instances that are already quadratic (b.20.03 and b.25.03) are solved by the method Q+Cplex in $7$ and $75$ seconds respectively. However, this method was not able to solve the other instances within the time limit.

We present in Table 2 a detailed comparison of PQCR with the direct submission to baron 17.1.4. For these experiences the total time limit was set to 5 hours, and we limit the CPU time for solving $(SDP)$ to 3 hours. Indeed, any feasible solution to the dual of $(SDP)$ can be used to get a convex objective function in the equivalent formulation. Thus, if the CPU time in column $tSdp$ is smaller than three hours it means that $(SDP)$ was solved to optimality. In the other case, we get a feasible dual solution and we can suppose that the initial gap of PQCR could be improved. For these instances PQCR is faster than baron since it solves $17$ instances out of $45$ within the time limit while baron solves only $13$ instances. Here, baron 17.1.4 solves $2$ instances that were stated as unsolved on minlplib. As expected, PQCR has an initial gap much smaller that baron (reduced by a factor $22$ on average). We also observe that the number of nodes visited by PQCR during the branch-and-bound is significantly larger than the the number of nodes of baron (increased by a factor of about $40000$ on average).

We present in Table 3 the values of the best solutions and of the final lower bounds obtained by PQCR within 5 hours of CPU time, and those available on the minlplib website. More precisely, we report in the column minlplib the best solution/final lower bound value obtained among the results of the solvers Antigone, Baron, Couenne, Lindo, and Scip. PQCR solves to optimality 6 unsolved instances (labeled as ∗∗). It also improves the best known solution values of 9 instances (labeled as #), and improves the final lower bound of all the unsolved instances (labeled as ∗). In this table, each line corresponds to one instance, and we only present results for instances that were stated as unsolved on minlplib.

To illustrate these results, we plot in Figure 1, for each instance reported in Table 3, the final gap of PQCR and minlplib. Clearly, the final gap of PQCR is much smaller than the final gap of minlplib (reduced by a factor $3$ on average).

A last remark concerns the CPU time necessary to solve $(SDP)$. Indeed, this time represents on average 75% of the total CPU time. A natural improvement is to identify the set of "important equalities" in a preprocessing step in order to improve the behavior of the solution of $(SDP)$. Obviously, this step should be dependent on the quadratization.

4.3 A short discussion on the impact of the chosen quadratization

In this section, we shortly explore the impact of the chosen quadratization on the tightness of the associated continuous relaxation bound. In Table 4, we report the continuous relaxation bound values obtained by convexification after applying the quadratization of Algorithm 2, and three quadratizations from [17], namely Pairwise Cover 1, 2 and 3 (PC1, PC2 and PC3). In Pairwise Cover 1, for each monomial of degree $d\geq 3$, the first two variables are linearized to obtain a monomial of degree $d-1$. The process is recursively reproduced until $d=2$. Pairwise Covers 2 and 3 try to minimize the number of additional variables. In PC2, the authors compute the sub-monomials of any degree that appear the most among all the intersection of pairs of monomials. Then they linearize these sub-monomials using the set $\mathcal{F_{E}}$ and they repeat the process until the objective function is quadratic. PC3 linearizes in priority the pair of variables that occurs the most frequently in all the monomials. For instance, if we consider the quadratization of the following monomial of degree $4$, $x_{1}x_{2}x_{3}x_{4}$, we will compute the most frequent pair of variables among the six possible products. If $x_{1}x_{2}$ is the most frequent, then the monomial will be quadratized using two variables, one for the reformulation of $x_{1}x_{2}$ and the other for $x_{3}x_{4}$.

We observe that the chosen quadratization impacts $N$, the number of variables of $(QP)$. It also impacts the quality of the associated semidefinite bound, $LB_{i}$. Indeed, the more variables are added, the more the size of sets $\mathcal{F_{E}}$ and $\mathcal{S_{E}}$ increases. Clearly, some equalities of $\mathcal{S_{E}}$ may be stronger than others. Interesting future research directions would be to identify, for a given quadratization, a set of "important" equalities in $\mathcal{S_{E}}$, and to determine which quadratization used in PQCR leads to faster solution time and/or sharper initial lower bound.

5 Conclusion

We consider the general problem $(P)$ of minimizing a polynomial function where the variables are binary. In this paper, we present PQCR a solution approach for $(P)$. PQCR can be split in 3 phases. We called the first phase quadratization, where we rewrite $(P)$ as an equivalent quadratic program $(QP)$. For this we have to add new variables and linear constraints. We get a linearly constrained quadratic program that still has a non-convex objective function and binary variables. Moreover, even for small instances of $(P)$, the existing convexification methods failed to solve the associate $(QP)$. This is why, we present a family of tailored quadratic convex reformulations of $(QP)$ that exploits its specific structure. For this, we introduce new valid quadratic equalities that vanish on the feasible domain of $(QP)$. We use these equalities to build a family of equivalent quadratic convex formulations to $(QP)$. Then, we focus on finding, within this family, the equivalent convex formulation that maximizes the continuous relaxation bound value. We show that we can compute this "best" convex reformulation using a new semidefinite relaxation of $(QP)$. Finally, we solve our optimal reformulation with a standard solver.

We present computational results on two applications and compare our algorithm with other convexification methods and the general solver Baron. In particular, we show that for the low auto-correlation binary sequence problem, PQCR is able to improve the best known solution of 10 instances out of 45. A future research direction would be to characterize which quadratization best fit with our convexification phase from the continuous relaxation value point of view.

Acknowledgment

The authors are thankful to Elisabeth Rodriguez-Heck and Yves Crama for sharing the executable of their quadratization code in order to compute the Pairwise Covers 1, 2 and 3 on LABS instances.