A unified approach to mixed-integer optimization problems with logical constraints

TL;DR

This paper introduces a unified, non-linear formulation for mixed-integer optimization with logical constraints, enabling faster and larger-scale solutions than traditional linear methods across various applications.

Contribution

It presents a novel non-linear reformulation and regularization approach that transforms complex logical mixed-integer problems into convex binary problems, improving scalability and solution quality.

Findings

Solves network design problems with hundreds of nodes faster and better.

Handles sparse portfolio selection with up to 3,200 securities.

Addresses sparse regression with up to 100,000 covariates.

Abstract

We propose a unified framework to address a family of classical mixed-integer optimization problems with logically constrained decision variables, including network design, facility location, unit commitment, sparse portfolio selection, binary quadratic optimization, sparse principal analysis and sparse learning problems. These problems exhibit logical relationships between continuous and discrete variables, which are usually reformulated linearly using a big-M formulation. In this work, we challenge this longstanding modeling practice and express the logical constraints in a non-linear way. By imposing a regularization condition, we reformulate these problems as convex binary optimization problems, which are solvable using an outer-approximation procedure. In numerical experiments, we establish that a general-purpose numerical strategy, which combines cutting-plane, first-order and local search methods, solves these problems faster and at a larger scale than state-of-the-art mixed-integer linear or second-order cone methods. Our approach successfully solves network design problems with 100s of nodes and provides solutions up to 40\% better than the state-of-the-art; sparse portfolio selection problems with up to 3,200 securities compared with 400 securities for previous attempts; and sparse regression problems with up to 100,000 covariates.

Tables

Table 1. Table 1 : Loss functions and Fenchel conjugates for ERM problems of interest.

Method	Loss function	Domain	Fenchel conjugate
OLS	$\frac{1}{2}{\left(y-u\right)}^2$	$y\in ℝ$	${\ell }^{\star }(y,\alpha )=\frac{1}{2}{\alpha }^2+\alpha y$
SVM	$\max (1-yu,0)$	$y\in \left\{\pm 1\right\}$	${\ell }^{\star }(y,\alpha )=\{\begin{array}{cc}\alpha y,\text{if}\alpha y\in [-1,0],& \\ \infty ,\text{otherwise}.& \end{array}$

Table 2. Table 2 : Summary of the advantages ( + + ) /disadvantages ( − - ) of both techniques.

Regularization	Characteristics
Big-$M$	($+$) Linear constraints
	($+$) Supplies the same objective if $M>M_1$, for some
	($-$) Leads to infeasible problem if , for some
Ridge	($+$) Strongly convex objective
	($-$) Systematically leads to a different objective for any $\gamma >0$
	($+$) Preserves the feasible set

Table 3. Table 3 : Best solution found after one hour on network design instances with m 𝑚 m nodes and ( 1 + p ) ​ m 1 𝑝 𝑚 (1+p)m initial edges. We report improvement, i.e., the relative difference between the solutions returned by CPLEX and the cutting-plane. Values are averaged over five randomly generated instances. For ridge regularization, we report the “unregularized” objective value, that is we fix 𝐳 𝐳 \bm{z} to the best solution found and resolve the corresponding sub-problem with big- M 𝑀 M regularization. A “ − - ” indicates that the solver could not finish the root node inspection within the time limit (one hour), and “Imp.” is an abbreviation of improvement.

			Big-$M$			Ridge			Overall
$m$	$p$	unit	CPLEX	Cuts	Imp.	CPLEX	Cuts	Imp.	Imp.
40	0	$\times {10}^9$	1.17	1.16	$0.86\%$	1.55	1.16	$24.38\%$	$1.74\%$
80	0	$\times {10}^9$	8.13	7.52	$6.99\%$	9.95	7.19	$26.74\%$	$10.85\%$
120	0	$\times {10}^{10}$	3.03	2.10	$29.94\%$	$-$	1.94	$-\%$	$35.30\%$
160	0	$\times {10}^{10}$	5.90	4.32	$26.69\%$	$-$	4.07	$-\%$	$30.91\%$
200	0	$\times {10}^{10}$	11.45	7.78	$31.45\%$	$-$	7.50	$-\%$	$32.32\%$
40	1	$\times {10}^8$	5.53	5.47	$1.07\%$	5.97	5.45	$8.74\%$	$1.41\%$
80	1	$\times {10}^9$	2.99	2.94	$1.81\%$	3.16	2.95	$6.78\%$	$1.89\%$
120	1	$\times {10}^9$	8.38	7.82	$6.69\%$	$-$	7.82	$-\%$	$6.86\%$
160	1	$\times {10}^{10}$	1.64	1.54	$5.98\%$	$-$	1.54	$-\%$	$6.03\%$
200	1	$\times {10}^{10}$	2.60	2.54	$2.33\%$	$-$	2.26	$-\%$	$12.98\%$
40	2	$\times {10}^8$	4.45	4.38	$1.62\%$	4.76	4.36	$8.27\%$	$2.06\%$
80	2	$\times {10}^9$	2.44	2.31	$5.39\%$	2.46	2.31	$5.97\%$	$5.40\%$
120	2	$\times {10}^9$	6.23	5.89	$5.55\%$	$-$	5.89	$-\%$	$5.75\%$
160	2	$\times {10}^{11}$	1.22	1.16	$4.74\%$	$-$	0.71	$-\%$	$19.33\%$
200	2	$\times {10}^{10}$	2.06	1.43	$30.46\%$	$-$	1.01	$-\%$	$73.43\%$
40	3	$\times {10}^8$	3.91	3.85	$1.58\%$	4.13	3.85	$6.73\%$	$1.78\%$
80	3	$\times {10}^9$	2.06	1.94	$5.76\%$	2.04	1.94	$5.44\%$	$5.85\%$
120	3	$\times {10}^9$	5.43	5.15	$5.31\%$	$-$	4.2	$-\%$	$12.35\%$
40	4	$\times {10}^8$	3.32	3.28	$1.35\%$	3.53	3.26	$7.71\%$	$1.85\%$
80	4	$\times {10}^9$	1.88	1.77	$5.59\%$	$-$	1.77	$-\%$	$5.64\%$

Table 4. Table 4 : Average runtime in seconds on binary quadratic optimization problems from the Biq-Mac library [ 53 , 15 ] . Values are averaged over 10 10 10 instances. A “ − - ” denotes an instance which was not solved because the approach did not respect the 32 32 32 GB peak memory budget.

Instance	$n$	Average runtime (s)/Average optimality gap $(\%)$
		CPLEX-M	CPLEX-M-Triangle	Cuts-M	Cuts-M-Triangle
bqp-$50$	$50$	$29.4$	$0.6$	$30.6$	0.4
bqp-$100$	$100$	$122.3$	$51.7$	$25.3\%$	38.6
bqp-$250$	$250$	$1108.1\%$	$83.5\%$	$87.0\%$	46.1%
bqp-$500$	$500$	$2055.8\%$	$1783.3\%$	157.3%	$410.7\%$
bqp-$1000$	$1000$	$-$	$-$	$\text{260.9}\%$	$-$
be$100$	$100$	79.7%	$208.0\%$	$249.4\%$	$201.2\%$
be$120.8$	$120$	146.4%	$225.8\%$	$264.1\%$	$220.3\%$

Table 5. Table 5 : Average incumbent objective value (higher is better) after 1 1 1 hour for medium-scale binary quadratic optimization problems from the Biq-Mac library [ 53 , 15 ] . “ − - ” denotes an instance which was not solved because the approach did not respect the 32 32 32 GB peak memory budget. Values are averaged over 10 10 10 instances. Cuts-Triangle includes an extended formulation in the master problem.

Instance	$n$	Average objective value
		CPLEX-M	CPLEX-M-Triangle	Cuts-M	Cuts-M-Triangle
bqp-$250$	$250$	$9920.8$	$41843.4$	43774.9	$43701.5$
bqp-$500$	$500$	$19417.1$	$19659.0$	122879.3	$122642.4$
bqp-$1000$	$1000$	$-$	$-$	351450.7	$-$
be$100$	$100$	$16403.0$	$16985.0$	$17152.1$	17178.5
be$120.8$	$120$	$17943.2$	$19270.3$	$19307.7$	19371.2

Table 6. Table 6 : Runtime in seconds per approach. We run all approaches on one thread, and impose a time limit of 600 600 600 s. If a solver fails to converge, we report the relative bound gap at termination in brackets, and the no. explored nodes and cuts at the time limit. For ridge regularization, we set γ = 100 / k 𝛾 100 𝑘 \gamma={100}/{k} .

Dataset	$p$	$k$	Big-$M$ regularization			Ridge regularization
			Time(s)	Nodes	Cuts	Time(s)	Nodes	Cuts
Pitprops	$13$	$5$	0.09	$45$	$22$	$0.42$	$42$	$16$
		$10$	0.08	$223$	$223$	$0.68$	$615$	$244$
Wine	$13$	$5$	0.04	$143$	$69$	$0.10$	$73$	$36$
		$10$	0.09	$364$	$232$	$0.61$	$394$	$230$
Miniboone	$50$	$5$	$0.03$	$3$	$6$	0.01	$0$	$2$
		$10$	0.04	$4$	$6$	$0.07$	$10$	$13$
Communities	$101$	$5$	0.15	$109$	$2$	$0.54$	$272$	$55$
		$10$	0.44	$373$	$76$	$2.20$	$1,800$	$328$
Arrhythmia	$274$	$5$	5.27	$1,080$	$192$	$6.75$	$1,242$	$282$
		$10$	(4.21%)	$61,000$	$11,600$	($4.63\%)$	$77,200$	$11,360$
Micromass	$1300$	$5$	131.3	$4,580$	$4$	$163.2$	$4$	$3,809$
		$10$	378.6	$321$	$16,090$	$510.3$	$21,700$	$566$

Table 7. Table 7 : Largest sparse portfolio instances reliably solved by each approach

Reference	Solution method	Largest instance size solved
		(no. securities)
[31]	Perspective cut+SDP	$400$
[16]	Nonlinear BB	$200$
[32]	Lagrangian relaxation BB	$300$
[18]	Lagrangian relaxation BB	$300$
[55]	SDP BB	$400$
[27]	Approx. Proj. Perspective Cut	$400$
[8]	Algorithm 1 with ridge regularization	$3,200$

Table 8. Table 8 : Proportion of wins and relative improvement over CPLEX in terms of computational time on the 112 instances from the OR-library [ 3 , 39 ] for different implementations of our method: an outer-approximation (OA) scheme with cuts generated at the root node using Kelley’s method (OA + Kelley), OA with the local search procedure (OA + Local search) and OA with a strategy for both the lower and upper bound (OA + Both). Relative improvement is averaged over all “win” instances.

	Big-$M$		Ridge
Algorithm	$\%$ wins	Relative improvement	$\%$ wins	Relative improvement
OA + Kelley	$1.8\%$	$36.6\%$	$30.1\%$	$91.6\%$
OA + Local search	$1.9\%$	$49.5\%$	$19.4\%$	$73.8\%$
OA + Both	$12.7\%$	$53.6\%$	$92.5\%$	$91.7\%$

Table 9. Table 9 : Average runtime in seconds per approach, on data from [ 29 ] where the quadratic cost are multiplied by a factor of α 𝛼 \alpha . If the method did not terminate in one hour, we report the bound gap. n 𝑛 n denotes the number of generators, each instances has 24 24 24 trade periods.

$\alpha$	$0.1$		$1$		$2$		$5$		$10$
$n$	Big-$M$	Ridge	Big-$M$	Ridge	Big-$M$	Ridge	Big-$M$	Ridge	Big-$M$	Ridge
$100$	93.6	$299.0$	16.2	$229.4$	$0.32\%$	47.9	$1.68\%$	4.6	$2.76\%$	6.0
$150$	35.6	$352.1$	6.2	$28.3$	$0.25\%$	33.4	$1.69\%$	6.4	$2.82\%$	8.0
$200$	56.3	$138.1$	3.3	$239.7$	$0.24\%$	112.9	$1.62\%$	16.7	$2.81\%$	21.2

Table 10. Table 10 : Optimality gap after one hour on network design instances with m 𝑚 m nodes and ( 1 + p ) ​ m 1 𝑝 𝑚 (1+p)m initial edges. We only report results for instances where the resulting gap was less than 100 % percent 100 100\% for at least one of the four approaches. A “ − - ” indicates that the solver could not finish the root node inspection within the time limit (one hour).

		Big-$M$		Ridge
$m$	$p$	CPLEX	Cuts	CPLEX	Cuts
40	0	$69.8\%$	$100\%$	$98.9\%$	$96.7\%$
80	0	$100\%$	$100\%$	$100\%$	$100\%$
40	1	$38.6\%$	$100\%$	$99.8\%$	$97.1\%$
80	1	$100\%$	$100\%$	$100\%$	$95.6\%$
120	1	$100\%$	$100\%$	$-$	$96.6\%$
40	2	$23.3\%$	$100\%$	$>100\%$	$97.7\%$
80	2	$100\%$	$100\%$	$100\%$	$96.3\%$
40	3	$74.6\%$	$>100\%$	$97.5\%$	$98.1\%$
40	4	$100\%$	$100\%$	$99.2\%$	$98.2\%$
80	4	$100\%$	$100\%$	$100\%$	$80.2\%$

Equations

$$\underset{\mathbf{z}\in \mathcal{Z},\mathbf{x}\in {\mathbb{R}}^n}{\min }$$

$$\underset{\mathbf{z}\in \mathcal{Z},\mathbf{x}\in {\mathbb{R}}_{+}^n}{\min }{\mathbf{c}}^{\top }\mathbf{z}+\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Q}\mathbf{x}+{\mathbf{d}}^{\top }\mathbf{x}\text{mbox}s.t.$$

$$x_e=0\text{mbox}if{z}_e=0\forall e\in E,$$

$$g(\mathbf{x})+\Omega (\mathbf{x}):=\{\begin{array}{ll}\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Q}\mathbf{x}+{\mathbf{d}}^{\top }\mathbf{x},& \text{mbox}if\mathbf{A}\mathbf{x}=\mathbf{b},\mathbf{x}\ge 0,\\ +\infty ,& \text{mbox}otherwise.\end{array}$$

$$\underset{\mathbf{z}\in \mathcal{Z}}{\min }$$

$$f(\mathbf{z}):=\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }$$

$$f(\mathbf{z}):=\underset{\mathbf{x}\in {\mathbb{R}}_{+}^n:\mathbf{A}\mathbf{x}=\mathbf{b}}{\min }$$

$$f(\mathbf{z}):=\underset{{\mathbf{f}}^j,\mathbf{x}\in {\mathbb{R}}_{+}^n}{\min }\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Q}\mathbf{x}+{\mathbf{d}}^{\top }\mathbf{x}\text{mbox}s.t.$$

$$\mathbf{x}=\sum _{j=1}^k{\mathbf{f}}^j,\mathbf{x}\le \mathbf{u},$$

$$x_e=0\text{mbox}if{z}_e=0\forall e\in E.$$

$$f(\mathbf{z}):=\underset{\mathbf{w}\in {\mathbb{R}}^p}{\min }$$

$$f(\mathbf{z}):=\underset{\mathbf{x}\in {\mathbb{R}}_{+}^n}{\min }$$

$$\mathbf{l}\le \mathbf{A}\mathbf{x}\le \mathbf{u},{\mathbf{e}}^{\top }\mathbf{x}=1,x_i=0\text{if}{z}_i=0\forall i\in [n].$$

$$\underset{\mathbf{z}}{\min }$$

$$\text{where: }f(\mathbf{z}):=\underset{\mathbf{x}}{\min }$$

$$x_{i,t}\in [0,u_{i,t}]\forall i\in [n],\forall t\in \mathcal{T},$$

$$x_{i,t}=0\text{if}{z}_{i,t}=0\forall i\in [n],\forall t\in \mathcal{T}.$$

$$f(\mathbf{z}):=\underset{\mathbf{X}\in {\mathbb{R}}_{+}^{n\times m}}{\min }$$

$$\sum _{i=1}^n{x}_{ij}=d_j\forall j\in [m],x_{ij}=0\text{ if }{z}_i=0\forall i\in [n],j\in [m].$$

$$f(\mathbf{z})=\underset{\mathbf{x}\in {\mathbb{R}}^p}{\max }{\mathbf{x}}^{\top }\mathbf{\Sigma }\mathbf{x}\text{s.t.}\Vert \mathbf{x}{\Vert }_2^2=1,x_i=0\text{if}{z}_i=0\forall i\in [p].$$

$$f(\mathbf{z})=\underset{\mathbf{X}\in S_{+}^p}{\max }⟨\mathbf{\Sigma },\mathbf{X}⟩\text{s.t.}\mathrm{tr}(\mathbf{X})=1,x_{i,j}=0\text{if}{z}_i=0\text{or}{z}_j=0\forall i,j\in [p].$$

$$f(\mathbf{z})={\mathbf{z}}^{\top }\mathbf{Q}\mathbf{z}.$$

$$f(\mathbf{z}):=\underset{\mathbf{Y}\in {\mathbb{R}}_{+}^{n\times n}}{\min }⟨\mathbf{Q},\mathbf{Y}⟩\text{s.t.}$$

$$y_{i,j}\ge z_i+z_j-1$$

$$y_{i,i}\ge z_i$$

$$y_{i,j}=0\text{mbox}if{z}_i=0$$

$$y_{i,j}=0\text{mbox}if{z}_j=0$$

$$\mathbf{x}=\sum _{i=1}^k{\mathbf{x}}_i,\sum _{i=1}^k{\delta }_i=1,$$

$${\mathbf{A}}_i{\mathbf{x}}_i\le {\mathbf{b}}_i+{\mathbf{\xi }}_i\forall i\in [k],$$

$${\mathbf{x}}_i^{\top }{\mathbf{Q}}_i{\mathbf{x}}_i+{\mathbf{h}}_i^{\top }{\mathbf{x}}_i+g_i\le {\rho }_i\forall i\in [k],$$

$${\mathbf{x}}_i=0\text{if}{\delta }_i=0\forall i\in [k],$$

$${\mathbf{\xi }}_i=0\text{if}(1-{\delta }_i)=0\forall i\in [k],$$

$${\rho }_i=0\text{if}(1-{\delta }_i)=0\forall i\in [k].$$

$$\underset{\mathbf{q}\ge 0}{\max }{\mathbf{e}}^{\top }\mathbf{q}$$

$$\underset{\mathbf{z}\in \mathcal{Z}}{\min }\underset{\mathbf{\alpha }\in {\mathbb{R}}^n}{\max }$$

$${\Omega }^{\star }(\beta )$$

$$f(\mathbf{z})$$

Full text

{AMS}

90C11, 90C57, 90C90

1 Introduction

Many important problems from the Operations Research literature exhibit a logical relationship between continuous variables $x$ and binary variables $z$ of the form “$x=0$ if $z=0$”. Among others, start-up costs in machine scheduling problems, financial transaction costs, cardinality constraints and fixed costs in facility location problems exhibit this relationship. Since the work of [34], this relationship is usually enforced through a “big-$M$” constraint of the form ${-Mz\leq x\leq Mz}$ for a sufficiently large constant $M>0$. Glover’s work has been so influential that big-$M$ constraints are now considered as intrinsic components of the initial problem formulations themselves, to the extent that textbooks in the field introduce facility location, network design or sparse portfolio problems with big-$M$ constraints by default, although they are actually reformulations of logical constraints.

In this work, we adopt a different perspective on the big-$M$ paradigm, viewing it as a regularization term, rather than a modeling trick. Under this lens, we show that regularization drives the computational tractability of problems with logical constraints, explore alternatives to the big-$M$ paradigm and propose an efficient algorithmic strategy which solves a broad class of problems with logical constraints.

1.1 Problem Formulation and Main Contributions

We consider optimization problems which unfold over two stages. In the first stage, a decision-maker activates binary variables, while satisfying resource budget constraints and incurring activation costs. Subsequently, in the second stage, the decision-maker optimizes over the continuous variables. Formally, we consider the problem

[TABLE]

where $\mathcal{Z}\subseteq\{0,1\}^{n}$, $\bm{c}\in\mathbb{R}^{n}$ is a cost vector, $g(\cdot)$ is a generic convex function which possibly models convex constraints $\bm{x}\in\mathcal{X}$ for a convex set $\mathcal{X}\subseteq\mathbb{R}^{n}$ implicitly—by requiring that $g(\bm{x})=+\infty$ if $\bm{x}\notin\mathcal{X}$, and $\Omega(\cdot)$ is a convex regularization function; we formally state its structure in Assumption 2.2.

In this paper, we provide three main contributions: First, we reformulate the logical constraint “$x_{i}=0$ if $z_{i}=0$” in a non-linear way, by substituting $z_{i}x_{i}$ for $x_{i}$ in Problem (1). Second, we leverage the regularization term $\Omega(\bm{x})$ to derive a tractable reformulation of (1). Finally, by invoking strong duality, we reformulate (1) as a mixed-integer saddle-point problem, which is solvable via outer approximation.

Observe that the structure of Problem (1) is quite general, as the feasible set $\mathcal{Z}$ can capture known lower and upper bounds on $\bm{z}$, relationships between different $z_{i}$’s, or a cardinality constraint $\bm{e}^{\top}\bm{z}\leq k$. Moreover, constraints of the form $\bm{x}\in\mathcal{X}$, for some convex set $\mathcal{X}$, can be encoded within the domain of $g$, by defining $g(\bm{x})=+\infty$ if $\bm{x}\notin\mathcal{X}$. As a result, Problem (1) encompasses a large number of problems from the Operations Research literature, such as the network design problem described in Example 1.1. These problems are typically studied separately. However, the techniques developed for each problem are actually different facets of a single unified story, and, as we demonstrate in this paper, can be applied to a much more general class of problems than is often appreciated.

{exbox}

Example 1.1**.**

Network design is an important example of problems of the form (1). Given a set of $m$ nodes, the network design problem consists of constructing edges to minimize the construction plus flow transportation cost. Let $E$ denote the set of all potential edges and let $n=|E|$. Then, the network design problem is given by:

[TABLE]

where $\mathcal{Z}=\{0,1\}^{n}$, $\bm{A}\in\mathbb{R}^{m\times n}$ is the flow conservation matrix, $\bm{b}\in\mathbb{R}^{m}$ is the vector of external demands and $\bm{Q}\in\mathbb{R}^{n\times n}$, $\bm{d}\in\mathbb{R}^{n}$ define the quadratic and linear costs of flow circulation. We assume that $\bm{Q}\succeq\bm{0}$ is a positive semidefinite matrix. Inequalities of the form $\bm{\ell}\leq\bm{z}\leq\bm{u}$ can be incorporated within $\mathcal{Z}$ to account for existing/forbidden edges in the network. Problem (2) is of the same form as Problem (1) with

[TABLE]

*We present a generalized model with edge capacities and multiple commodities in Section 2.1.1. *

1.2 Background and Literature Review

Our work falls into two areas of the mixed-integer optimization literature which are often considered in isolation: $(a)$ modeling forcing constraints which encode whether continuous variables are active and can take non-zero values or are inactive and forced to [math], and $(b)$ decomposition algorithms for mixed-integer optimization problems.

Formulations of forcing constraints

The most popular way to impose forcing constraints on continuous variables is to introduce auxiliary discrete variables which encode whether the continuous variables are active, and relate the discrete and continuous variables via the big-$M$ approach of [34]. This approach was first applied to mixed-integer non-linear optimization (MINLO) in the context of sparse portfolio selection by [14]. With the big-$M$ approach, the original MINLO admits bounded relaxations and can therefore be solved via branch-and-bound. Moreover, because the relationship between discrete and continuous variables is enforced via linear constraints, a big-$M$ reformulation has a theoretically low impact on the tractability of the MINLOs continuous relaxations. However, in practice, high values of $M$ lead to numerical instability and provide low-quality bounds [see 4, Section 5].

This observation led [28] to propose a class of cutting-planes for MINLO problems with indicator variables, called perspective cuts, which often provide a tighter reformulation of the logical constraints. Their approach was subsequently extended by [1], who, building upon the work of [5, pp. 88, item 5], proved that MINLO problems with indicator variables can often be reformulated as mixed-integer second-order cone problems (see [37] for a survey). More recently, a third approach for coupling the discrete and the continuous in MINLO was proposed independently for sparse regression by [47] and [13]: augmenting the objective with a strongly convex term of the form $\|\bm{x}\|_{2}^{2}$, called a ridge regularizer.

In the present paper, we synthesize the aforementioned and seemingly unrelated three lines of research under the unifying lens of regularization. Notably, our framework includes big-$M$ and ridge regularization as special cases, and provides an elementary derivation of perspective cuts.

Numerical algorithms for mixed-integer optimization

A variety of “classical” general-purpose decomposition algorithms have been proposed for general MINLOs. The first such decomposition method is known as Generalized Benders Decomposition, and was proposed by [33] as an extension of [6]. A similar method, known as outer-approximation was proposed by [22], who proved its finite termination. The outer-approximation method was subsequently generalized to account for non-linear integral variables by [25]. These techniques decompose MINLOs into a discrete master problem and a sequence of continuous separation problems, which are iteratively solved to generate valid cuts for the master problem.

Though slow in their original implementation, decomposition schemes have benefited from recent improvements in mixed-integer linear solvers in the past decades, beginning with the branch-and-cut approaches of [45, 48], which embed the cut generation process within a single branch-and-bound tree, rather than building a branch-and-bound tree before generating each cut. We refer to [23, 24] for recent successful implementations of “modern” decomposition schemes. From a high-level perspective, these recent successes require three key ingredients: First, a fast cut generation strategy. Second, as advocated by [23], a rich cut generation process at the root node. Finally, a cut selection rule for degenerate cases where multiple valid inequalities exist (e.g., the Pareto optimality criteria of [43]).

In this paper, we connect the regularization used to reformulate logical constraints with the aforementioned key ingredients for modern decomposition schemes. Hence, instead of considering a MINLO formulation as a given and subsequently attempt to solve it at scale, our approach view big-$M$ constraints as one of many alternatives. We argue that regularization is a modeling choice that impacts the tractability of the formulation and should be made accordingly.

1.3 Structure

We propose a unifying framework to address mixed-integer optimization problems, and jointly discuss modeling choice and numerical algorithms.

In Section 2, we identify a general class of mixed-integer optimization problems, which encompasses sparse regression, sparse portfolio selection, sparse principal component analysis, unit commitment, facility location, network design and binary quadratic optimization as special cases. For this class of problems, we discuss how imposing either big-$M$ or ridge regularization accounts for non-linear relationships between continuous and binary variables in a tractable fashion. We also establish that regularization controls the convexity and smoothness of Problem (1)’s objective function.

In Section 3, we propose a conjunction of general-purpose numerical algorithms to solve Problem (1). The backbone of our approach is an outer approximation framework, enhanced with first-order methods to solve the Boolean relaxations and obtain improved lower bounds, certifiably near-optimal warm-starts via randomized rounding, and a discrete local search procedure. We also connect our approach to the perspective cut approach [28] from a theoretical and implementation standpoint.

Finally, in Section 4, we demonstrate empirically that algorithms derived from our framework can outperform state-of-the-art solvers. On network design problems with $100$s of nodes and binary quadratic optimization problems with $100$s of variables, we improve the objective value of the returned solution by $5$ to $40\%$ and $5$ to $85\%$ respectively, and our edge increases as the problem size increases. On empirical risk minimization problems, our method with ridge regularization is able to accurately select features among $100,000$s (resp. $10,000$s) of covariates for regression (resp. classification) problems, with higher accuracy than both Lasso and non-convex penalties from the statistics literature. For sparse portfolio selection, we solve to provable optimality problems one order of magnitude larger than previous attempts. We then analyze the benefits of the different ingredients in our numerical recipe on facility location problems, and discuss the relative merits of different regularization approaches on unit commitment instances.

Notation

We use nonbold face characters to denote scalars and components of matrices, lowercase bold faced characters such as $\bm{x}$ to denote vectors, uppercase bold faced characters such as $\bm{X}$ to denote matrices, and calligraphic characters such as $\mathcal{X}$ to denote sets. We let $\mathbf{e}$ denote a vector of all $1$’s, and $\bm{0}$ denote a vector of all [math]’s, with dimension implied by the context. If $\bm{x}$ is a $n$-dimensional vector then $\mathrm{Diag}(\bm{x})$ denotes the $n\times n$ diagonal matrix whose diagonal entries are given by $\bm{x}$. If $f(\bm{x})$ is a convex function then its perspective function $\varphi(\bm{x},t)$, defined as $\varphi(\bm{x},t)=tf(\bm{x}/t)$ if $t>0$, $\varphi(\mathbf{0},0)=0$, and $\infty$ elsewhere, is also convex [17, Chapter 3.2.6.]. Finally, we let $\mathbb{R}_{+}^{n}$ denote the $n$-dimensional nonnegative orthant.

2 Framework and Examples

In this section, we present the family of problems to which our analysis applies, discuss the role played by regularization, and provide some examples from the Operations Research literature.

2.1 Examples

Problem (1) has a two-stage structure which comprises first “turning on” some indicator variables $\bm{z}$, and second solving a continuous optimization problem over the active components of $\bm{x}$. Precisely, Problem (1) can be viewed as a discrete optimization problem:

[TABLE]

where the inner minimization problem

[TABLE]

yields a best choice of $\bm{x}$ given $\bm{z}$. As we illustrate in this section, a number of problems of practical interest exhibit this structure.{exbox}

Example 2.1**.**

For the network design example (2), we have

[TABLE]

2.1.1 Network Design

Example 1.1 illustrates that the single-commodity network design problem is a special case of Problem (1). We now formulate the $k$-commodity network design problem with directed capacities as minimizing over $\mathcal{Z}=\{0,1\}^{n}$ the function:

[TABLE]

2.1.2 Sparse Empirical Risk Minimization

Given a matrix of covariates $\bm{X}\in\mathbb{R}^{n\times p}$ and a response vector $\bm{y}\in\mathbb{R}^{n}$, the sparse empirical risk minimization problem seeks a vector $\bm{w}$ which explains the response in a compelling manner, i.e., minimizes over $\mathcal{Z}:=\{\bm{z}\in\{0,1\}^{p}:\bm{e}^{\top}\bm{z}\leq k\}$ the function:

[TABLE]

where $\ell$ is an appropriate convex loss function; we provide examples of suitable loss functions in Table 1.

2.1.3 Sparse Portfolio Selection

Given an expected marginal return vector $\bm{\mu}\in\mathbb{R}^{n}$, estimated covariance matrix ${\bm{\Sigma}\in\mathcal{S}^{n}_{+}}$, uncertainty budget parameter $\sigma>0$, cardinality budget parameter $k\in\{2,\ldots,n-1\}$, linear constraint matrix $\bm{A}\in\mathbb{R}^{n\times m}$, and right-hand-side bounds $\bm{l},\bm{u}\in\mathbb{R}^{m}$, investors determine an optimal allocation of capital between assets by minimizing over $\mathcal{Z}=\left\{\bm{z}\in\{0,1\}^{n}:\ \bm{e}^{\top}\bm{z}\leq k\right\}$ the function

[TABLE]

2.1.4 Unit Commitment

In the DC-load-flow unit commitment problem, each generation unit $i$ incurs a cost given by a quadratic cost function $f^{i}(x)=a_{i}x^{2}+b_{i}x+c_{i}$ for its power generation output $x\in[0,u_{i}]$. Let $\mathcal{T}$ denote a finite set of time periods covering a time horizon (e.g., $24$ hours). At each time period $t\in\mathcal{T}$, there is an estimated demand $d_{t}$. The objective is to generate sufficient power to satisfy demand at minimum cost, while respecting minimum time on/time off constraints.

By introducing binary variables $z_{i,t}$, which denote whether generation unit $i$ is active in time period $t$, requiring that $\bm{z}\in\mathcal{Z}$, i.e., $\bm{z}$ obeys physical constraints such as minimum time on/off, the unit commitment problem admits the formulation:

[TABLE]

[TABLE]

2.1.5 Facility Location

Given a set of $n$ facilities and $m$ customers, the facility location problem consists of constructing facilities $i\in[n]$ at cost $c_{i}$ to satisfy demand at minimal cost, i.e., minimizing over $\mathcal{Z}=\{0,1\}^{n}$ the function:

[TABLE]

In this formulation, ${x}_{ij}$ corresponds to the quantity produced in facility $i$ and shipped to customer $j$ at a marginal cost of ${c}_{ij}$. Moreover, each facility $i$ has a maximum output capacity of ${u}_{i}$ and each customer $j$ has a demand of $d_{j}$. In the uncapacitated case where ${u}_{i}=\infty$, the inner minimization problems decouple into independent knapsack problems for each customer $j$.

2.1.6 Sparse Principal Component Analysis (PCA)

Given a $p\times p$ positive semidefinite covariance matrix $\bm{\Sigma}$, $\bm{\Sigma}\in S^{p}_{+}$ in short, the sparse PCA problem is to select a vector $\bm{z}$ which maximizes over $\mathcal{Z}=\left\{\bm{z}\in\{0,1\}^{p}:\ \bm{e}^{\top}\bm{z}\leq k\right\}$ the function

[TABLE]

This function is apparently non-concave in $\bm{z}$, because $f(\bm{z})$ is the optimal value of a non-convex quadratic optimization problem. Fortuitously however, this problem admits an exact mixed-integer semidefinite reformulation, namely

[TABLE]

Indeed, for any fixed $\bm{z}$, Problem (12) maximizes a linear function in $\bm{X}$ and therefore admits a rank-one optimal solution. Thus, we prove that sparse PCA admits an exact mixed-integer semidefinite optimization reformulation.

2.1.7 Binary Quadratic Optimization

Given a symmetric cost matrix $\bm{Q}$, the binary quadratic optimization problem consists of selecting a vector of binary variables $\bm{z}$ which minimizes over $\mathcal{Z}=\{0,1\}^{n}$ the function:

[TABLE]

This formulation is non-convex and does not include continuous variables. However, introducing auxiliary continuous variables yields the equivalent formulation [26] of minimizing over $\mathcal{Z}=\{0,1\}^{n}$ the function:

[TABLE]

2.1.8 Union of Ellipsoidal Constraints

We now demonstrate that an even broader class of problems than MIOs with logical constraints can be cast within our framework. Concretely, we demonstrate that constraints $\bm{x}\in\mathcal{S}:=\bigcup_{i=1}^{k}(Q_{i}\cap P_{i}),$ where $Q_{i}:=\{\bm{x}\in\mathbb{R}^{n}:\bm{x}^{\top}\bm{Q}_{i}\bm{x}+\bm{h}_{i}^{\top}\bm{x}+g_{i}\leq 0\}$, with $\bm{Q}_{i}\succeq\bm{0}$, is an ellipsoid and $P_{i}:=\{\bm{x}:\bm{A}_{i}\bm{x}\leq\bm{b}_{i}\}$ is a polytope, can be reformulated as a special case of our framework. We remark that the constraint $\bm{x}\in\mathcal{S}$ is very general. Indeed, if we were to omit the quadratic constraints then we obtain a so-called ideal union of polyhedra formulation, which essentially all mixed-binary linear feasible regions admit [see 52].

To derive a mixed-integer formulation with logical constraints of $\mathcal{S}$ that fits within our framework, we introduce $\bm{x}_{i}\in\mathbb{R}^{n}$ and $\delta_{i}\in\{0,1\}^{n}$, such that $\bm{x}_{i}\in Q_{i}\cap P_{i}$ if $\delta_{i}=1$, $\bm{x}_{i}=\bm{0}$ otherwise, and $\bm{x}=\sum_{i}\bm{x}_{i}$. We enforce $\bm{x}_{i}\in Q_{i}\cap P_{i}$ by introducing slack variables $\bm{\xi_{i}}$, $\rho_{i}$ for the linear and quadratic constraints respectively, and forcing them to be zero whenever $\delta_{i}=1$. Formally, $\mathcal{S}$ admits the following formulation

[TABLE]

2.2 A Regularization Assumption

When we stated Problem (1), we assumed that its objective function consists of a convex function $g(\bm{x})$ plus a regularization term $\Omega(\bm{x})$. We now formalize this assumption: {assumption} In Problem (1), the regularization term $\Omega(\bm{x})$ is one of:

• •

a big-$M$ penalty function, $\Omega(\bm{x})=0$ if $\ \|\bm{x}\|_{\infty}\leq M$ and $\infty$ otherwise,

• •

a ridge penalty, $\Omega(\bm{x})=\dfrac{1}{2\gamma}\|\bm{x}\|_{2}^{2}$.

This decomposition often constitutes a modeling choice in itself. We now illustrate this idea via the network design example.{exbox}

Example 2.2**.**

In the network design example (2), given the flow conservation structure $\bm{A}\bm{x}=\bm{b}$, we have that $\bm{x}\leq M\bm{e}$, where $M=\sum_{i:b_{i}>0}b_{i}$. In addition, if $\bm{Q}\succ\mathbf{0}$ then the objective function naturally contains a ridge regularization term with $1/\gamma$ equal to the smallest eigenvalue of $\bm{Q}$. Moreover, it is possible to obtain a tighter natural ridge regularization term by solving the following auxiliary semidefinite optimization problem a priori

[TABLE]

*and using $q_{i}$ as the ridge regularizer for each index $i$ [30]. *

Big-$M$ constraints are often considered to be a modeling trick. However, our framework demonstrates that imposing either big-$M$ constraints or a ridge penalty is a regularization method, rather than a modeling trick. Interestingly, ridge regularization accounts for the relationship between the binary and continuous variables just as well as big-$M$ regularization, without performing an algebraic reformulation of the logical constraints111Specifically, ridge regularization enforces logical constraints through perspective functions, as is made clear in Section 3.4..

Conceptually, both regularization functions are equivalent to a soft or hard constraint on the continuous variables $\bm{x}$. However, they admit practical differences: For big-$M$ regularization, there usually exists a finite value ${M}_{0}$, typically unknown a priori, such that if $M<{M}_{0}$, the regularized problem is infeasible. Alternatively, for every value of the ridge regularization parameter $\gamma$, if the original problem is feasible then the regularized problem is also feasible. Consequently, if there is no natural choice of $M$ then imposing ridge regularization may be less restrictive than imposing big-$M$ regularization. However, for any $\gamma>0$, the objective of the optimization problem with ridge regularization is different from its unregularized limit as $\gamma\rightarrow\infty$, while for big-$M$ regularization, there usually exists a finite value $M_{1}$ above which the two objective values match. We illustrate this discussion numerically in Section 4.3.

2.3 Duality to the Rescue

In this section, we derive Problem (4)’s dual and reformulate $f(\bm{z})$ as a maximization problem. This reformulation is significant for two reasons: First, as shown in the proof of Theorem 2.3, it leverages a non-linear reformulation of the logical constraints “$x_{i}=0$ if $z_{i}=0$” by introducing additional variables $v_{i}$ such that $v_{i}=z_{i}x_{i}$. Second, it proves that the regularization term $\Omega(\bm{x})$ drives the convexity and smoothness of $f(\bm{z})$, and thereby drives the computational tractability of the problem. To derive Problem (4)’s dual, we require:{assumption} For each subproblem generated by $f(\bm{z})$, where $\bm{z}\in\mathcal{Z}$, either the optimization problem is infeasible, or strong duality holds.

Note that all seven problems stated in Section 2.1 satisfy Assumption 2.3, as their inner problems are convex quadratics with linear or semidefinite constraints [17, Section 5.2.3]. Under Assumption 2.3, the following theorem reformulates Problem (3) as a saddle-point problem:

Theorem 2.3**.**

Under Assumption 2.3, Problem (3) is equivalent to:

[TABLE]

where $h(\bm{\alpha}):=\inf_{\bm{v}}\>g(\bm{v})-\bm{v}^{\top}\bm{\alpha}$ is, up to a sign, the Fenchel conjugate of $g$ [see 17, Chap. 3.3], and

[TABLE]

Proof 2.4**.**

Let us fix some $\bm{z}\in\{0,1\}^{n}$, and suppose that strong duality holds for the inner minimization problem which defines $f(\bm{z})$. Then, after introducing additional variables $\bm{v}\in\mathbb{R}^{n}$ such that $v_{i}=z_{i}x_{i}$, we have

[TABLE]

Let $\bm{\alpha}$ denote the dual variables associated with the coupling constraint $\bm{v}=\mathrm{Diag}(\bm{z})\bm{x}$. The minimization problem is then equivalent to its dual problem, which is given by:

[TABLE]

Since $\Omega(\cdot)$ is decomposable, i.e., $\Omega(\bm{x})=\sum_{i}\Omega_{i}({x}_{i})$, we obtain:

[TABLE]

where the last equality holds as ${z}_{i}>0$ for the big-$M$ and ${z}_{i}^{2}={z}_{i}$ for the ridge penalty.

*Alternatively, if the inner minimization problem defining $f(\bm{z})$ is infeasible, then its dual problem is unbounded by weak duality222Weak duality implies that the dual problem is either unfeasible or unbounded. Since the feasible set of the maximization problem does not depend on $\bm{z}$, it is always feasible, unless the original problem (1) is itself infeasible. Therefore, we assume without loss of generality that it is unbounded.. *

Remark 2.5**.**

Without regularization, i.e., $\Omega(\bm{x})=0$, a similar proof shows that Problem (3) admits an interesting saddle-point formulation:

[TABLE]

*since $\Omega^{\star}(\alpha)={\min_{{x}}\left[{x}{\alpha}-\Omega({x})\right]=}0$ if $\alpha=0$, and $+\infty$ otherwise. Consequently, the regularized formulation can be regarded as a relaxation of the original problem where the hard constraint $\alpha_{i}=0\ \text{if}\ z_{i}=1$ is replaced with a soft penalty term $-z_{i}\Omega^{\star}(\alpha_{i})$. *

Remark 2.6**.**

*The proof of Theorem 2.3 exploits three attributes of the regularizer $\Omega(\bm{x})$. Namely, (1) decomposability, i.e., $\Omega(\bm{x})=\sum_{i}\Omega_{i}(x_{i})$, for appropriate scalar functions $\Omega_{i}$, (2) the convexity of $\Omega(\bm{x})$ in $\bm{x}$, and (3) the fact that $\Omega(\cdot)$ regularizes333Importantly, the third attribute allows us to strengthen the formulation by not associating $\bm{z}$ with $\bm{x}$ in $\Omega(\bm{x})$, since $\bm{x}_{i}=0$ is a feasible, indeed optimal choice of $\bm{x}$ for minimizing the regularizer when $z_{i}=0$; this issue is explored in more detail in [8, Lemma 1]; see also [11, Appendix A.1]. towards [math], i.e., $\bm{0}\in\arg\min_{\bm{x}}\Omega(\bm{x})$. However, the proof does not explicitly require that $\Omega(\bm{x})$ is either a big-$M$ or a ridge regularizer. This suggests that our framework could be extended to other regularization functions. *

{exbox}

Example 2.7**.**

For the network design problem (2), we have

[TABLE]

Introducing $\bm{\xi}=\bm{Q}^{-1/2}\left(\bm{A}^{\top}\bm{p}-\bm{d}+\bm{\alpha}+\bm{\beta}_{0}\right)$, we can further write

[TABLE]

Hence, Problem (2) is equivalent to minimizing over $\bm{z}\in\mathcal{Z}$ the function

[TABLE]

Theorem 2.3 reformulates $f(\bm{z})$ as an inner maximization problem, namely

[TABLE]

for any feasible binary $\bm{z}\in\mathcal{Z}$. The regularization term $\Omega$ will be instrumental in our numerical strategy for it directly controls both the convexity and smoothness of $f$. Note that (17) extends the definition of $f(\bm{z})$ to the convex set $\mathrm{Bool}(\mathcal{Z})$, obtained by relaxing the constraints $\bm{z}\in\{0,1\}^{p}$ to $\bm{z}\in[0,1]^{p}$ in the definition of $\mathcal{Z}$.

Convexity

$f(\bm{z})$* is convex in $\bm{z}$ as a point-wise maximum of linear function of $\bm{z}$. In addition, denoting $\bm{\alpha}^{\star}(\bm{z})$ a solution of (17), we have the lower-approximation:*

[TABLE]

where $[\nabla f({\bm{z}})]_{i}:=-\Omega^{\star}(\alpha^{\star}(\bm{z})_{i})$ is a sub-gradient of $f$ at $\bm{z}$.

We remark that if the maximization problem in $\bm{\alpha}$ defined by $f(\bm{z})$ admits multiple optimal solutions then the corresponding lower-approximation of $f$ at $\bm{z}$ may not be unique. This behavior can severely hinder the convergence of outer-approximation schemes such as Benders’ decomposition. Since the work of **[43]** on Pareto optimal cuts, many strategies have been proposed to improve the cut selection process in the presence of degeneracy [see 23, Section 4.4 for a review]. However, the use of ridge regularization ensures that the objective function in (16) is strongly concave in ${\alpha}_{i}$ such that $z_{i}>0$, and therefore guarantees that there is a unique optimal choice of ${\alpha}_{i}^{\star}(\bm{z})$. In other words, ridge regularization naturally inhibits degeneracy.

Smoothness

$f(\bm{z})$* is smooth, in the sense of Lipschitz continuity, which is a crucial property for deriving bounds on the integrality gap of the Boolean relaxation, and designing local search heuristics in Section 3. Formally, the following proposition follows from Theorem 2.3:*

Proposition 2.8**.**

For any $\bm{z},\bm{z}^{\prime}\in\mathrm{Bool}\left(\mathcal{Z}\right)$,

1. (a)

With big-$M$* regularization, $f(\bm{z}^{\prime})-f(\bm{z})\leq M\displaystyle\sum_{i=1}^{n}(z_{i}-z^{\prime}_{i})|\alpha^{\star}(\bm{z}^{\prime})_{i}|$.* 2. (b)

With ridge regularization, $f(\bm{z}^{\prime})-f(\bm{z})\leq\dfrac{\gamma}{2}\displaystyle\sum_{i=1}^{n}(z_{i}-z^{\prime}_{i})\alpha^{\star}(\bm{z}^{\prime})_{i}^{2}$.

Proof 2.9**.**

By Equation (16),

[TABLE]

*where the inequality holds because an optimal choice of $\bm{\alpha}^{\prime}$ is a feasible choice of $\bm{\alpha}$. *

Proposition 2.8 demonstrates that, when the coordinates of $\bm{\alpha}^{\star}(\bm{z})$ are uniformly bounded444Such a uniform bound always exists, as $f(\bm{z})$ is only supported on a finite number of binary points. Moreover, the strong concavity of $h$ can yield stronger bounds (see Appendix B).* with respect to $\bm{z}$, $f(\bm{z})$ is Lipschitz-continuous, with a constant $L$ proportional to $M$ (resp. $\gamma$) in the big-$M$ (resp. ridge) case. We provide explicit bounds on the magnitude of $L$ in Appendix B.*

2.4 Merits of Ridge, Big-$M$ Regularization: Theoretical Perspective

In this section, we propose a framework to reformulate MINLOs with logical constraints, which comprises regularizing MINLOs via either the widely used big-$M$* modeling paradigm or the less popular ridge regularization paradigm. We summarize the advantages and disadvantages of each regularizer in Table 2. However, note that we have not yet established how these characteristics impact the numerical tractability and quality of the returned solution; this is the topic of the next two sections.*

3 An Efficient Numerical Approach

We now present an efficient numerical approach to solve Problem (16). The backbone is an outer-approximation strategy, embedded within a branch-and-bound procedure to solve the problem exactly. We also propose local search and rounding heuristics to find good feasible solutions, and use information from the Boolean relaxation to improve the duality gap.

3.1 Overall Outer-Approximation Scheme

Theorem 2.3 reformulates the function $f(\bm{z})$ as an inner maximization problem, and demonstrates that $f(\bm{z})$ is convex in $\bm{z}$, meaning a linear outer approximation provides a valid underestimator of $f(\bm{z})$, as outlined in Equation (18). Consequently, a valid numerical strategy for minimizing $f(\bm{z})$ is to iteratively minimize a piecewise linear lower-approximation of $f$ and refining this approximation at each step until some approximation error $\varepsilon$ is reached, as described in Algorithm 1. This scheme was originally proposed for continuous decision variables by **[40]**, and later extended to binary decision variables by **[22]**, who provide a proof of termination in a finite, yet exponential in the worst case, number of iterations.

To avoid solving a mixed-integer linear optimization problem at each iteration, as suggested in the pseudo-code, this strategy can be integrated within a single branch-and-bound procedure using lazy callbacks, as originally proposed by **[48]**. Lazy callbacks are now standard tools in commercial solvers such as Gurobi and CPLEX and provide significant speed-ups for outer-approximation algorithms. With this implementation, the commercial solver constructs a single branch-and-bound tree and generates a new cut at a feasible solution $\bm{z}$.

We remark that the second-stage minimization problem may be infeasible at some $\bm{z}^{t}$. In this case, we generate a feasibility cut rather than outer-approximation cut. In particular, the constraint $\sum_{i}z^{t}_{i}(1-z_{i})+\sum_{i}(1-z^{t}_{i})z_{i}\geq 1$ excludes the iterate $\bm{z}^{t}$ from the feasible set. Stronger feasibility cuts can be obtained by leveraging problem specific structure. For instance, when the feasible set satisfies $\bm{z}^{t}\notin\mathcal{Z}\implies\forall\bm{z}\leq\bm{z}^{t},\ \bm{z}\notin\mathcal{Z}$, $\sum_{i}(1-z^{t}_{i})z_{i}\geq 1$ is a valid feasibility cut. Alternatively, one can invoke conic duality if $g(\bm{x})$ generates a conic feasibility problem. Formally, assume

[TABLE]

where $\mathcal{K}$ is a closed convex cone. This assumption gives rise to some loss of generality. Note, however, that all the examples in the previous section admit conic reformulations by taking appropriate Cartesian products of the linear, second-order and semidefinite cones **[5]**. Assuming that $g(\bm{x})$ is of the prescribed form, we have the dual conjugate

[TABLE]

where $\mathcal{K}^{\star}$ is the dual cone to $\mathcal{K}$. In this case, if some binary vector $\bm{z}$ gives rise to an infeasible subproblem, i.e., $f(\bm{z})=+\infty$, then the conic duality theorem implies555We should note that this statement is, strictly speaking, not true unless we impose regularization. Indeed, the full conic duality theorem [5, Theorem 2.4.1] allows for the possibility that a problem is infeasible but asymptotically feasible, i.e.,

$\nexists\bm{x}:\bm{A}\bm{x}=\bm{b},\bm{x}\in\mathcal{K}\ \text{but}\ \exists\{\bm{x}_{t}\}_{t=1}^{\infty}:\bm{x}_{t}\in\mathcal{K}\ \forall t\ \text{with}\ \|\bm{A}\bm{x}_{t}-\bm{b}\|\rightarrow 0.$

Fortunately, the regularizer $\Omega(\bm{x})$ alleviates this issue, because it is coercive (i.e., “blows up” to $+\infty$ as $\|\bm{x}\|\rightarrow\infty$) and therefore renders all unbounded solutions infeasible and ensures the compactness of the level sets of $g(\bm{x})+\Omega(\bm{x})$.* that there is a certificate of infeasibility $(\bm{\alpha},\bm{\pi}$) such that*

[TABLE]

Therefore, to restore feasibility, we can simply impose the cut $\langle\bm{b},\bm{\pi}\rangle\leq\sum_{i=1}^{n}z_{i}\Omega^{\star}(\alpha_{i}).$

As mentioned in Section 1.2, the rate of convergence of outer-approximation schemes depends heavily on three criterion. We now provide practical guidelines on how to meet these criterion:

Fast cut generation strategy: ** To generate a cut, one solves the second-stage minimization problem (4) (or its dual) in $\bm{x}$, which contains no discrete variables and is usually orders of magnitude faster to solve than the original mixed-integer problem (1). Moreover, the minimization problem in $\bm{x}$ needs to be solved only for the coordinates $x_{i}$ such that $z_{i}=1$. In practice, this approach yields a sequence of subproblems of much smaller size than the original problem, especially if $\mathcal{Z}$ contains a cardinality constraint. For instance, for the sparse empirical risk minimization problem (6), each cut is generated by solving a subproblem with $n$ observations and $k$ features, where $k\ll p$. For this reason, we recommend generating cuts at binary $\bm{z}$’s, which are often sparser than continuous $\bm{z}$’s. This recommendation can be relaxed in cases where the separation problem can be solved efficiently even for dense $\bm{z}$’s; for instance, in uncapacitated facility location problems, each subproblem is a knapsack problem which can be solved by sorting **[24]**. If possible, we recommend theoretically analyzing the sparsity of the optimal solution a priori, to derive an explicit cardinality or budget constraint on $\bm{z}$ and ensure the sparsity of each incumbent solution. 2. 2.

Cut selection rule in presence of degeneracy: ** In the presence of degeneracy, selection criteria, such as Pareto optimality **[43]**, have been proposed to accelerate convergence. However, these criteria are numerous, computationally expensive and all in all, can do more harm than good **[46]**. In an opposite direction, we recommend alleviating the burden of degeneracy by design, by imposing a ridge regularizer whenever degeneracy hinders convergence. 3. 3.

Rich root node analysis: ** As suggested in **[23]**, providing the solver with as much information as possible at the root node can drastically improve convergence of outer-approximation methods. This is the topic of the next two sections. Restarting mechanisms, as described in **[23, Section 5.2]**, could also be useful, although we do not implement them in the present paper.

These ingredients, and especially the ability to generate cuts efficiently, dictate which types of problems could benefit the most from our approach and which regularizer to use. Problems with an explicit cardinality constraint, for instance, would require a small subproblem to be solved at each iteration. For network design problems, the network flow structure of the feasible set is a key numerical asset so we intuit that ridge regularization, which leaves the feasible set unchanged, would be very efficient. On the other hand, for uncapacitated facility location, sub-problems with big-$M$* regularization boils down to a knapsack problem and can be solved efficiently via sorting, as discussed in [24, Section 3.1].*

3.2 Improving the Lower-Bound: A Boolean Relaxation

To certify optimality, high-quality lower bounds are of interest and can be obtained by relaxing the integrality constraint $\bm{z}\in\{0,1\}^{n}$ in the definition of $\mathcal{Z}$ to $\bm{z}\in[0,1]^{n}$. In this case, the Boolean relaxation of (3) is:

[TABLE]

which can be solved using Kelley’s algorithm **[40]**, which is a continuous analog of Algorithm 1. Stabilization strategies have been empirically successful to accelerate the convergence of Kelley’s algorithm, as recently demonstrated on uncapacitated facility location problems by **[24]**. However, for Boolean relaxations, Kelley’s algorithm computes $f(\bm{z})$ and $\nabla f(\bm{z})$ at dense vectors $\bm{z}$, which is (sometimes substantially) more expensive than for sparse binary vectors $\bm{z}$’s, unless each subproblem can be solved efficiently as in **[24]**.

Alternatively, the continuous minimization problem admits a reformulation

[TABLE]

analogous to Problem (16). Under Assumption 2.3, we can further write the min-max relaxation formulation (19) as a non-smooth maximization problem

[TABLE]

and apply a projected sub-gradient ascent method as in **[12]**. We refer to **[7, Chapter 7.5.]** for a discussion on implementation choices regarding step-size schedule and stopping criteria, and **[50]** for recent enhancements using restarting.

The benefit from solving the Boolean relaxation with these algorithms is threefold. First, it provides a lower bound on the objective value of the discrete optimization problem (3). Second, it generates valid linear lower approximations of $f(\bm{z})$ to initiate the cutting-plane algorithm with. Finally, it supplies a sequence of continuous solutions that can be rounded and polished to obtain good binary solutions. Indeed, the Lipschitz continuity of $f(\bm{z})$ suggests that high-quality feasible binary solutions can be found in the neighborhood of a solution to the Boolean relaxation. We formalize this observation in the following theorem:

Theorem 3.1**.**

Let $\bm{z}^{\star}$ denote a solution to the Boolean relaxation (19), $\mathcal{R}$ denote the indices of $\bm{z}^{\star}$ with fractional entries, and $\alpha^{\star}(\bm{z})$ denote a best choice of $\bm{\alpha}$ for a given $\bm{z}$. Suppose that for any $\bm{z}\in\mathcal{Z}$, $|\alpha^{\star}(\bm{z})_{j}|\leq L$. Then, a random rounding $\bm{z}$ of $\bm{z}^{\star}$, i.e., $z_{j}\sim Bernoulli(z^{\star}_{j})$, satisfies $0\leq f(\bm{z})-f(\bm{z}^{\star})\leq\epsilon$ with probability at least $p=1-|\mathcal{R}|\exp\left(-\tfrac{\epsilon^{2}}{\kappa}\right),$ where

[TABLE]

We provide a formal proof of this result in Appendix A.1. This result calls for multiple remarks:

• •

For $\varepsilon>\sqrt{\kappa\ln(|\mathcal{R}|)}$, we have that $p>0$, which implies the existence of a binary $\varepsilon$-optimal solution in the neighborhood of $\bm{z}^{\star}$, which in turn bounds the integrality gap by $\varepsilon$. As a result, lower values of $M$ or $\gamma$ typically make the discrete optimization problem easier.

• •

A solution to the Boolean relaxation often includes some binary coordinates, i.e., $|\mathcal{R}|<n$. In this situation, it is tempting to fix $z_{i}=z^{\star}_{i}$ for $i\notin\mathcal{R}$ and solve the master problem (3) over coordinates in $\mathcal{R}$. In general, this approach provides sub-optimal solutions. However, Theorem 3.1 quantifies the price of fixing variables and bounds the optimality gap by $\sqrt{\kappa\ln(|\mathcal{R}|)}$.

• •

In the above high-probability bound, we do not account for the feasibility of the randomly rounded solution $\bm{z}$. Accounting for $\bm{z}$’s feasibility marginally reduces the probability given above, as shown for general discrete optimization problems by **[49]**.

Under specific problem structure, other strategies might be more efficient than Kelley’s method or the subgradient algorithm. For instance, if $\mathrm{Bool}\left(\mathcal{Z}\right)$ is a polyhedron, then the inner minimization problem defining $q(\bm{\alpha})$ is a linear optimization problem that can be rewritten as a maximization problem by invoking strong duality. Although we only consider linear relaxations here, tighter bounds could be attained by taking a higher-level relaxation from a relaxation hierarchy, such as the **[41]** hierarchy [see 42, for a comparison]. The main benefit of such a relaxation is that while the aforementioned Boolean relaxation only controls the first moment of the probability measure studied in Theorem 3.1, higher level relaxations control an increasing sequence of moments of the probability measure and thereby provide non-worsening probabilistic guarantees for randomized rounding methods. However, the additional tightness of these bounds comes at the expense of solving relaxations with additional variables and constraints666$n^{2}$ additional variables and $n^{2}$ additional constraints for empirical risk minimization, versus $n+1$ additional variables and $n$ additional constraints for the linear relaxation.; yielding a sequence of ever-larger semidefinite optimization problems. Indeed, even the SDP relaxation which controls the first two moments of a randomized rounding method is usually intractable when $n>300$, with current technology. For an analysis of higher-level relaxations in sparse regression problems, we refer the reader to **[2]**.

3.3 Improving the Upper-Bound: Local Search and Rounding

To improve the quality of the upper-bound, i.e., the cost associated with the best feasible solution found so far, we implement two rounding and local-search strategies.

Our first strategy is a randomized rounding strategy, which is inspired by Theorem 3.1. Given $\bm{z}_{0}\in\mathrm{Bool}(\mathcal{Z})$, we generate randomly rounded vectors $\bm{z}$ by sampling $\bm{z}$ according to $z_{i}\sim\mathrm{Bernoulli}(z_{0i})$ until $\bm{z}\in\mathcal{Z}$, which happens with high probability since $\mathbb{E}[\bm{z}]=\bm{z}_{0}$ satisfies all the constraints which describe $\mathcal{Z}$, besides integrality **[49]**.

*Our second strategy is a sequential rounding procedure, which is informed by the lower-approximation on $f(\bm{z})$, as laid out in Equation (18). Observing that the *$i$th coordinate $\nabla f(\bm{z}_{0})_{i}$ provides a first-order indication of how a change in $z_{i}$ might impact the overall cost, we proceed in two steps. We first round down all coordinates such that $\nabla f(\bm{z}_{0})_{i}(0-\bm{z}_{0i})<0$. Once the linear approximation of $f$ only suggests rounding up, we round all coordinates of $\bm{z}$ to $1$ and iteratively bring some coordinates to [math] to restore feasibility.

*If $\bm{z}_{0}$ is binary, we implement a comparable local search strategy. If $\bm{z}_{0i}=0$, then switching the *$i$th coordinate to one increases the cost by at least $\nabla f(\bm{z}_{0})_{i}$. Alternatively, if $\bm{z}_{0i}=1$, then switching it to zero increases the cost by at least $-\nabla f(\bm{z}_{0})_{i}$. We therefore compute the one-coordinate change which provides the largest potential cost improvement. However, as we only have access to a lower approximation of $f$, we are not guaranteed to generate a cost-decreasing sequence. Therefore, we terminate the procedure as soon as it cycles. A second complication is that, due to the constraints defining $\mathcal{Z}$, the best change sometimes yields an infeasible $\bm{z}$. In practice, for simple constraints such as $\bm{\ell}\leq\bm{z}\leq\bm{u}$, we forbid switches which break feasibility; for cardinality constraints, we perform the best switch and then restore feasibility at minimal cost when necessary.

3.4 Relationship With Perspective Cuts

In this section, we connect the perspective cuts introduced by **[28]** with our framework and discuss the merits of both approaches, in theory and in practice. To the best of our knowledge, a connection between Boolean relaxations of the two approaches has only been made in the context of sparse regression, by **[54]**. That is, the general connection we make here between the discrete problems, as well as their respective cut generating procedures, is novel.

We first demonstrate that imposing the ridge regularization term $\Omega(\bm{x})=\tfrac{1}{2\gamma}\|\bm{x}\|_{2}^{2}$ naturally leads to the perspective formulation of **[28]**:

Theorem 3.2**.**

Suppose that $\Omega(\bm{x})=\tfrac{1}{2\gamma}\|\bm{x}\|_{2}^{2}$ and that Assumption 2.3 holds. Then, Problem (16) is equivalent to the following optimization problem:

[TABLE]

Theorem 3.2 follows from taking the dual of the inner-maximization problem in Problem (17); see Appendix A.2 for a formal proof. Note that the equivalence stated in Theorem 3.2 also holds for $\bm{z}\in\mathrm{Bool}(\mathcal{Z})$. As previously observed in **[5, 1]**, Problem (20) can be formulated as a second-order cone problem (SOCP)

[TABLE]

and solved by linearizing the SOCP constraints into so-called perspective cuts, i.e., ${\theta_{i}\geq\tfrac{1}{2\gamma}\bar{x}_{i}(2x_{i}-\bar{x}_{i}z_{i}),\forall\bar{x}\in\bar{\mathcal{X}}}$, which have been extensively studied in the literature in the past fifteen years **[28, 36, 20, 27, 2]**. Observe that by separating Problem (20) into master and subproblems, an outer approximation algorithm yields the same cut (18) as in our scheme. In this regard, our approach supplies a new and insightful derivation of the perspective cut approach. It is worth noting that our proposal can easily be implemented within a standard integer optimization solver such as CPLEX or Gurobi using callbacks, while existing implementations of the perspective cut approach have required tailored branch-and-bound procedures [see, e.g., 28, Section $3.1$].

3.5 Merits of Ridge, Big-$M$ Regularization: Algorithmic Perspective

We now summarize the relative merits of applying either ridge or big-$M$* regularization from an algorithmic perspective:*

• •

As noted in our randomized rounding guarantees in Section 3.2, the two regularization methods provide comparable bound gaps when $2M\approx\gamma L$, while if $2M\ll\gamma L$, big-$M$* regularization provides smaller gaps, and if $2M\gg\gamma L$, ridge regularization provides smaller gaps.*

• •

For linear problems, ridge regularization limits dual degeneracy, while big-$M$* regularization does not. This benefit, however, has to be put in balance with the extra runtime and memory requirements needed for solving a quadratic, instead of linear, separation problem.*

In summary, the benefits of applying either big-$M$* or ridge regularization are largely even and depend on the specific instance to be solved. In the next section, we perform a sequence of numerical experiments on the problems studied in Section 2.1, to provide empirical guidance on which regularization approach works best when.*

4 Numerical Experiments

In this section, we evaluate our single-tree cutting-plane algorithm, implemented in Julia 1.0 using CPLEX $12.8.0$ and the Julia package JuMP.jl version $0.18.4$ **[21]**. We compare our method against solving the natural big-$M$* or MISOCP formulations directly, using CPLEX $12.8.0$. All experiments were performed on one Intel Xeon E*$5-2690$* v4 *$2.6$GHz CPU core and using $32$ GB RAM.

4.1 Overall Empirical Performance Versus State-of-the-Art

In this section, we compare our approach to state-of-the-art methods, and demonstrate that our approach outperforms the state-of-the-art for several relevant problems.

4.1.1 Network Design

We begin by evaluating the performance of our approach for the multi-commodity network design problem (5). We adapt the methodology of **[36]** and generate instances where each node $i\in[m]$ is the unique source of exactly one commodity ($k=m$). For each commodity $j\in[m]$, we generate demands according to $b^{j}_{j^{\prime}}=\lfloor\mathcal{U}(5,25)\rceil$ for $j^{\prime}\neq j$ and $b^{j}_{j}=-\sum_{j^{\prime}\neq j}b^{j}_{j^{\prime}}$, where $\lfloor x\rceil$ is the closest integer to $x$ and $\mathcal{U}(a,b)$ is a uniform random variable on $[a,b]$. We generate edge construction costs, $c_{e}$, uniformly on $\mathcal{U}(1,4)$, and marginal flow circulation costs proportionally to each edge length777Nodes are uniformly distributed over the unit square $[0,1]^{2}$. We fix the cost to be ten times the Euclidean distance.. The discrete set $\mathcal{Z}$ contains constraints of the form $\bm{z}_{0}\leq\bm{z}$, where $\bm{z}_{0}$ is a binary vector which encodes existing edges. We generate graphs which contain a spanning tree plus $pm$ additional randomly picked edges, with $p\in[4]$, so that the initial network is connected with $O(m)$ edges. We also impose a cardinality constraint $\bm{e}^{\top}\bm{z}\leq(1+5\%)\bm{z}_{0}^{\top}\bm{e}$, which ensures that the network size increases by no more than $5\%$. For each edge, we impose a capacity $u_{e}\sim\lfloor\mathcal{U}(0.2,1)B/A\rceil$, where $B=-\sum_{j=1}^{m}b_{j}^{j}$ is the total demand and $A=(1+p)m$. We penalize the constraint $\bm{x}\leq\bm{u}$ with a penalty parameter $\lambda=1,000$888We do so to allow for a fair comparison between big-$M$ and ridge regularization. By penalizing the capacity constraint, we remove a natural big-$M$ regularization term and no regularization can be considered as more natural than the other.. For big-$M$* regularization, we set $M=\sum_{j}|b_{j}^{j}|$, and take $\gamma=\tfrac{2}{m(m-1)}$ for ridge regularization.*

We apply our approach to large networks with $100$s nodes, i.e., $10,000$s edges, which is ten times larger than the state-of-the-art **[38, 36]**, and compare the quality of the incumbent solutions after an hour, since no approach could terminate up to a satisfiable optimality gap within this time limit. Note that we define the quality of a solution as its cost in absence of regularization, although we might have augmented the original formulation with a regularization term to compute the solution. As a result, we can compare the performance big-$M$ and ridge regularization directly, despite the fact that the optimization problems they solve are actually different. On the other hand, performance metrics that depend on the function being minimized, such as the optimality gap, would not permit such a comparison. In $100$ instances, our cutting plane algorithm with big-$M$* regularization provides a better solution $94\%$ of the time, by $9.9\%$ on average, and by up to $40\%$ for the largest networks. For ridge regularization, the cutting plane algorithm scales to higher dimensions than plain mixed-integer SOCP, returns solutions systematically better than those found by CPLEX (in terms of unregularized cost), by $11\%$ on average. Also, ridge regularization usually outperforms big-$M$ regularization, as reported in Table 3. Given how numerically challenging these optimization problems are, the optimality gaps returned by all methods are often uninformative ($>100\%$) - see Section C Table 10. Still, we observe that, with big-$M$ regularization, CPLEX systematically returns tighter optimality gaps that the cutting-plane approach, while with ridge regularization, the gaps obtained by the cutting-plane algorithm are tighter $86\%$ of the times. All in all, even artificially added, ridge regularization improves the tractability of outer approximation.*

4.1.2 Binary Quadratic Optimization

We study some of the binary quadratic optimization problems collated in the BQP library by **[53]**. Specifically, the bqp-$\{50,100,250,500,1000\}$* instances generated by [3], which have a cost matrix density of $0.1$, and the be-$100$ and be-$120.8$ instances generated by [15], which respectively have cost matrix densities of $1.0$ and $0.8$. Note that these instances were generated as maximization problems, and therefore we consider a higher objective value to be better. We warm-start the cutting-plane approach with the best solution found after $10,000$ iterations of Goemans-Williamson rounding [see 35]. We also consider imposing triangle inequalities [19] via lazy callbacks, for they substantially tighten the continuous relaxations.*

Within an hour, only the bqp-$50$* and bqp-$100$ instances could be solved by any approach considered here, in which case cutting-planes with big-$M$ regularization is faster than CPLEX (see Table 4). For instances which cannot be solved to optimality, although CPLEX has an edge in producing tighter optimality gaps for denser cost matrices, as depicted in Table 4, the cutting-plane method provides tighter optimality gaps for sparser cost matrices, and provides higher-quality solutions than CPLEX for all instances, especially as $n$ increases (see Table 5).*

We remark that the cutting plane approach has low peak memory usage compared with the other methods: For the bqp-$1000$* instances, cutting-planes without triangle inequalities was the only method which respected the *$32$GB memory budget. This is another benefit of decomposing Problem (1) into master and sub-problems.

4.1.3 Sparse Empirical Risk Minimization

For sparse empirical risk minimization, our method with ridge regularization scales to regression problems with up $p=100,000$s features and classification problems with $p=10,000$s of features **[12]**. This constitutes a three-order-of-magnitude improvement over previous attempts using big-$M$ regularization **[10]**. We also select features more accurately, as shown in Figure 1, which compares the accuracy of the features selected by the outer-approximation algorithm (in green) with those obtained from the Boolean relaxation (in blue) and other methods.

4.1.4 Sparse Principal Component Analysis

We applied our approach to sparse principal component analysis problems in **[9]**, and by (a) introducing either big-$M$* or ridge regularization and (b) introducing additional valid inequalities into the master problem, which we derived from the Gershgorin Circle Theorem [see 9, Section 2.3, for details] successfully solved problems where $p=100$s to certifiable optimality, and problems where $p=1000$s to certifiable near optimality, as reported in Table 6; we refer to **[9] for descriptions of the datasets studied and more extensive numerical experiments. This constitutes an order-of-magnitude improvement over existing certifiably near-optimal approaches, which rely on semidefinite techniques and therefore cannot scale to *$p=1000$s.

4.1.5 Sparse Portfolio Selection

We applied our approach to sparse portfolio selection problems in **[8]**. By introducing a ridge regularization term, we successfully solved instances to optimality at a scale of one order of magnitude larger than previous attempts as summarized in Table 7. Specifically, we optimized over the securities in the Wilshire $5000$, which contains around $3,200$ securities, an improvement upon existing techniques, which cannot currently scale beyond the securities in the $S\&P$ 500. Moreover, at smaller scales which existing techniques have been benchmarked on—including the set of synthetic instances generated by **[28]** with $200-400$ securities—our approach is as fast as and often faster than existing state-of-the-art approaches including **[55, 27]** among others [see 8, Section 5.2, for details].

4.2 Evaluation of Different Ingredients in Our Numerical Recipe

We now consider the capacitated facility problem (10) on 112 real-world instances available from the OR-Library **[3, 39]**, with the natural big-$M$* and the ridge regularization with $\gamma=1$. In both cases, the algorithms return the true optimal solution. Compared to CPLEX with big-$M$ regularization, our cutting plane algorithm with big-$M$ regularization is faster in $12.7\%$ of instances (by $53.6\%$ on average), and in $23.85\%$ of instances (by $54.5\%$ on average) when using a ridge penalty. This observation suggests that ridge regularization is better suited for outer-approximation, most likely because, as discussed in Section 3.1, a strongly convex ridge regularizer breaks the degeneracy of the separation problems. Note that our approach could benefit from multi-threading and restarting.*

We take advantage of these instances to breakdown the independent contribution of each ingredient in our numerical recipe in Table 8. Although each ingredient contributes independently, jointly improving the lower and upper bounds provides the greatest improvement.

4.3 Big-$M$ Versus Ridge Regularization

In this section, our primary interest is in ascertaining conditions under which it is advantageous to solve a problem using big-$M$* or ridge regularization, and argue that ridge regularization is preferable over big-$M$ regularization as soon as the objective is sufficiently strongly convex.*

To illustrate this point, we consider large instances of the thermal unit commitment problem originally generated by **[29]**, and multiply the quadratic coefficient $a_{i}$ for each generator $i$ by a constant factor $\alpha\in\{0.1,1,2,5,10\}$. Table 9 depicts the average runtime for CPLEX to solve both formulations to certifiable optimality, or provides the average bound-gap whenever CPLEX exceeds a time limit of $1$ hour. Observe that when $\alpha\leq 1$, the big-$M$* regularization is faster, but, when $\alpha>1$ the MISOCP approach converges fast while the big-$M$ approach does not converge within an hour. Consequently, ridge regularization performs more favorably whenever the quadratic term is sufficiently strong.*

We also compare big-$M$* and ridge regularization for the sparse portfolio selection problem (7). Figure 2 depicts the relationship between the optimal allocation of funds $\bm{x}^{\star}$ and the regularization parameter $M$ (left) and $\gamma$ (right), and Figure 3 depicts the magnitude of the gap between the optimal objective and the Boolean relaxation’s objective, normalized by the unregularized objective. The two investment profiles are comparable, selecting the same stocks. Yet, we observe two main differences: First, setting $M<\frac{1}{k}$ renders the entire problem infeasible, while the problem remains feasible for any $\gamma>0$. This is a serious practical concern in cases where a lower bound on the value of $M$ is not known a priori. Second, the profile for ridge regularization seems smoother than its equivalent with big-$M$.*

4.4 Relative Merits of Big-$M$, Ridge Regularization: An Experimental Perspective

We now conclude our comparison of big-$M$* and ridge regularization, as initiated in Sections 2.4 and 3.5, by indicating the benefits of big-$M$ and ridge regularization, from an experimental perspective:*

• •

As observed in Section 4.3, big-$M$* and ridge regularization play fundamentally the same role in reformulating logical constraints. This observation echoes our theoretical analysis in Section 2.*

• •

As observed in the unit commitment and sparse portfolio selection problems studied in Section 4.3, ridge regularization should be the method of choice whenever the objective function contains a naturally occurring strongly convex term, which is sufficiently large.

• •

As observed for network design and capacitated facility location problems in sections 4.1.1-4.2, ridge regularization is usually more amenable to outer-approximation than big-$M$* regularization, because it eliminates most degeneracy issues associated with outer-approximating MINLOs.*

• •

The efficiency of outer-approximation schemes relies on the speed at which separation problems are solved. In this regard, special problem-structure or cardinality constraints on the discrete variable $\bm{z}$ drastically help. This has been the case in network design, sparse empirical risk minimization and sparse portfolio selection problems in Section 4.1.1.

5 Conclusion

In this paper, we proposed a new interpretation of the big-$M$* method, as a regularization term rather than a modeling trick. By expanding this regularization interpretation to include ridge regularization, we considered a wide family of relevant problems from the Operations Research literature and derived equivalent reformulations as mixed-integer saddle-point problems, which naturally give rise to theoretical analysis and computational algorithms. Our framework provides provably near-optimal solutions in polynomial time via solving Boolean relaxations and performing randomized rounding999By “polynomial time”, we mean with respect to the dimensionality of the relaxation, assuming the relaxation is solved to a fixed and finite precision and can be described using the symmetric cones described by [44], as occurs for all examples discussed in this paper. as well as certifiably optimal solutions through an efficient branch-and-bound procedure, and indeed frequently outperforms the state-of-the-art in numerical experiments.*

We believe our framework, which decomposes the problem into a discrete master problem and continuous subproblems, could be extended more generally to mixed-integer semidefinite optimization, as developed in **[11, 9]**.

Acknowledgments

*We thank the associate editor and the two anonymous referees for their valuable comments which improved the paper.

---

Appendix A Omitted Proofs

A.1 Proof of Theorem 3.1: Quality of the Random Rounding Strategy

Proof A.1**.**

We only detail the proof for the big-$M$ regularization case, as the ridge regularization case follows mutatis mutandis. From Proposition 2.8,

[TABLE]

The polyhedron $\{\bm{\alpha}:\bm{\alpha}\geq\bm{0},\|\bm{\alpha}\|_{1}\leq 1\}$ admits $|\mathcal{R}|+1$ extreme points. However, if

[TABLE]

for some $t>0$, then the maximum can only occur at some $\alpha>\mathbf{0}$ so that we can restrict our attention to the $|\mathcal{R}|$ positive extreme points. Applying tail bounds on the maximum of sub-Gaussian random variables over a polytope [see 51, Theorem 1.16], since $\|\bm{\alpha}\|_{2}\leq\|\bm{\alpha}\|_{1}\leq 1$, we have for any $t>0$,

[TABLE]

so that

[TABLE]

A.2 Proof of Theorem 3.2: Relationship With Perspective Cuts

Proof A.2**.**

Let us fix $\bm{z}\in\mathcal{Z}$. Then, we have that:

[TABLE]

Finally, observing that

[TABLE]

*concludes the proof. *

Appendix B Bounding the Lipschitz Constant

In our results, we relied on the observation that there exists some constant $L>0$ such that, for any $\bm{z}\in\mathcal{Z}$, $\|\bm{\alpha}^{\star}(\bm{z})\|\leq L$. Such an $L$ always exists, since $\mathcal{Z}$ is a finite set. However, as our randomized rounding results depend on $L$, explicit bounds on $L$ are desirable.

We remark that while our interest is in the Lipschitz constant with respect to “$\bm{\alpha}$” in a generic setting, we have used different notation for some of the problems which fit in our framework, in order to remain consistent with the literature. In this sense, we are also interested in obtaining a Lipschitz constant with respect to $\bm{w}$ for the portfolio selection problem (7), among others.

In this appendix, we bound the magnitude of $L$ in a less conservative manner. Our first result provides a bound on $L$ which holds whenever the function $h(\bm{\alpha})$ in Equation (16) is strongly concave in $\bm{\alpha}$, which occurs for the sparse ERM problem (6) with ordinary least-squares loss, the unit commitment problem (9), the portfolio selection (7), and network design problems whenever $\bm{\Sigma}$ (resp. $\bm{Q}$) is full-rank:

Lemma B.1**.**

Let $h(\cdot)$ be a strongly concave function with parameter $\mu>0$ [see 17, Chapter 9.1.2 for a general theory of strong convexity], and suppose that $\bm{0}\in\text{dom}(g)$ and $\bm{\alpha}^{\star}:=\arg\max_{\bm{\alpha}}\,h(\bm{\alpha})$. Then, for any choice of $\bm{z}$, we have

[TABLE]

*i.e., $\|\bm{\alpha}^{\star}(\bm{z})\|_{\infty}\leq L$, where $L:=2\sqrt{2\frac{h(\bm{\alpha}^{\star})-h(\bm{0})}{\mu}}$. *

Proof B.2**.**

By the definition of strong concavity, for any $\bm{\alpha}$ we have

[TABLE]

where $\nabla h(\bm{\alpha}^{\star})^{\top}(\bm{\alpha}-\bm{\alpha}^{\star})\leq{0}$ by the first-order necessary conditions for optimality, leading to

[TABLE]

In particular for $\bm{\alpha}=\bm{0}$, we have

[TABLE]

and for $\bm{\alpha}=\bm{\alpha}^{\star}(\bm{z})$,

[TABLE]

since

[TABLE]

*The result then follows by the triangle inequality. *

An important special case of the above result arises for the sparse ERM problem, as we demonstrate in the following corollary to Lemma B.1:

Corollary B.3**.**

For the sparse ERM problem (6) with an ordinary least squares loss function and a cardinality constraint $\bm{e}^{\top}\bm{z}\leq k$, a valid bound on the Lipschitz constant is given by

[TABLE]

*where $\bm{X}_{i,[k]}$ is the sum of the $k$ largest entries in the column $\bm{X}_{i,[k]}$. *

Proof B.4**.**

Applying Lemma B.1 yields the bound

[TABLE]

after observing that we can parameterize this problem in terms of $\bm{\alpha}$, and for this problem:

Setting $\bm{\alpha}=0$ yields $h(\bm{\alpha})=0$. 2. 2.

$0\leq h(\bm{\alpha}^{\star})\leq\bm{y}^{\top}\bm{\alpha}^{\star}-\frac{1}{2}\bm{\alpha}^{\star\top}\bm{\alpha}^{\star}\leq\frac{1}{2}\bm{y}^{\top}\bm{y}$. 3. 3.

$h(\cdot)$* is strongly concave in $\bm{\alpha}$, with concavity constant $\mu\geq 1$.*

*The result follows by applying the definition of the operator norm, and pessimizing over $\bm{z}$. *

Appendix C Supplementary material for the numerical experiments

In this section, we report additional performance metrics for the network design experiments presented in Section 4.1.1. There, we reported the quality of the solution returned by all methods within one hour, and compared two regularization strategies (big-$M$* vs. ridge) and two algorithms (CPLEX vs. Cuts). Indeed, given the sizes of problems considered, the network design instances are computationally very challenging to solve. At such scales, finding a good feasible solution is already a very difficult task. In practice, this translates into optimality gaps that are often irrelevant (i.e., higher than $100\%$) in most of the instances. Table 10 reports the optimality gaps returned by each method after one hour, on instances where at least one of the four gaps was less than $100\%$.*