Théo Guyardtg
Cédric Herzetch
Clément Elvirace
Ayşe-Nur Arslanaa

tgInria and Insa, Univ Rennes, CNRS, IRMAR - UMR 6625, Rennes, France
chEnsai, Univ Rennes, CNRS, CREST - UMR 9194, Rennes, France
ceCentraleSupélec, Univ Rennes, CNRS, IETR - UMR 6164, Rennes, France
aaInria, Univ Bordeaux, CNRS, IMB - UMR 5251, Bordeaux, France

Théo<EMAIL_ADDRESS>
L0-regularization, Branch-and-Bound, Pruning test

We consider the resolution of learning problems involving $\ell_0$-regularization via bnb algorithms. These methods explore regions of the feasible space of the problem and check whether they do not contain solutions through “pruning tests”.
In standard implementations, evaluating a pruning test requires to solve a convex optimization problem, which may result in computational bottlenecks.
In this paper, we present an alternative to implement pruning tests for some generic family of $\ell_0$-regularized problems. Our proposed procedure allows the simultaneous assessment of several regions and can be embedded in standard bnb implementations with a negligible computational overhead. We show through numerical simulations that our pruning strategy can improve the solving time of bnb procedures by several orders of magnitude for typical problems encountered in machine-learning applications.

§ INTRODUCTION

This paper focuses on optimization problems of the form:

\begin{equation}
\stepcounter{equation}
\tag{$\mathscr{P}$}
\label{prob:prob}
\opt{\pobj} = \textstyle\inf_{\pv \in \kR^{\pdim}} \lossfunc(\dic\pv)
+ \regfunc(\pv)
\end{equation}

where $\lossfunc: \kR^{\ddim} \mapsto \kR\cup\{+\infty\}$ is a loss function, $\dic\in \kR^{\ddim\times\pdim}$ is a given matrix and $\regfunc: \kR^{\pdim} \mapsto \kR\cup\{+\infty\}$ is a regularization function expressed as

\begin{equation}
\label{eq:regfunc}
\textstyle
\regfunc(\pv) = \reg\norm{\pv}{0} + \sum_{\idxentry=1}^\pdim\pertfunc(\pvi{\idxentry})
\end{equation}

for some $\reg > 0$ and $\pertfunc: \kR \mapsto \kR \cup \{+\infty\}$.
On the one hand, the so-called “$\ell_0$-norm” is defined as

\begin{equation*}
\norm{\pv}{0} = \mathrm{card}(\kset{\idxentry \in \intervint{1}{\pdim}}{\pvi{\idxentry}\neq 0})
\end{equation*}

and promotes sparsity in the optimizers of (<ref>) by counting the number of non-zero elements in its argument.
On the other hand, the term $\pertfunc(\cdot)$ allows to enforce additional application-specific properties, see , [Bruer et al., 2015, Bertsimas et al., 2021].
Solving problem (<ref>) is of interest in many fields including machine learning, high-dimensional statistics or signal processing.
This problem is for instance linked to feature selection [Bertsimas et al., 2016], compressive sensing [Candes et al., 2007], principal component analysis [Bertsimas & Cory-Wright, 2022], sparse SVM [Tan et al., 2010] or neural network pruning [Carreira-Perpinán & Idelbayev, 2018], among others.
The reader can refer to [Tillmann et al., 2021, Bertsimas et al., 2021] for an extensive review of related applications.

Since (<ref>) is NP-hard [Nguyen et al., 2019], the main trends of work in the last decades have focused on addressing relaxed instances of this problem or inferring its solutions through heuristic procedures [Tropp & Wright, 2010].
However, it has recently been emphasized that the solutions of the original problem (<ref>) may enjoy much better statistical properties than those obtained by these sub-optimal strategies [Bertsimas & Van Parys, 2020, Zhong et al., 2022].
Consequently, there has recently been a revived interest in solving (<ref>) exactly and several studies have emphasized that discrete-optimization tools can sometimes provide tractable solutions, see  [Bertsimas et al., 2016].
In this vein, state-of-the-art procedures are mostly based on bnb algorithms whose process can be specialized to exploit the structure of (<ref>) and achieve competitive running times [Ben Mhenni et al., 2022, Hazimeh et al., 2022].

In a nutshell, bnb algorithms solve an optimization problem by successively:
i) dividing the feasible space into regions and
ii) trying to detect regions that cannot contain a minimizer.
This second step is commonly referred to as “pruning test” and is based on the construction of some lower bounds on the value that the objective function can take.
A standard approach to constructing these lower bounds is based on the minimization of a convex lower-approximation of the objective function of (<ref>), called “relaxation”.
This operation usually dominates the complexity of bnb algorithms and can lead to tractability issues for some problem instances.
In this work, we propose a strategy to alleviate this computational bottleneck.

§.§ Contributions

We present a novel methodology to implement pruning tests in bnb algorithms that does not leverage the solution of a convex optimization problem and show that it can be embedded at virtually no cost in standard bnb implementations.
Numerical simulations reveal that our strategy can improve the solving time by several orders of magnitude.
Our contribution follows some recent lines of work [Atamturk & Gómez, 2020, Guyard et al., 2022] and exploits Fenchel-Rockafellar duality [Rockafellar, 1967].
In contrast to these prior contributions which focused on specific[
Prior contributions focused on instances defined by specific choices of the functions $\lossfunc(\cdot)$ and $\pertfunc(\cdot)$.
These works referred to their methodology as “screening”.
In this paper, we rather use the terminology “pruning” which is standard in the bnb literature.
] instances of problem (<ref>), we introduce a general framework encompassing a large family of problems typically encountered in machine learning.
Specifically, our framework applies under the following set of hypotheses:

* The function $\lossfunc(\cdot)$ is proper, closed and convex.

* The function $\pertfunc(\cdot)$ is proper, closed and convex.

* $\pvi{}=0$ is an accumulation point[
\(\pvi{}\) is said to be an accumulation point of a set \(\proofset\subseteq\kR\) if for all neighborhoods \(\mathcal{N}\) of \(\pvi{}\), the set \(\mathcal{N}\cap\proofset\setminus\{\pvi{}\}\) is nonempty.
] of $\dom(\pertfunc)$.

* \(0\in\dom(\pertfunc)\) and \(\pertfunc(0) = 0\).

Hypotheses <ref>-<ref> are verified by many functions encountered in standard machine-learning problems.
For example, this includes least-squares, logistic or hinge losses [Wang et al., 2020] and terms $\pertfunc(\cdot)$ constructed as mixed-norms [Dedieu et al., 2021] or as the logarithm of Bayesian priors [Polson & Sokolov, 2019].

Finally, the complexity analysis of our method (see <Ref>) is discussed in view of the following assumption which holds for a wide range of problem instances encountered in practice, see , Section 4.4.16 in [Beck, 2017]:

* The evaluation complexity of the convex conjugate of \(\lossfunc(\cdot)\) scales as \(\bigO(\ddim)\).

§.§ Outline

The rest of the paper is organized as follows.
In <Ref>, we introduce the main ingredients of standard bnb algorithms.
In <Ref>, we then present our new pruning strategy and discuss its impact on the bnb algorithm.
Finally, our method is assessed numerically in <Ref>.
To ease our exposition, all the proofs are deferred to <Ref>.

§.§ Notational Conventions

Classical letters (, $x$), boldface lowercase letters (, $\pv$) and boldface uppercase letters (, $\dic$) represent scalars, vectors and matrices, respectively.
$\0$ and $\1$ denote the all-zero and all-one vectors whose dimension is usually clear from the context.
Vectorial operations involving equalities or inequalities have to be understood coordinate-wise.
We note $\pvi{\idxentry}$  the $\idxentry$-th entry of a vector $\pv$ and $\pv_{\setidx}$ its restriction to the entries indexed by some set of indices $\setidx$.
Similarly, we note $\atom{\idxentry}$ the $\idxentry$-th column of a matrix $\dic$ and $\dic_{\setidx}$ its restriction to the columns indexed by $\setidx$.
$\intervint{a}{b}$ corresponds to the set of integers ranging from $a$ to $b$.
The notations $\card{\cdot}$ and $\cdot\setminus\cdot$ are used to denote respectively the cardinality of a set and the difference between two sets.
$\icvx(\cdot)$ stands for the convex indicator function defined as $\icvx(\cdot) = 0$ if the condition in the parentheses is fulfilled and $\icvx(\cdot) = +\infty$ otherwise.
We let $\pospart{\pvi{}} = \max(\pvi{},0)$.
Given some proper function $\genericfunc(\cdot)$, we note $\dom(\genericfunc)$ its domain, $\conj{\genericfunc}(\cdot)$ its convex conjugate, $\biconj{\genericfunc}(\cdot)$ its convex  biconjugate
and $\subdiff\genericfunc(\cdot)$ its subdifferential.
We refer to [Beck, 2017] for a precise definition of these notions.
Finally, we employ the notational convention $\genericfunc(\pv) = \sum_{\idxentry}\separable{\genericfunc}{\idxentry}(\pvi{\idxentry})$ when
$\genericfunc(\cdot)$ is separable.

§ BRANCH-AND-BOUND ALGORITHMS

In this section, we outline the main ingredients of bnb procedures.
We focus on the elements necessary to present our contribution and refer the reader to Chap. 5 in [Locatelli & Schoen, 2013] for a thorough description.

Illustration of the bnb decision-tree exploration. We note that a relaxation has to be solved at each node of the tree to evaluate the lower bound (<ref>) involved in pruning test (<ref>).
Here, the pruning test is passed for nodes $\nodeSymb_2$ and $\nodeSymb_4$.

§.§ Constructing and Pruning Regions

bnb algorithms partition the feasible space of an optimization problem into regions and try to detect those that do not contain any optimizer.
In the specific context of (<ref>), standard bnb implementations consider regions of the form

\begin{equation} \label{eq:region}
\node{\pset} = \kset{\pv \in \kR^{\pdim}}{\pv_{{\setzero}} = \0, \pv_{{\setone}} \neq \0, \pv_{{\setnone}} \in \kR^{\card{{\setnone}}}}
\end{equation}

where $\nodeSymb=(\setzero, \setone, \setnone)$ is a partition of $\intervint{1}{\pdim}$, see [Ben Mhenni et al., 2020, Hazimeh et al., 2022].

\begin{equation}
\stepcounter{equation}
\tag{$\node{\ppb}$}
\label{prob:node-prob}
\textstyle
\node{\pobj}
= \inf_{\pv \in \kR^{\pdim}} \lossfunc(\dic\pv) + \regfunc(\pv) + \icvx(\pv \in \node{\pset})
\end{equation}

one can therefore deduce that no minimizer of (<ref>) is contained in region $\node{\pset}$ when the inequality

\begin{equation}
\label{eq:ideal-pruning-test}
\node{\pobj} > \opt{\pobj}
\end{equation}

is verified.
Unfortunately, condition (<ref>) is of little interest in practice since evaluating $\node{\pobj}$ and $\opt{\pobj}$ is an NP-hard task.
A workaround to this issue consists in relaxing (<ref>) as

\begin{equation}
\label{eq:pruning-test}
\node{\LB{\pobj}} > \UB{\pobj}
\end{equation}

where $\node{\LB{\pobj}}$ and $\UB{\pobj}$ are some tractable lower and upper bounds on $\node{\pobj}$ and $\opt{\pobj}$, respectively.
Inequality (<ref>) is often referred to as “pruning test” since if it is verified, $\node{\pset}$ does provably not contain any minimizer of (<ref>) and can therefore be safely pruned from the optimization problem.

From the point of view of pruning efficiency, solving (<ref>) at the root node \(\nodeSymb_0\) is unnecessary as the pruning test (<ref>) is never passed.
It is nevertheless common practice to solve it as a wide range of branching rules used in practice leverage its solutions, (see e.g., <Ref>).

§.§ Standard Bounding Strategy

The standard strategy to construct the bounds involved in (<ref>) is as follows.
First, an upper bound $\UB{\pobj}$ on $\opt{\pobj}$ can be computed by evaluating the objective function of (<ref>) at any feasible point.[Many methods allow to construct relevant candidates at a reasonable cost, see , [Wolsey, 1980].]
Second, the computation of $\node{\LB{\pobj}}$ is generally done by minimizing some convex lower bound on the objective function of (<ref>).
More specifically, a standard choice consists in replacing the term

\begin{equation}
\label{eq:node-regfunc}
\node{\regfunc}(\pv) = \regfunc(\pv) + \icvx(\pv \in \node{\pset})
\end{equation}

in problem (<ref>) by its convex biconjugate denoted $\biconj{(\node{\regfunc})}(\cdot)$ in this paper, see Item (i) of Proposition 13.16 in [Bauschke & Combettes, 2017].
Hence, a valid choice for the lower bound in pruning test (<ref>) reads

\begin{equation}
\label{eq:std-lb}
\node{\LB{\pobj}} = \node{\robj}
\end{equation}

\begin{equation} \label{prob:relax-node}
\stepcounter{equation}
\tag{$\node{\rpb}$}
\node{\robj} = \textstyle\inf_{\pv \in \kR^{\pdim}} \lossfunc(\dic\pv) +
\biconj{(\node{\regfunc})}(\pv).
\end{equation}

Problem (<ref>) is called a “relaxation” of (<ref>) and is usually addressed by first-order convex optimization methods [Beck, 2017].
The complexity of these algorithms typically scales as $\bigO(\ddim\pdim \kappa)$ where $\kappa$ denotes the number of iterations performed by the numerical procedure.

§.§ Feasible Space Exploration

In bnb procedures, the partitioning of the feasible set into regions can be identified with the expansion of a decision tree where each node corresponds to some region \(\node{\pset}\) defined as in (<ref>).
As illustrated in <Ref>, the exploration starts at the root node  \(\nodeSymb_0 = (\emptyset,\emptyset,\intervint{1}{\pdim})\) which corresponds to $\pset^{\nodeSymb_0}=\kR^\pdim$.
For each leaf node $\nodeSymb=(\setzero,\setone,\setnone)$ examined by the bnb procedure, problem (<ref>) is solved and pruning test (<ref>) is evaluated using lower bound (<ref>).
If the test is passed, the corresponding region $\pset^{\nodeSymb}$ is pruned from the problem and the exploration of the tree is stopped below this node.
If the pruning test is not passed, $\pset^{\nodeSymb}$ is partitioned into two new regions
as follows.
An index \(\idxentry\in\setnone\) is selected and the following two child nodes of $\nodeSymb$ are created:

\begin{align}
\label{eq:direct-zero}
\succnode{0}{\idxentry}&\triangleq ({\setzero}\cup\{\idxentry\},{\setone},{\setnone}\backslash\{\idxentry\})\\
\label{eq:direct-one}
\succnode{1}{\idxentry}&\triangleq({\setzero},{\setone}\cup\{\idxentry\},{\setnone}\backslash\{\idxentry\})
\end{align}

This process is repeated until all the leaf nodes of the tree are such that $\setnone=\emptyset$.
In the latter case, (<ref>) reduces to a convex optimization problem which can be commonly solved to machine precision.

§.§ Pruning Effectiveness and Method Complexity

The computational burden associated to some node $\nodeSymb$ is generally dominated by the resolution of (<ref>).
Since such a relaxation has to be solved at each node of the decision tree, the overall complexity of bnb procedures typically roughly scales linearly with the total number of nodes explored by the algorithm.
On the one hand, we note that passing pruning test (<ref>) often requires to compute tight lower bounds on \(\node{\pobj}\), in the sense that \(\node{\pobj}- \node{\robj}\) should be small.
On the other hand, achieving the prescribed tightness usually imposes exploring deep nodes in the decision tree.[This claim is discussed more precisely in <Ref>.]
As a consequence, the number of nodes explored by bnb procedures may become prohibitively large and thus lead to impractical solving times.
In the next section, we present a strategy to alleviate this computational bottleneck.

§ A NEW PRUNING STRATEGY

In this section, we present our proposed new strategy to implement pruning test (<ref>).
Our exposition is organized as follows.
In <Ref>, we first discuss another lower bound on \(\node{\pobj}\) and show that the latter can be evaluated at low cost for a set of nodes in <Ref>. Using this observation, we explain in <Ref> how the bnb decision tree can be expanded when several pruning tests exploiting this lower bound are passed simultaneously. Finally, in <Ref> we describe how to embed the proposed methodology in standard bnb implementations at virtually no cost.

§.§ Exploiting Duality to Construct Lower Bounds

We propose to construct a lower bound on $\node{\pobj}$ that does not require to solve (<ref>). To this end, we follow another line of work in the bnb literature [Sarin et al., 1988] and leverage the Fenchel-Rockafellar dual problem [Rockafellar, 1967] associated to (<ref>).
Under hypotheses <ref>-<ref>, the latter is given by[
Here, we use the fact that \(\conj{(\biconj{(\node{\regfunc})})}(\cdot) = \conj{(\node{\regfunc})}(\cdot)\) under <ref> by Proposition 13.16.

\begin{equation}
\stepcounter{equation}
\tag{$\node{\mathscr{D}}$}
\label{prob:dual-node}
\textstyle
\node{\dobj} = \sup_{\dv \in \kR^{\ddim}} \underbrace{-\conj{\lossfunc}(-\dv) - \conj{(\node{\regfunc})}(\transp{\dic}\dv)}_{\triangleq~\node{\dfunc}(\dv)}
\end{equation}

and verifies the inequality

\begin{equation}
\label{eq:bound-ordering}
\node{\dobj} \leq \node{\robj},
\end{equation}

which is tight under mild assumptions, see , Proposition 15.22 in [Bauschke & Combettes, 2017].
Hence, setting

\begin{equation}
\label{eq:dual-bound}
\node{\LB{\pobj}} = \node{\dfunc}(\dv)
\end{equation}

leads to a valid lower bound to implement (<ref>) for any \(\dv \in \kR^{\ddim}\).
Regarding inequality (<ref>), we note that lower bound (<ref>) may not be as tight as the standard one given in (<ref>).
Nevertheless, we show in the sequel that (<ref>) can be evaluated for several “successors” of node $\nodeSymb$ at virtually no cost.

Let us first precise the notion of “successors” of a node \(\nodeSymb\):

The node $\nodeSymb'=(\setzero',\setone',\setnone')$ is said to be a successor of $\nodeSymb=(\setzero,\setone,\setnone)$ if it verifies

\begin{equation}
\label{eq:successor}
\setzero \subseteq \setzero'
\quad\text{and}\quad
\setone \subseteq \setone'
\end{equation}

Moreover, $\nodeSymb'$ is said to be a direct successor of $\nodeSymb$ if it fulfills property (<ref>) and verifies

\begin{equation}
\label{eq:direct-successor}
(\setzero' \setminus \setzero) \cup (\setone' \setminus \setone) = \{\idxentry\}
\end{equation}

for some \(\idxentry \in \setnone\).

From a bnb tree perspective, direct successors correspond to the nodes \(\succnode{0}{\idxentry}\) and \(\succnode{1}{\idxentry}\) described in (<ref>)-(<ref>) for some \(\idxentry\in\setnone\).

Interestingly, the objective functions of the dual problems (<ref>) at some node $\nodeSymb$ and any of its successors share a similar mathematical structure.
To reveal this link, we first establish the following result in <Ref>:

Let $\nodeSymb=(\setzero,\setone,\setnone)$ be a node.
Under <ref>-<ref>, the function
$\conj{(\node{\regfunc})}(\cdot)$ is separable and defined coordinate-wise for all $\vvi{} \in \kR$ as

\begin{equation}
\label{eq:dual-relaxregfunc}
\conj{(\separable{\node{\regfunc}}{\idxentry})}(\vvi{}) =
\begin{cases}
0 &\text{if} \ \idxentry \in \setzero \\
\conj{\pertfunc}(\vvi{}) - \reg &\text{if} \ \idxentry \in \setone \\
\pospart{\conj{\pertfunc}(\vvi{}) - \reg} &\text{if} \ \idxentry \in \setnone.
\end{cases}
\end{equation}

As a consequence of <Ref>, we prove in <Ref> that the objective function of the dual problem (<ref>) at some node $\nodeSymb$ verifies a notable relation with that of its successors.
More precisely, letting

\begin{align}
\label{eq:pivot-zero}
\pivot{0}(\vvi{}) &\triangleq \pospart{\conj{\pertfunc}(\vvi{}) - \reg} \\
\label{eq:pivot-one}
\pivot{1}(\vvi{}) &\triangleq \pospart{\reg - \conj{\pertfunc}(\vvi{})}
\end{align}

for all \(\vvi{} \in \kR\), we obtain the following property:

Let $\nodeSymb' = (\setzero',\setone',\setnone')$ be a successor of $\nodeSymb = (\setzero,\setone,\setnone)$.
Under <ref>-<ref>, we have for all \(\dv \in \kR^{\ddim}\):

\begin{equation} \label{eq:dual-function}
\begin{split}
\dfunc^{\nodeSymb'}(\dv) = \node{\dfunc}(\dv) &+ \textstyle\sum_{\idxentry \in \setzero'\setminus\setzero} \pivot{0}(\transp{\atom{\idxentry}}\dv) \\
&+ \textstyle\sum_{\idxentry \in \setone'\setminus\setone} \pivot{1}(\transp{\atom{\idxentry}}\dv)
\end{split}
\end{equation}

From <Ref>, we note that the term $\node{\dfunc}(\dv)$ is common to the objective functions of all the dual problems associated to the successors of $\nodeSymb$.
As a consequence, lower bound (<ref>) can be computed for several successors $\nodeSymb'$ of $\nodeSymb$ through a single evaluation of this term.
In the next section, we leverage this observation to jointly evaluate lower bound (<ref>) at all the direct successors of \(\nodeSymb\) with a complexity scaling as \(\bigO(\ddim\pdim)\).

§.§ Evaluation of ( <ref>) for All the Direct Successors

In this section, we consider the scenario where one wants to evaluate lower bound (<ref>) for all the direct successors of some node $\nodeSymb=(\setzero,\setone,\setnone)$.
In this case, letting \(\nodeSymb'=(\setzero',\setone',\setnone')\) denote a direct successor of \(\nodeSymb\), we have that (<ref>) simplifies to

\begin{equation} \label{eq:dual-function-direct}
\hspace*{-3pt}
\dfunc^{\nodeSymb'}(\dv) =
\node{\dfunc}(\dv) +
\begin{cases}
\pivot{0}(\transp{\atom{\idxentry}}\dv) & \hspace*{-5pt}\mbox{if \(\setzero' \setminus \setzero=\{\idxentry\}\)}
\\
\pivot{1}(\transp{\atom{\idxentry}}\dv) & \hspace*{-5pt}\mbox{if \(\setone' \setminus \setone=\{\idxentry\}\)}.
\end{cases}
\end{equation}

We thus remark that the evaluation of (<ref>) for all the direct successors $\nodeSymb'$ of $\nodeSymb$ amounts to computing:

the inner products \(\{\transp{\atom{\idxentry}}{\dv}\}_{\idxentry=1}^{\pdim}\),

* the quantity $\node{\dfunc}(\dv)$,

* the terms $\pivot{0}(\transp{\atom{\idxentry}}\dv)$ and $\pivot{1}(\transp{\atom{\idxentry}}\dv)$ for all $\idxentry \in \setnone$.

The first task can obviously be done with a complexity $\bigO(\ddim\pdim)$.
Moreover, using the definition of $\node{\dfunc}(\dv)$ in (<ref>) and <Ref>,
we have that the last two tasks can be implemented with a complexity scaling as \(\bigO(\ddim+\pdim)\) and \(\bigO(\card{\setnone})\), respectively, under hypothesis <ref>.[This is achieved by re-using the computations of task $i)$.]
Overall, the computational burden induced by the evaluation of (<ref>) for all the direct successors of node $\nodeSymb$ scales as $\bigO(\ddim\pdim)$.
This is a substantial decrease with respect to the complexity required to perform the same task for standard lower bound (<ref>).
More specifically, applying $\kappa$ iterations of some first-order method to (<ref>) leads to a complexity of $\mathcal{O}(\ddim\pdim\kappa)$ (see <Ref>). Repeating this operation for all $2|\setnone|$ direct successors of $\nodeSymb$ then leads to a degradation of the computation burden by a factor $2|\setnone|\kappa$.

In fact, we will show in <Ref> that the complexity overhead needed to evaluate lower bound (<ref>) for all the direct successors of $\nodeSymb$ can be further reduced from $\bigO(\ddim\pdim)$ to $\bigO(\ddim+\pdim)$ within a standard bnb implementation.

§.§ Simultaneous Tests: Effect on the Tree Exploration

The evaluation at low cost of lower bound (<ref>) for all the direct successors $\nodeSymb'$ of \(\nodeSymb\) opens the way to the practical implementation of simultaneous pruning tests for several regions $\pset^{\nodeSymb'}$.
In this section, we investigate how the bnb decision tree can be expanded if pruning test (<ref>) is passed simultaneously for several direct successors of some node \(\nodeSymb=(\setzero,\setone,\setnone)\).

We first note that
\(
\kset{\succnode{0}{\idxentry}}{\idxentry\in\setnone}
\cup
\kset{\succnode{1}{\idxentry}}{\idxentry\in\setnone}
\), where \(\succnode{0}{\idxentry}\) and \(\succnode{1}{\idxentry}\) have been defined in (<ref>)-(<ref>),
corresponds to all the direct successors of \(\nodeSymb\).
Hence, the two sets

\begin{align}
\node{\setidxzero} &= \kset{\idxentry\in\setnone}{\dfunc^{\succnode{0}{\idxentry}}(\dv) > \UB{\pobj}}\label{def:I0}\\
\node{\setidxone} &= \kset{\idxentry\in\setnone}{\dfunc^{\succnode{1}{\idxentry}}(\dv) > \UB{\pobj}}\label{def:I1}
\end{align}

characterize the direct successors of $\nodeSymb$ satisfying pruning condition (<ref>) when implemented with the proposed lower bound (<ref>) at some dual point $\dv\in\kR^\ddim$.[
As discussed in <Ref>, the dual point $\dv$ is constructed from the iterates of the solving procedure addressing (<ref>).]
By definition of the sets $\node{\setidxzero}$ and $\node{\setidxone}$, the region

\begin{equation}
(\cup_{\idxentry\in\node{\setidxzero}} \pset^{\succnode{0}{\idxentry}})
\cup
(\cup_{\idxentry\in\node{\setidxone}} \pset^{\succnode{1}{\idxentry}})
\end{equation}

does not contain any minimizer of problem (<ref>).

In the rest of this paragraph, we show how the simultaneous success of several pruning tests can be translated in terms of deployment of the bnb decision tree.
To guide our reasoning, we examine two cases and then present a generic procedure to expand the bnb decision tree in <Ref>.

First, suppose that $\node{\setidxzero}\cap\node{\setidxone} \neq \emptyset$, that is there exists some $\idxentry\in{\setnone}$ such that both $\succnode{0}{\idxentry}$ and $\succnode{1}{\idxentry}$ pass the pruning test (<ref>) with lower bound (<ref>).
Since $\node{\pset}=\pset^{\succnode{1}{\idxentry}}\cup \pset^{\succnode{0}{\idxentry}}$, one concludes that no minimizers to problem (<ref>) can be found in region \(\node{\pset}\) which can thus be pruned from the decision tree.

Second, suppose that \(\node{\setidxzero}\cap\node{\setidxone} = \emptyset\) but \(\node{\setidxzero}\cup\node{\setidxone} \neq \emptyset\).
We first examine the case where only one direct successor of \(\nodeSymb\) passes the pruning test (<ref>) with lower bound (<ref>), i.e., $\node{\setidxzero}\cup\node{\setidxone}=\{\idxentry\}$, say $\node{\setidxzero}=\{\idxentry\}$ and \(\node{\setidxone}=\emptyset\) for instance.
Then, one concludes that region $\pset^{\succnode{0}{\idxentry}}$ does not contain any solution to problem (<ref>) and can therefore be pruned from the problem's feasible set without altering its solution.
From a decision-tree perspective, this information can be taken into account by: i) creating two new nodes $\succnode{0}{\idxentry}$ and $\succnode{1}{\idxentry}$ below $\nodeSymb$ and ii) immediately pruning the region \(\pset^{\succnode{0}{\idxentry}}\).
If \(\node{\setidxzero}\cup\node{\setidxone}\) contains more than one element, the next proposition suggests that this procedure can be applied recursively.

Under <ref>-<ref>, we have for all successors \(\nodeSymb'=(\setzero',\setone',\setnone')\) of \(\nodeSymb\) and indices \(\idxentry\in\setnone'\):

\begin{align}
\dfunc^{\succnode{0}{\idxentry}}(\dv) > \UB{\pobj}
&\;\implies\; \dfunc^{\nodeSymb'_{0,\idxentry}}(\dv) > \UB{\pobj}\\
\dfunc^{\succnode{1}{\idxentry}}(\dv) > \UB{\pobj}
&\;\implies\; \dfunc^{\nodeSymb'_{1,\idxentry}}(\dv) > \UB{\pobj}
\end{align}

A proof of this result can be found in <Ref>.
<Ref> states that if pruning test (<ref>) is passed using lower bound (<ref>)
for some direct successor of node $\nodeSymb$, the result of the test can be propagated to any successor $\nodeSymb'$ of $\nodeSymb$ compliant with the condition $\idxentry\in\setnone'$.
This observation leads to <Ref> which describes how the bnb decision tree can be expanded upon the knowledge of $\node{\setidxzero}$ and $\node{\setidxone}$.
<Ref> illustrates the output of <Ref> when $\nodeSymb=\nodeSymb_0$, $\setidxzero^{\nodeSymb_0}=\emptyset$ and $\setidxone^{\nodeSymb_0}=\{\idxentry_0,\idxentry_1\}$.
In comparison to <Ref>, we observe that the proposed pruning procedure does not require to solve the relaxations at the nodes $\nodeSymb_1$, $\nodeSymb_2$ and $\nodeSymb_4$, although ultimately leading to the same expanded tree.

node $\nodeSymb$, sets $\node{\setidxzero}$, $\node{\setidxone}$ defined in (<ref>)-(<ref>)

$\node{\setidxzero}\cap \node{\setidxone}\neq\emptyset$Prune \(\pset^{\nodeSymb}\) from the bnb tree

Set $\nodeSymb' \leftarrow \nodeSymb$

$\idxentry \in \node{\setidxzero} \cup \setidxone$
Create the two direct successors $\nodeSymb'_{0,\idxentry}$ and $\nodeSymb'_{1,\idxentry}$ to $\nodeSymb'$

Prune \(\pset^{\nodeSymb'_{0,\idxentry}}\) and set $\nodeSymb' \leftarrow \nodeSymb'_{1,\idxentry}$
Prune \(\pset^{\nodeSymb'_{1,\idxentry}}\) and set $\nodeSymb' \leftarrow \nodeSymb'_{0,\idxentry}$

Tree expansion based on simultaneous pruning

§.§ Implementation in Branch-and-Bound Methods

In this final section, we discuss how the proposed pruning strategy, leveraging lower bound (<ref>), can be efficiently integrated into standard bnb implementations.

We consider the following strategy.
Given some node $\nodeSymb$ processed by the bnb procedure, we test simultaneously all the direct successors of $\nodeSymb$ at each iteration of the solving process of (<ref>).
More specifically, letting $\hat{\pv} \in \kR^{\pdim}$ be the current iterate constructed by the first-order method addressing (<ref>), we evaluate the lower bound (<ref>) with some dual point verifying

\begin{equation} \label{eq:dual-point}
\dv \in -\subdiff\lossfunc(\dic\hat{\pv})
\end{equation}

and construct the sets $\node{\setidxzero}$ and $\node{\setidxone}$ according to (<ref>)-(<ref>).
We defer the discussion on our motivations for choosing \(\dv\) as in (<ref>) to the next paragraph.
We stop (prematurely) the resolution of (<ref>) as soon as $\node{\setidxzero}\cup\node{\setidxone}\neq\emptyset$ for some $\dv \in \kR^{\ddim}$ and expand the decision tree as described in <Ref>.
The bnb algorithm is then continued and this process is repeated.
If none of the pruning tests are passed during the resolution of (<ref>), the standard pruning test using lower bound (<ref>) is applied and the decision tree is expanded according to standard bnb operations as described in <Ref>.

We devote the rest of the section to motivate the choice of dual point \(\dv \in \kR^{\ddim}\) described in (<ref>).
From an effectiveness point of view, our rationale is to (try to) maximize the first term in (<ref>).
If strong duality holds between (<ref>) and (<ref>), by virtue of Theorem 19.1 from [Bauschke & Combettes, 2017], this can be achieved by choosing \(\dv \in - \partial \lossfunc(\dic\opt{\pv})\) where \(\opt{\pv}\) denotes any minimizer of (<ref>).
Since such a minimizer is not available, we use the current iterate (denoted \(\hat{\pv}\) in (<ref>)) of the numerical procedure solving (<ref>) as a surrogate.

From a complexity point of view, this proposed pruning methodology can be integrated within standard bnb implementations at virtually no cost.
Indeed, we notice that a dual point $\dv \in \kR^{\ddim}$ verifying (<ref>) and the corresponding vector $\transp{\dic}\dv \in \kR^{\pdim}$ are already computed during the iterations of many standard first-order methods tailored to solve (<ref>) such as proximal gradient, coordinate descent or ADMM [Beck, 2017].
According to our discussion in <Ref>, the complexity overhead required to compute lower bound (<ref>) associated to all the direct successors of $\nodeSymb$ then drops from $\bigO(\ddim\pdim)$ to $\mathcal{O}(\ddim+\pdim)$ since the inner products $\{\ktranspose{\atom{}}_i{\dv}\}_{i=1}^{\pdim}$ involved in task <ref> are already available.
We note that this additional computational burden is negligible as compared to the complexity of standard first-order methods which typically scales as $\mathcal{O}(\ddim\pdim)$ per iteration.

Impact of simultaneous pruning tests on the bnb tree exploration. Output of <Ref> when applied with $\nodeSymb=\nodeSymb_0$, $\setidxzero^{\nodeSymb_0}=\emptyset$ and $\setidxone^{\nodeSymb_0}=\{\idxentry_0,\idxentry_1\}$.

§ NUMERICAL EXPERIMENTS

In this final section, we assess numerically the proposed pruning strategy to accelerate bnb algorithms addressing problem (<ref>).
We do not focus on the statistical characterization of the solutions but refer to [Bertsimas et al., 2020, Hastie et al., 2020] for a thorough discussion on this topic.

The research presented in this paper is reproducible.
The associated code is open-sourced[<https://github.com/TheoGuyard/El0ps>] and all the datasets used in our simulations are publicly available.
Computations were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by INRIA and including CNRS, RENATER and several universities as well as other organizations.[<https://www.grid5000.fr>]
Experiments were run on a Debian 10 operating system, featuring one Intel Xeon E5-2660 v3 CPU clocked at 2.60 GHz with 16 GB of RAM.

Solver specifications
In our comparisons, we consider different methods solving problem (<ref>) exactly, that is returning the value of (at least) one minimizer to machine precision.
First, we use ,  and  which are off-the-shelf mip solvers [Anand et al., 2017].
Second, we consider the  solver [Hazimeh et al., 2022] which is dedicated to some specific[Namely with a quadratic function $\lossfunc(\cdot)$ and a term $\pertfunc(\cdot)$ corresponding to an $\ell_2$-norm and/or a bound constraint.] instances of problem (<ref>).
We compare these procedures to a standard bnb implementation enhanced with the simultaneous pruning tests described in this paper, noted .

For the sake of reproducibility, the mip formulations of the problem considered are specified in <Ref>.
<Ref> details our bnb implementation choices and <Ref> gives the expression of the function $\biconj{(\node{\regfunc})}(\cdot)$ involved in the relaxations (<ref>) for the instances of problem (<ref>) considered in this section.
Finally, the tuning procedure for $\reg$ and the hyperparameters involved in $\pertfunc(\cdot)$ is described in <Ref>.

§.§ Performance on Synthetic Data

In this section, we analyze the performance of different solvers on synthetic data.
We consider instances of problem (<ref>) defined by

\begin{align}
\label{eq:leastsquares}
\lossfunc(\cdot) &= \tfrac{1}{2}\norm{\obs-\cdot}{2}^2 \\
\label{eq:bigm}
\pertfunc(\cdot) &= \icvx(\abs{\cdot} \leq \bigM)
\end{align}

for some $\obs \in \kR^{\ddim}$ and $\bigM > 0$.
This choice is motivated by various applications, see , [Tillmann et al., 2021, Bertsimas et al., 2021, Bertsimas & Johnson, 2023].
Our results are averaged over 100 instances independently generated.

Instance generation

For each problem instance, we generate the rows of $\dic \in \kR^{\ddim\times\pdim}$ as independent realizations of a multivariate normal distribution with zero mean and covariance matrix $\corrmat \in \kR^{\pdim\times \pdim}$ where each entry $(i,j)$ is defined as $\corrmatel{ij}=\corrparam^{\abs{i-j}}$ for some $\corrparam \in [0,1)$.
Moreover, we set $\obs = \dic\groundtruth + \noise$ where $\groundtruth \in \kR^{\pdim}$ has $\sparsitylevel$ evenly-spaced non-zero entries of unit amplitude in absolute value and where $\noise \in \kR^{\ddim}$ is a zero-mean Gaussian noise with a variance tuned to obtain some signal-to-noise ratio $\snr = 10 \log_{10}(\norm{\dic\groundtruth}{2}^2/\norm{\noise}{2}^2)$.

Performance profiles

We generate each problem instance as described above with the parameters $\sparsitylevel=5$, $\ddim=500$, $\pdim=1000$, $\corrparam=0.9$ and $\snr=10$.
<Ref> represents the percentage of instances solved (to machine precision) by each method within a given time budget.

Performance profiles of different solvers.

We notice that regardless of the considered time budget, our method can solve a larger proportion of instances than its competitors.
In particular, all the problem instances are solved within a time budget for which no instances have been solved by the other methods.
More specifically, we observe that our methodology enables an acceleration of at least one order of magnitude with respect to the other procedures to solve all the problem instances.
In particular, we mention that the bnb implementation choices of  are similar to those of .
This suggests that the observed improvement in terms of computation time is essentially due to our simultaneous pruning strategy.

Sensibility study
To study more finely the gains permitted by the contribution proposed in this paper, we compare two different versions of a bnb algorithm.
The first one implements the standard pruning strategy with the lower bound (<ref>) whereas the second one implements both this standard strategy and the proposed simultaneous pruning methodology involving lower bound (<ref>).
We generate synthetic problem instances as described above by varying one parameter at a time to cover different working regimes.
In <Ref>, we represent the acceleration factor in terms of solving time obtained by our novel pruning strategy.
More precisely, it is defined as the ratio of the solving times obtained with  where the simultaneous pruning tests (see <Ref>) have been disabled and .

Acceleration factor when implementing the simultaneous pruning tests in addition to the standard pruning strategy during the bnb algorithm.

Time to construct solutions with a given sparsity level. The black dotted line represents the maximum time of one hour allowed.

We remark that the gains obtained by our methodology increase with the dimension parameter $\pdim$ and the correlation parameters $\corrparam$.
In contrast, the gain does not seem to be significantly impacted by the sparsity level $\sparsitylevel$ or the signal-to-noise ratio $\snr$.
This suggests that depending on the characteristics of the problem, the proposed pruning methodology can lead to different gains in terms of running time.
We note nonetheless that in all the tested scenarios, the solving time is improved by at least a factor $5$.

§.§ Performance on Real-World Datasets

In this section, we assess the proposed pruning methodology on six real-world datasets.

Problems and datasets
We address feature selection problems, which correspond to typical machine learning tasks.
Each task corresponds to an instance of the loss function $\lossfunc$ and two choices of $\pertfunc$, as described below.
Each dataset provides a matrix \(\dic \in \kR^{\ddim \times \pdim}\) and a vector \(\obs \in \kR^{\ddim}\).
The dimension of each dataset is specified in <Ref>.

First, we consider linear regression with the least-squares loss function:

\begin{equation} \label{eq:leastsquares-loss}
\lossfunc(\cdot) = \tfrac{1}{2}\norm{\obs - \cdot}{2}^2
\end{equation}

We use the instances of \(\obs\) and \(\dic\) provided by the Riboflavin [Bühlmann et al., 2014] and Bctcga [Liu et al., 2018] datasets which are related to vitamin production and cancer screening, respectively.

Second, we consider binary classification tasks with the logistic loss function:

\begin{equation} \label{eq:logistic-loss}
\lossfunc(\cdot) = \transp{\1}\log(\1 + \exp(-\obs \odot \cdot))
\end{equation}

where $\odot$ denotes the Hadamard product and the functions $\log(\cdot)$ and $\exp(\cdot)$ are taken component-wise.
We use instances of \(\obs\) and \(\dic\) from the Colon cancer [Alon et al., 1999] and Leukemia [Golub et al., 1999] datasets related to cancer screening.

Finally, we consider binary classification tasks with the squared hinge loss:

\begin{equation} \label{eq:squaredhinge-loss}
\lossfunc(\cdot) = \norm{\pospart{\1 - \obs \odot \cdot}}{2}^2
\end{equation}

where $\pospart{\cdot}$ is taken component-wise.
We use instances of \(\obs\) and \(\dic\) from the datasets Breast cancer and Arcene [Chang & Lin, 2011] related to DNA analysis and tumor categorization, respectively.

Dimensions of the datasets and data fidelity term \(\lossfunc\).

For each dataset, we consider the two following choices:

(·)    = α·       + (· ≤)

(·)    = α·^2       + (· ≤)

where $\alpha>0$ and $\bigM>0$.
These choices are motivated by the statistical properties of the solutions that can be obtained [Dedieu et al., 2021].

We mention that the  and  solvers can only handle linear and quadratic functions. Hence, they cannot address problem instances involving the logistic loss (<ref>).
Moreover,  can only handle instances combining functions (<ref>) and (<ref>).

Performance profiles

For each problem instance, we first calibrate the hyperparameters $\alpha>0$ and $\bigM>0$ as explained in <Ref>.
We then fit a regularization path [Friedman et al., 2010], that is, we vary the value of $\reg$ to construct solutions with different sparsity levels.
We start at some $\reg$ so that the all-zero vector is a solution of (<ref>) and sequentially decrease its value as long as at least one solver can solve the problem within one hour.
The solution obtained for each value of $\reg$ is used as a warm-start for the problem with the next value of $\reg$ considered in the regularization path.

<Ref> represents the time needed by each solver to construct a solution with a given sparsity level.
We observe that the implementation of the pruning methodology proposed in this paper allows for significant gains in terms of running time.
More precisely, our method outperforms the other solvers in all the considered scenarios.
In comparison with off-the-shelf mip solvers such as , , and , the time savings can reach up to four orders of magnitude in the most favorable cases.
Regarding the specialized solver , improvements of up to two orders of magnitude are achievable by our method in the best-case scenarios.

§ CONCLUSION

In this paper, we introduce a new methodology to perform pruning tests in Branch-and-Bound algorithms addressing $\ell_0$-regularized optimization problems.
Our method is only grounded on a few hypotheses and can thus be applied to a large variety of problems, notably in machine learning.
Our numerical results demonstrate that the proposed methodology significantly reduces solving time compared to other state-of-the-art methods.
It therefore allows to address some problem instances that were out of computational reach so far.

§ IMPACT STATEMENT

The goal of our work is to accelerate the solving time of some particular optimization problems.
Our contribution is primarily methodological and any of its potential societal impact would only be indirectly related to our work.

[Alon et al., 1999]
Alon, U., Barkai, N., Notterman, D. A., Gish, K., Ybarra, S., Mack, D., and
Levine, A. J.
Broad patterns of gene expression revealed by clustering analysis of
tumor and normal colon tissues probed by oligonucleotide arrays.
Proceedings of the National Academy of Sciences, 960
(12):0 6745–6750, 1999.

[Anand et al., 2017]
Anand, R., Aggarwal, D., and Kumar, V.
A comparative analysis of optimization solvers.
Journal of Statistics and Management Systems, 200
(4):0 623–635, 2017.

[Atamturk & Gómez, 2020]
Atamturk, A. and Gómez, A.
Safe screening rules for l0-regression from perspective relaxations.
In International Conference on Machine Learning, pp.  421–430. PMLR, 2020.

[Bauschke & Combettes, 2017]
Bauschke, H. H. and Combettes, P. L.
Convex Analysis and Monotone Operator Theory in Hilbert
Springer International Publishing, 2017.

[Beck, 2017]
Beck, A.
First-order methods in optimization.
SIAM, 2017.

[Ben Mhenni et al., 2020]
Ben Mhenni, R., Bourguignon, S., Mongeau, M., Ninin, J., and Carfantan, H.
Sparse branch-and-bound for exact optimization of l0-norm penalized
least squares.
In International Conference on Acoustics, Speech and Signal
Processing, pp.  5735–5739. IEEE, 2020.

[Ben Mhenni et al., 2022]
Ben Mhenni, R., Bourguignon, S., and Ninin, J.
Global optimization for sparse solution of least squares problems.
Optimization Methods and Software, 370 (5):0
1740–1769, 2022.

[Bertsimas & Cory-Wright, 2022]
Bertsimas, D. and Cory-Wright, R.
Solving large-scale sparse pca to certifiable (near) optimality.
The Journal of Machine Learning Research, 230
(1):0 566–600, 2022.

[Bertsimas & Johnson, 2023]
Bertsimas, D. and Johnson, N.
Compressed sensing: A discrete optimization approach.
arXiv preprint arXiv:2306.04647, 2023.

[Bertsimas & Van Parys, 2020]
Bertsimas, D. and Van Parys, B.
Sparse high-dimensional regression: Exact scalable algorithms and
phase transitions.
The Annals of Statistics, 480 (1):0 300–323,

[Bertsimas et al., 2016]
Bertsimas, D., King, A., and Mazumder, R.
Best subset selection via a modern optimization lens.
The Annals of Statistics, 440 (2):0 813–852,

[Bertsimas et al., 2020]
Bertsimas, D., Pauphilet, J., and Van Parys, B.
Sparse regression.
Statistical Science, 350 (4):0 555–578, 2020.

[Bertsimas et al., 2021]
Bertsimas, D., Cory-Wright, R., and Pauphilet, J.
A unified approach to mixed-integer optimization problems with
logical constraints.
SIAM Journal on Optimization, 310 (3):0
2340–2367, 2021.

[Bruer et al., 2015]
Bruer, J. J., Tropp, J. A., Cevher, V., and Becker, S. R.
Designing statistical estimators that balance sample size, risk, and
computational cost.
Journal of Selected Topics in Signal Processing, 90
(4):0 612–624, 2015.

[Bühlmann et al., 2014]
Bühlmann, P., Kalisch, M., and Meier, L.
High-dimensional statistics with a view toward applications in
Annual Review of Statistics and Its Application, 1:0
255–278, 2014.

[Candes et al., 2007]
Candes, E., Braun, N., and Wakin, M.
Sparse signal and image recovery from compressive samples.
In International Symposium on Biomedical Imaging: From Nano to
Macro, pp.  976–979. IEEE, 2007.

[Carreira-Perpinán & Idelbayev, 2018]
Carreira-Perpinán, M. A. and Idelbayev, Y.
“learning-compression” algorithms for neural net pruning.
In Conference on Computer Vision and Pattern Recognition, pp.   8532–8541. IEEE, 2018.

[Chang & Lin, 2011]
Chang, C.-C. and Lin, C.-J.
Libsvm: a library for support vector machines.
Transactions on Intelligent Systems and Technology, 20
(3):0 1–27, 2011.

[Dedieu et al., 2021]
Dedieu, A., Hazimeh, H., and Mazumder, R.
Learning sparse classifiers: Continuous and mixed integer
optimization perspectives.
The Journal of Machine Learning Research, 220
(1):0 6008–6054, 2021.

[Friedman et al., 2010]
Friedman, J., Hastie, T., and Tibshirani, R.
Regularization paths for generalized linear models via coordinate
Journal of Statistical Software, 330 (1):0 1,

[Golub et al., 1999]
Golub, T. R., Slonim, D. K., Tamayo, P., Huard, C., Gaasenbeek, M., Mesirov,
J. P., Coller, H., Loh, M. L., Downing, J. R., Caligiuri, M. A., et al.
Molecular classification of cancer: class discovery and class
prediction by gene expression monitoring.
Science, 2860 (5439):0 531–537, 1999.

[Guyard et al., 2022]
Guyard, T., Herzet, C., and Elvira, C.
Node-screening tests for the l0-penalized least-squares problem.
In International Conference on Acoustics, Speech and Signal
Processing, pp.  5448–5452. IEEE, 2022.

[Hastie et al., 2020]
Hastie, T., Tibshirani, R., and Tibshirani, R.
Best subset, forward stepwise or lasso? analysis and recommendations
based on extensive comparisons.
Statistical Science, 350 (4):0 579–592, 2020.

[Hazimeh et al., 2022]
Hazimeh, H., Mazumder, R., and Saab, A.
Sparse regression at scale: Branch-and-bound rooted in first-order
Mathematical Programming, 1960 (1-2):0
347–388, 2022.

[Liu et al., 2018]
Liu, J., Lichtenberg, T., Hoadley, K. A., Poisson, L. M., Lazar, A. J.,
Cherniack, A. D., Kovatich, A. J., Benz, C. C., Levine, D. A., Lee, A. V.,
et al.
An integrated tcga pan-cancer clinical data resource to drive
high-quality survival outcome analytics.
Cell, 1730 (2):0 400–416, 2018.

[Locatelli & Schoen, 2013]
Locatelli, M. and Schoen, F.
Global optimization: theory, algorithms, and applications.
SIAM, 2013.

[Nguyen et al., 2019]
Nguyen, T. T., Soussen, C., Idier, J., and Djermoune, E.-H.
Np-hardness of l0 minimization problems: revision and extension to
the non-negative setting.
In International Conference on Sampling Theory and
Applications, pp.  1–4. IEEE, 2019.

[Polson & Sokolov, 2019]
Polson, N. G. and Sokolov, V.
Bayesian regularization: From tikhonov to horseshoe.
Wiley Interdisciplinary Reviews: Computational Statistics,
110 (4):0 e1463, 2019.

[Rockafellar, 1967]
Rockafellar, R.
Duality and stability in extremum problems involving convex
Pacific Journal of Mathematics, 210 (1):0
167–187, 1967.

[Sarin et al., 1988]
Sarin, S., Karwan, M. H., and Rardin, R. L.
Surrogate duality in a branch-and-bound procedure for integer
European Journal of Operational Research, 330
(3):0 326–333, 1988.

[Tan et al., 2010]
Tan, M., Wang, L., and Tsang, I. W.
Learning sparse svm for feature selection on very high dimensional
In International Conference on Machine Learning, pp.  1047–1054, 2010.

[Tillmann et al., 2021]
Tillmann, A. M., Bienstock, D., Lodi, A., and Schwartz, A.
Cardinality minimization, constraints, and regularization: a survey.
arXiv preprint arXiv:2106.09606, 2021.

[Tropp & Wright, 2010]
Tropp, J. A. and Wright, S. J.
Computational methods for sparse solution of linear inverse problems.
Proceedings of the IEEE, 980 (6):0 948–958,

[Wang et al., 2020]
Wang, Q., Ma, Y., Zhao, K., and Tian, Y.
A comprehensive survey of loss functions in machine learning.
Annals of Data Science, pp.  1–26, 2020.

[Wolsey, 1980]
Wolsey, L. A.
Heuristic analysis, linear programming and branch and bound.
Combinatorial Optimization II, pp.  121–134, 1980.

[Zhong et al., 2022]
Zhong, P., Chen, T., Barroso-Luque, L., Xie, F., and Ceder, G.
An l0l2-norm regularized regression model for construction of robust
cluster expansions in multicomponent systems.
Physical Review B, 1060 (2):0 024203, 2022.

§ SUPPLEMENTARY MATERIAL RELATED TO SECTIONS <REF> AND <REF>

This section gathers discussions and proofs of the results presented in <Ref> of the paper.

§.§ Discussion on Section <ref>

In this paragraph, we give additional details on the relation between the tightness of lower-bound (<ref>) and the depth of the node at which it is computed.
Our claim is grounded on the two results stated in the following lemma:

Let \(\nodeSymb=(\setone,\setzero,\setnone)\) be a node of the bnb tree.

* For all successors \(\nodeSymb'\) of \(\nodeSymb\), we have \(\robj^{\nodeSymb'} \geq \robj^{\nodeSymb}\).

* If <ref>-<ref> hold and \(\setnone=\emptyset\) then \(\node{\robj} = \node{\pobj}\).

We prove the two items separately.

* Let \(\nodeSymb'=(\setone',\setzero',\setnone')\) be a successor of \(\nodeSymb\).
By definition of a successor, we observe from (<ref>) that \(\pset^{\nodeSymb'} \subseteq \node{\pset}\).
Therefore, using the definitions of \({\regfunc}^{\nodeSymb}\) and \({\regfunc}^{\nodeSymb'}\) in (<ref>), we have

\begin{equation} \label{eq:proof bounds get tighter:central inequality}
\geq
\end{equation}

Item (ii) of Proposition 13.16 in [Bauschke & Combettes, 2017] then leads to

\begin{equation}
\biconj{(\regfunc^{\nodeSymb'})}(\cdot)
\geq
\biconj{(\node{\regfunc})}(\cdot)
\end{equation}

Thus, we have the inequality \(\lossfunc(\dic\cdot) + \biconj{(\regfunc^{\nodeSymb'})}(\cdot) \geq \lossfunc(\dic\cdot) + \biconj{(\regfunc^{\nodeSymb})}(\cdot)\).
By taking the infimum on both sides of this inequality, we obtain the desired result.

* If <ref>-<ref> hold and \(\setnone=\emptyset\), one has \(\biconj{(\node{\regfunc})} = \node{\regfunc}\) from <Ref> and Theorem 4.8 in [Beck, 2017]. Hence, the result directly follows.

We now motivate our claim “achieving the prescribed tightness usually imposes exploring deep nodes in the decision tree” in light of the following analysis.

Let \(\UB{\pobj}\) be an upper bound on \(\opt{\pobj}\) and \(\initialnode=(\setone^{(0)},\setzero^{(0)},\setnone^{(0)})\) be a node[
Contrary to the notational convention used in <Ref>, \(\initialnode\) does not necessarily refers (here) to the root node \((\emptyset,\emptyset,\intervint{1}{\pdim})\).
] of the bnb tree satisfying \(\setnone^{(0)}\neq\emptyset\) and

\begin{equation} \label{eq:appendix discussion:node passes ideal test}
\pobj^{\initialnode}
\UB{\pobj}
\geq
\robj^{\initialnode}
\end{equation}

Denote \(\lengthnone=\card{\setnone^{(0)}}\).
The first inequality guarantees that \(\pset^{\nodeSymb_0}\) definitively excludes any minimizer of (<ref>), while the second indicates that the pruning test (<ref>) with the lower bound (<ref>) is unsuccessful.
We nevertheless demonstrate that  for any sequence \(\{\sequencenode{\idxdepth}\}_{\idxdepth=1}^{\lengthnone}\)

\begin{equation}\label{eq:cond-node-sequence}
\forall\idxdepth\in\intervint{0}{\lengthnone-1}:\quad
\sequencenode{\idxdepth+1} \text{ is a direct successor of } \sequencenode{\idxdepth}
\end{equation}

there exists an index \(\criticalidx\in\intervint{1}{\lengthnone}\)
such that \(\sequencenode{\criticalidx}\) passes test (<ref>) with the lower bound (<ref>).

Let \(\{\sequencenode{\idxdepth}\}_{\idxdepth=1}^{\lengthnone}\) be a sequence verifying (<ref>).
First, using item <ref> of <Ref>, we have that $\delta_\idxdepth \triangleq \robj^{\sequencenode{\idxdepth}} - \UB{\pobj}$ is a non-decreasing function of $\idxdepth$.
Note that the second inequality in (<ref>) leads to \(\delta_0 \leq 0\).
Moreover, we have that \(\delta_{\lengthnone} > 0\) by virtue of the following arguments: first, the definition of a successor implies that \(\pset^{\sequencenode{\lengthnone}} \subsetneq \pset^{\initialnode}\) and therefore \(\pobj^{\sequencenode{\lengthnone}} \geq \pobj^{\initialnode}\); second, we have from item <ref> of <Ref> that \(\pobj^{\nodeSymb_{\lengthnone}} = \robj^{\sequencenode{\lengthnone}}\); finally the result follows from the first inequality in (<ref>).

In conclusion, since \(\delta_{\lengthnone} > 0\), the pruning test (<ref>) with the lower bound (<ref>) will ultimately succeed if enough successors are visited by the bnb procedure.
Furthermore, the likelihood of the success increases as the bnb procedure explores deeper nodes in the tree because \(\delta_\idxdepth\) is a non-decreasing function of \(\idxdepth\).

§.§ Proof of Proposition <ref>

We first prove the following technical lemma:

Let $\kfuncdef{\prooffun}{\kR}{\kR\cup\{+\infty\}}$ be a closed, convex, proper function and
let \(\pvi{0}\in\dom(\prooffun)\) be an accumulation point of \(\dom(\prooffun)\).
Then we have

\begin{equation}
\inf_{\pvi{}\in\kR\setminus\{\pvi{0}\}} \prooffun(\pvi{})
\inf_{\pvi{}\in\kR} \prooffun(\pvi{}).
\end{equation}

We obviously have

\begin{equation}
\inf_{\pvi{}\in\kR\setminus\{\pvi{0}\}} \prooffun(\pvi{})
\geq
\inf_{\pvi{}\in\kR} \prooffun(\pvi{})
\end{equation}

so we concentrate on the reverse inequality hereafter.

First, using the fact that \(\pvi{0}\) is an accumulation point of \(\dom(\prooffun)\), we have that \(\dom(\prooffun)\setminus\{\pvi{0}\}\neq \emptyset\) and therefore

\begin{align}
\inf_{\pvi{}\in\kR} \prooffun(\pvi{}) &=  \inf_{\pvi{}\in\dom(\prooffun)} \prooffun(\pvi{})
\\
\inf_{\pvi{}\in\kR\setminus\{\pvi{0}\}} \prooffun(\pvi{}) &=  \inf_{\pvi{}\in\dom(\prooffun)\setminus\{\pvi{0}\}} \prooffun(\pvi{})
\end{align}

It is thus sufficient to prove that

\begin{equation} \label{eq:proof lemma infima=:target inequality}
\inf_{\pvi{}\in\dom(\prooffun)\setminus\{\pvi{0}\}} \prooffun(\pvi{})
\leq
\inf_{\pvi{}\in\dom(\prooffun)} \prooffun(\pvi{})
\end{equation}

Second, using the fact that \(\pvi{0}\in\dom(\prooffun)\) by hypothesis, we also have

\begin{equation} \label{eq:proof lemma infima=:infima in not +infty}
\inf_{\pvi{} \in\dom(\prooffun)} \prooffun(\pvi{})
\leq \prooffun(\pvi{0})
< +\infty
\end{equation}

We then prove (<ref>) by considering two separate cases.

* Assume first that

\begin{equation} \label{eq:proof lemma infima=:infima in not +infty:case 1}
\inf_{\pvi{} \in\dom(\prooffun)} \prooffun(\pvi{})
\prooffun(\pvi{0})
\end{equation}

Since \(\pvi{0}\) is an accumulation point of \(\dom(\prooffun)\), there exists a sequence \(\{\pvi{}^{(\idxsequence)}\}_{\idxsequence\in\kN}\subset\dom(\prooffun)\setminus\{\pvi{0}\}\).[

Such a sequence can be constructed as follows: for all \(\idxsequence\in\kN\), one chooses \(\pvi{}^{(\idxsequence)}\in(\proofset\cap B(\pvi{0}, \kinv{(\idxsequence+1)}))\setminus\{\pvi{0}\}\) which is nonempty by definition of an accumulation point.
From Th. 2.22 in [Beck, 2017], we have that closedness, convexity and properness of \(\prooffun\) implies that it is continuous on its domain.[
That is for any $\{x^{(\idxsequence)}\}_{\idxsequence\in\kN} \subset \dom(\prooffun)$ converging to some limit point $\pvi{0}\in\dom(\prooffun)$, we have $\lim_{\idxsequence\rightarrow+\infty} \prooffun(\pvi{}^{(\idxsequence)})=\prooffun(\pvi{0})$.
Hence \(\lim_{\idxsequence\rightarrow+\infty} \prooffun(\pvi{}^{(\idxsequence)})=\prooffun(\pvi{0})\) and therefore

\begin{equation} \label{eq:proof lemma infima=:infima in not +infty:case 1 buf}
\inf_{\pvi{}\in\dom(\prooffun)\setminus\{\pvi{0}\}}  \prooffun(\pvi{})
\leq
\prooffun(\pvi{0})
\end{equation}

Inequality (<ref>) immediately follows by combining (<ref>) and (<ref>).

* Assume now that

\begin{equation} \label{eq:proof lemma infima=:infima in not +infty:case 2}
\inf_{\pvi{} \in\dom(\prooffun)} \prooffun(\pvi{})
\prooffun(\pvi{0})
\end{equation}

and denote \(\prooffun_{\inf}\in\kR\cup\{-\infty\}\) the latter infimum.
By definition of an infimum, there exists a sequence \(\{\proofsequence{\idxsequence}\}_{\idxsequence\in\kN}\subset\dom(\prooffun)\) such that \(\prooffun_{\inf} = \lim_{\idxsequence\rightarrow+\infty} \prooffun(\proofsequence{\idxsequence})\).
Let \(\varepsilon>0\) be such that \(\varepsilon<\prooffun(\pvi{0}) - \prooffun_{\inf}\).
The hypothesis case (<ref>) therefore implies that there exists \(\idxsequence_{\varepsilon}\) such that

\begin{equation}
\forall\idxsequence\in\kN,\qquad
\idxsequence\geq\idxsequence_\varepsilon \Longrightarrow
\prooffun(\proofsequence{\idxsequence}) < \prooffun(\pvi{0}) - \varepsilon
\end{equation}

and therefore \(\proofsequence{\idxsequence}\neq\pvi{0}\) for all \(\idxsequence\geq\idxsequence_\varepsilon\).
On can thus construct a subsequence \(\{\proofsubsequence{\idxsequence}\}_{\idxsequence\in\kN}\) such that
\(\proofsubsequence{\idxsequence}\neq\pvi{0}\) for all \(\idxsequence\in\kN\).
This implies that \(\{\proofsubsequence{\idxsequence}\}_{\idxsequence\in\kN}\subset\dom(\prooffun)\setminus\{\pvi{0}\}\) and one immediately deduces that

\begin{equation}
\inf_{\pvi{}\in\dom(\prooffun)\setminus\{\pvi{0}\}}
\leq
\lim_{\idxsequence\to+\infty} \prooffun(\proofsubsequence{\idxsequence}) =
\prooffun_{\inf}
\end{equation}

where the equality holds since \(\{\prooffun(\proofsubsequence{\idxsequence})\}_{\idxsequence\in\kN}\) is a subsequence of a converging sequence.
This leads to (<ref>).

We are now ready to give a proof to <Ref>.
By definition of a convex conjugate function, we have:

\begin{equation}
\forall\vv \in \kR^{\pdim}:\ \conj{(\node{\regfunc})}(\vv) = \textstyle\sup_{\pv \in \kR^{\pdim}} \transp{\vv}\pv - \node{\regfunc}(\pv)
\end{equation}

Observing from (<ref>) that $\node{\regfunc}(\cdot)$ is separable since both $\regfunc(\cdot)$ and $\icvx(\cdot \in \node{\pset})$ are separable (see their definitions in (<ref>) and (<ref>), respectively), we deduce that $\conj{(\node{\regfunc})}(\cdot)$ is also separable and is given coordinate-wise by

\begin{align}
\conj{(\separable{\node{\regfunc}}{\idxentry})}(\vvi{}) &= \textstyle\sup_{\pvi{} \in \kR} \vvi{}\pvi{} - \separable{\node{\regfunc}}{\idxentry}(\pvi{})
\end{align}

Imposing explicitly the constraints defined in $\node{\pset}$, we obtain:

\begin{align}
\conj{(\separable{\node{\regfunc}}{\idxentry})}(\vvi{}) &=
\begin{cases}
\textstyle\sup_{\pvi{} = 0} \, \vvi{}\pvi{} - \pertfunc(\pvi{}) - \reg\norm{\pvi{}}{0} &\text{if} \ \idxentry \in \setzero \\
\textstyle\sup_{\pvi{} \neq 0} \, \vvi{}\pvi{} - \pertfunc(\pvi{}) - \reg\norm{\pvi{}}{0} &\text{if} \ \idxentry \in \setone \\
\textstyle\sup_{\pvi{} \in \kR} \vvi{}\pvi{} - \pertfunc(\pvi{}) - \reg\norm{\pvi{}}{0} &\text{if} \ \idxentry \in \setnone
\end{cases}
\label{eq:proof:dual-relaxregfunc:2}
\end{align}

We next address the three above cases separately.

If \(\idxentry \in \setzero\), the first case in (<ref>) simplifies to

\begin{equation} \label{eq:proof:dual-relaxregfunc:3}
\conj{(\separable{\node{\regfunc}}{\idxentry})}(\vvi{}) = \textstyle\sup_{\pvi{} = 0} \, \vvi{}\pvi{} - \pertfunc(\pvi{}) = 0
\end{equation}

where the first equality holds since \(\norm{0}{0}=0\) and the second since $\pertfunc(0)=0$ in virtue of hypothesis <ref>.

If $\idxentry \in \setone$, then

\begin{align}
\label{eq:proof:dual-relaxregfunc:4}
\conj{(\separable{\node{\regfunc}}{\idxentry})}(\vvi{}) &= \textstyle\sup_{\pvi{} \neq 0} \ \, \vvi{}\pvi{} - \pertfunc(\pvi{}) - \reg
\end{align}

since \(\norm{\pvi{}}{0}=1\) for all \(\pvi{}\neq0\).
As $\pertfunc(\cdot)$ is closed, convex and proper, the function $\prooffun(\pvi{})\triangleq -\vvi{}\pvi{} + \pertfunc(\pvi{}) + \reg$ inherits from these properties. On the one hand, it is easy to see that $\dom (\prooffun)=\dom(\pertfunc)$.
On the other hand, since $0$ is an accumulation point of $\dom(\pertfunc)$ from <ref>, we have from <Ref> that

\begin{align}
\conj{(\separable{\node{\regfunc}}{\idxentry})}(\vvi{})
\label{eq:proof:dual-relaxregfunc:5}
&= \textstyle\sup_{\pvi{} \in \kR} \, \vvi{}\pvi{} - \pertfunc(\pvi{}) - \reg.
\end{align}

We finally obtain the result by using the definition of the convex conjugate of $\pertfunc$.

If $\idxentry \in \setnone$, then

()()    = sup_ ∈    - ()  - 0

= max{
sup_ = 0    - () - 0        ,             sup_ ≠0    - () - 0

= max{
sup_ = 0    - ()        ,             sup_ ≠0    - () -       }

= max{0, () - }

= () -

where (<ref>) is obtained by definition of the $\ell_0$-norm, (<ref>) follows from the same reasoning as for the case “$\idxentry\in\setone$”.

§.§ Proof of Proposition <ref>

Let \(\nodeSymb'=(\setzero',\setone',\setnone')\) be a successor of \(\nodeSymb=(\setzero,\setone,\setnone)\) and \(\dv\in\kR^\ddim\).
<Ref> is a direct consequence of the following two inequalities

\begin{align}
\dfunc^{\nodeSymb'_{0,\idxentry}}(\dv) \,\geq\,& \dfunc^{\succnode{0}{\idxentry}}(\dv) \label{eq:proof propagation test:corner inequality a}\\
\dfunc^{\nodeSymb'_{1,\idxentry}}(\dv) \,\geq\,& \dfunc^{\succnode{1}{\idxentry}}(\dv) \label{eq:proof propagation test:corner inequality b}.
\end{align}

We thus establish (<ref>) and (<ref>) in the remaining of the section.

Note first that a direct consequence of <Ref> is

\begin{equation}\label{eq:increase-dual-fun}
\dfunc^{\nodeSymb'}(\dv)\geq \dfunc^{\nodeSymb}(\dv)
\end{equation}

since all terms \(\{\pivot{0}(\transp{\atom{\idxentry}}\dv)\}_{\idxentry\in\setzero'\setminus \setzero}\) and \(\{\pivot{1}(\transp{\atom{\idxentry}}\dv)\}_{\idxentry\in\setone'\setminus \setone}\) are nonnegative.

Let \(\binarynumber\in\{0, 1\}\).
Particularizing <Ref> to \(\nodeSymb=\nodeSymb'\) and $\nodeSymb'=\nodeSymb'_{\binarynumber,\idxentry}$ –the direct successor of $\nodeSymb'$ defined in (<ref>)– one obtains:

= ^'()    +

≥^()    +

where the inequality follow from (<ref>) (since $\nodeSymb'$ is a successor of $\nodeSymb$) and the last equality from <Ref>. This establishes the result.

§ SUPPLEMENTARY MATERIALS RELATED TO SECTION <REF>

This section gives supplementary materials to our numerical experiments.

§.§ Mixed-integer Programming Formulations

In our experiments, problem (<ref>) is formulated as a mip so that it can be handled by commercial solvers like ,  and .
For the problem considered in <Ref> where $\lossfunc(\cdot)$ and $\pertfunc(\cdot)$ are given by (<ref>)-(<ref>), we use the following formulation

\begin{equation}
\left\{
\begin{array}{rl}
\min & \tfrac{1}{2}\norm{\obs-\dic\pv}{2}^2 + \reg\transp{\1}\bv \\
\text{s.t.} & -\bigM \bv \leq \pv \leq \bigM \bv \\
& \pv \in \kR^{\pdim}, \ \bv \in \{0,1\}^{\pdim}
\end{array}
\right.
\end{equation}

where an additional binary variable $\bv \in \{0,1\}^{\pdim}$ is used to encode the nullity of the entries of the continuous variable $\pv \in \kR^{\pdim}$.
A similar approach can be used to reformulate the problems treated in <Ref>.
We refer the reader to <Ref> for the mip formulation of each of the problems considered in our numerical simulations.

1c|Penalty (<ref>)     1cPenalty (<ref>)

[origin=c]90Loss (<ref>)

\left\{
\begin{array}{rl}
\min & \tfrac{1}{2}\norm{\obs-\dic\pv}{2}^2 + \reg\transp{\1}\bv + \alpha\transp{\1}\mathbf{s} \\
\text{s.t.} & \pv \geq -\mathbf{s} \\
& \pv \leq \mathbf{s} \\
& \pv \geq -\bigM \bv \\
& \pv \leq \bigM \bv \\
& \pv \in \kR^{\pdim}, \ \bv \in \{0,1\}^{\pdim}, \ \mathbf{s} \in \kR^{\pdim}
\end{array}
\right.

\left\{
\begin{array}{rl}
\min & \tfrac{1}{2}\norm{\obs-\dic\pv}{2}^2 + \reg\transp{\1}\bv + \alpha\transp{\1}\mathbf{s} \\
\text{s.t.} & \pv\odot\pv \leq \mathbf{s}\odot\bv \\
& \pv \geq -\bigM \bv \\
& \pv \leq \bigM \bv \\
& \pv \in \kR^{\pdim}, \ \bv \in \{0,1\}^{\pdim}, \ \mathbf{s} \in \kR^{\pdim}
\end{array}
\right.

[origin=c]90Loss (<ref>)
\left\{
\begin{array}{rll}
\min & \transp{\1}\mathbf{u} + \reg\transp{\1}\bv + \alpha\transp{\1}\mathbf{s} \\
\text{s.t.} & \1 \geq \mathbf{v} + \mathbf{w} \\
& \mathbf{u} \geq -\log(\mathbf{v}) + \obs \odot \dic\pv \\
& \mathbf{u} \geq -\log(\mathbf{w}) \\
& \pv \geq -\mathbf{s} \\
& \pv \leq \mathbf{s} \\
& \pv \geq -\bigM \bv \\
& \pv \leq \bigM \bv \\
& \pv \in \kR^{\pdim}, \ \bv \in \{0,1\}^{\pdim}, \ \mathbf{s} \in \kR^{\pdim} \\
& \mathbf{u} \in \kR^{\ddim}, \ \mathbf{v} \in \kR^{\ddim}, \ \mathbf{w} \in \kR^{\ddim}
\end{array}
\right.

\left\{
\begin{array}{rll}
\min & \transp{\1}\mathbf{u} + \reg\transp{\1}\bv + \alpha\transp{\1}\mathbf{s} \\
\text{s.t.} & \mathbf{v} + \mathbf{w} \leq \1 \\
& \mathbf{u} \geq -\log(\mathbf{v}) + \obs \odot \dic\pv \\
& \mathbf{u} \geq -\log(\mathbf{w}) \\
& \pv \odot \pv \leq \mathbf{s} \odot \bv \\
& -\bigM \bv \leq \pv \leq \bigM \bv \\
& \pv \in \kR^{\pdim}, \ \bv \in \{0,1\}^{\pdim}, \ \mathbf{s} \in \kR^{\pdim} \\
& \mathbf{u} \in \kR^{\ddim}, \ \mathbf{v} \in \kR^{\ddim}, \ \mathbf{w} \in \kR^{\ddim}
\end{array}
\right.

[origin=c]90Loss (<ref>)
\left\{
\begin{array}{rll}
\min & \norm{\mathbf{w}}{2}^2 + \reg\transp{\1}\bv + \alpha\transp{\1}\mathbf{s} \\
\text{s.t.} & \mathbf{w} \geq \1 - \obs \odot \dic\pv \\
& \mathbf{w} \geq \0 \\
& \pv \geq -\mathbf{s} \\
& \pv \leq \mathbf{s} \\
& \pv \geq -\bigM \bv \\
& \pv \leq \bigM \bv \\
& \pv \in \kR^{\pdim}, \ \bv \in \{0,1\}^{\pdim} \\
& \mathbf{s} \in \kR^{\pdim}, \mathbf{w} \in \kR^{\pdim}
\end{array}
\right.

\left\{
\begin{array}{rll}
\min & \norm{\mathbf{w}}{2}^2 + \reg\transp{\1}\bv + \alpha\transp{\1}\mathbf{s} \\
\text{s.t.} & \mathbf{w} \geq \1 - \obs \odot \dic\pv \\
& \mathbf{w} \geq \0 \\
& \pv \odot \pv \leq \mathbf{s} \odot \bv  \\
& \pv \geq -\bigM \bv \\
& \pv \leq \bigM \bv \\
& \pv \in \kR^{\pdim}, \ \bv \in \{0,1\}^{\pdim} \\
& \mathbf{s} \in \kR^{\pdim}, \mathbf{w} \in \kR^{\pdim}
\end{array}
\right.

MIP formulations used Section 4. The vectorial inequalities as well as the function $\log(\cdot)$ are taken component-wise and $\odot$ denotes the Hadamard product.

§.§ Implementation Choices

Our bnb solver follows the standard implementation specified in <Ref> and explores the tree in a “depth-first” fashion, as presented in Sec. 5.2.2 by [Locatelli & Schoen, 2013].
Given some node $\nodeSymb=(\setzero,\setone,\setnone)$ where the decision tree must be expanded, we select an index $\idxentry \in \setnone$ to create new nodes (<ref>)-(<ref>) as

\begin{equation}
\label{eq:node-selection}
\idxentry \in \argmax_{\idxentry' \in \setnone} \abs{\node{\pvi{\idxentry'}}}
\end{equation}

where $\node{\pv}$ is the final iterate of the numerical procedure addressing (<ref>).
We use an approach similar to that considered in [Hazimeh et al., 2022] to solve numerically problem (<ref>).
More precisely, we solve a sequence of sub-problems defined as

\begin{equation}
\stepcounter{equation}
\tag{$\node{\rpb}_{\mathcal{W}}$}
\label{prob:relax-node-subproblem}
\left\{
\begin{array}{rl}
\inf &\lossfunc(\dic\pv) +
\biconj{(\node{\regfunc})}(\pv)
\\
\text{s.t.} & \pvi{\idxentry} = 0 \quad \forall \idxentry \notin \mathcal{W}
\end{array}
\right.
\end{equation}

where $\mathcal{W} \subseteq \intervint{1}{\pdim}$ is some  working set.
Each sub-problem (<ref>) is solved by using a classical coordinate descent procedure.
The size of the working set is $\mathcal{W}$ increased based on Fermat's optimality condition violation until no more violations occur.
More specifically, letting $\pv^{\mathcal{W}} \in \kR^{\pdim}$ be the output of the numerical procedure addressing (<ref>), we let $\mathcal{W}_{\mathrm{new}} = \mathcal{W} \cup \kset{\idxentry}{0 \notin \transp{\atom{\idxentry}}\subdiff\lossfunc(\dic\pv^{\mathcal{W}}) + \subdiff{\biconj{(\node{\separable{\regfunc}{\idxentry}})}}(\pvi{\idxentry}^{\mathcal{W}})}$ and stop the procedure when $\mathcal{W}_{\mathrm{new}} = \mathcal{W}$.

The simultaneous pruning procedure proposed in <Ref> is performed during the resolution of problems (<ref>).
In view of our discussion in <Ref>, we consider that one iteration of the solving process of (<ref>) corresponds to the resolution of one sub-problem (<ref>).
Stated otherwise, the iterates $\hat{\pv} \in \kR^{\pdim}$ used to implement our pruning methodology in (<ref>) are the output of the numerical procedure addressing each sub-problem (<ref>).

§.§ Technical Implementation Details

Given some node $\nodeSymb = (\setzero, \setone, \setnone)$, the bnb algorithm requires characterizing the convex biconjugate $\biconj{(\node{\regfunc})}(\cdot)$ associated with $\node{\regfunc}(\cdot)$ to construct relaxation (<ref>).
We derive its expression from the parametrization of the convex conjugate $\conj{(\node{\regfunc})}(\cdot)$ given in <Ref>.
With equation (<ref>), we first observe that $\biconj{(\node{\regfunc})}(\pv)=\sum_{\idxentry=1}^{\pdim}\biconj{(\node{\separable{\regfunc}{\idxentry}})}(\pvi{\idxentry})$ where

\begin{align}
\label{eq:supp_numerics:implementation:1:1}
\biconj{(\separable{\node{\regfunc}}{\idxentry})}(\pvi{}) &=
\begin{cases}
\textstyle\sup_{\vvi{} \in \kR} \pvi{}\vvi{} - 0 &\text{if} \ \idxentry \in \setzero \\
\textstyle\sup_{\vvi{} \in \kR} \pvi{}\vvi{} - (\conj{\pertfunc}(\vvi{}) - \reg) &\text{if} \ \idxentry \in \setone \\
\textstyle\sup_{\vvi{} \in \kR} \pvi{}\vvi{} - \pospart{\conj{\pertfunc}(\vvi{}) - \reg} &\text{if} \ \idxentry \in \setnone
\end{cases} \\
\label{eq:supp_numerics:implementation:1:2}
\begin{cases}
\icvx(\pvi{} = 0) &\text{if} \ \idxentry \in \setzero \\
\biconj{\pertfunc}(\pvi{}) + \reg &\text{if} \ \idxentry \in \setone \\
\textstyle\sup_{\vvi{} \in \kR} \pvi{}\vvi{} - \pospart{\conj{\pertfunc}(\vvi{}) - \reg} &\text{if} \ \idxentry \in \setnone
\end{cases} \\
\label{eq:supp_numerics:implementation:1:3}
\begin{cases}
\icvx(\pvi{} = 0) &\text{if} \ \idxentry \in \setzero \\
\pertfunc(\pvi{}) + \reg &\text{if} \ \idxentry \in \setone \\
\textstyle\sup_{\vvi{} \in \kR} \pvi{}\vvi{} - \pospart{\conj{\pertfunc}(\vvi{}) - \reg} &\text{if} \ \idxentry \in \setnone
\end{cases}
\end{align}

Equalities (<ref>)-(<ref>) follow from Definition 4.1 in [Beck, 2017].
Equality (<ref>) follows from Theorem 4.8 in [Beck, 2017] since $\pertfunc(\cdot)$ is proper, closed and convex under hypothesis <ref>.
We next provide a closed-form expression of the case $\idxentry \in \setnone$ in (<ref>) for the three expressions of function $\pertfunc(\cdot)$ considered in our numerical experiments.

Function $\pertfunc(\cdot)$ given by (<ref>)
With this parametrization, we have

\begin{align}
\conj{\pertfunc}(\vvi{}) &= \textstyle\sup_{\pvi{} \in \kR} \vvi{}\pvi{} - \icvx(\abs{\pvi{}} \leq \bigM) \\
&= \bigM\abs{\vvi{}}
\end{align}

Hence, we deduce that $\conj{\pertfunc}(\vvi{}) \leq \reg \iff \abs{\vvi{}} \leq \reg/\bigM$, which gives

sup_ ∈  - () -     = max{
sup_ ≤/         ,
sup_ ≥/  - ( - )

= max{
sup_ ≥/  -  +       }

Further, we remark that the supremum in the right member of (<ref>) lies in $\kset{\vvi{} \in \kR \cup \{+\infty\}}{\vvi{} \geq 0}$ when $\pvi{} \geq 0$ and in $\kset{\vvi{} \in \kR \cup \{-\infty\}}{\vvi{} \leq 0}$ when $\pvi{} \leq 0$.
Therefore, we obtain

sup_ ≥/  -  +    = sup_ ≥/ ( - ) +

+∞   if   >

(/)    if   ≤

= (/) + ( ≤)

Overall, combining (<ref>)-(<ref>) with (<ref>)-(<ref>) gives

\begin{align}
\biconj{(\node{\separable{\regfunc}{\idxentry}})}(\pvi{}) &= \textstyle\sup_{\vvi{} \in \kR} \pvi{}\vvi{} - \pospart{\conj{\pertfunc}(\vvi{}) - \reg} \\
&= \max\{(\reg/\bigM)\abs{\pvi{}}, \ (\reg/\bigM)\abs{\pvi{}} + \icvx(\abs{\pvi{}} \leq \bigM)\} \\
&= (\reg/\bigM)\abs{\pvi{}} + \icvx(\abs{\pvi{}} \leq \bigM)
\end{align}

when $\idxentry \in \setnone$ in (<ref>).

Function $\pertfunc(\cdot)$ given by (<ref>)

With this parametrization, we have

\begin{align}
\conj{\pertfunc}(\vvi{}) &= \textstyle\sup_{\abs{\pvi{}} \leq \bigM} \vvi{}\pvi{} - \alpha\abs{\pvi{}} \\
&= \textstyle\sup_{\abs{\pvi{}} \leq \bigM} \abs{\pvi{}}(\abs{\vvi{}} - \alpha) \\
&= \begin{cases}
0 &\text{if} \ \abs{\vvi{}} \leq \alpha \\
\bigM(\abs{\vvi{}} - \alpha) &\text{if} \ \abs{\vvi{}} > \alpha
\end{cases} \\
\label{eq:proof:relaxregfunc:6}
&= \bigM\pospart{\abs{\vvi{}} - \alpha}
\end{align}

Hence, we deduce that $\conj{\pertfunc}(\vvi{}) \leq \reg \iff \abs{\vvi{}} \leq \alpha+\reg/\bigM$, which gives

sup_ ∈  - () -
= max{
sup_ ≤α+ /         ,
sup_ ≥α+ /  -  + α+       }

= max{
(α+ /)        ,
sup_ ≥α+ /  -  + α+       }

Further, we remark that the supremum of the right member in (<ref>) lies in $\kset{\vvi{} \in \kR \cup \{+\infty\}}{\vvi{} \geq 0}$ when $\pvi{} \geq 0$ and in $\kset{\vvi{} \in \kR \cup \{-\infty\}}{\vvi{} \leq 0}$ when $\pvi{} \leq 0$.
Therefore, we obtain

\begin{align}
\label{eq:proof:relaxregfunc:7:4}
\textstyle\sup_{\abs{\vvi{}} \geq \alpha + \reg/\bigM} \pvi{}\vvi{} - \bigM\abs{\vvi{}} + \bigM\alpha + \reg &= \textstyle\sup_{\abs{\vvi{}} \geq \alpha + \reg/\bigM} \abs{\vvi{}}(\abs{\pvi{}} - \bigM) + \bigM\alpha + \reg \\
\begin{cases}
+\infty &\text{if} \ \abs{\pvi{}} > \bigM \\
(\alpha + \reg/\bigM)\abs{\pvi{}} &\text{if} \ \abs{\pvi{}} \leq \bigM
\end{cases} \\
\label{eq:proof:relaxregfunc:7:5}
&= (\alpha + \reg/\bigM)\abs{\pvi{}} + \icvx(\abs{\pvi{}} \leq \bigM)
\end{align}

Overall, combining (<ref>)-(<ref>) with (<ref>)-(<ref>) gives

\begin{align}
\biconj{(\node{\separable{\regfunc}{\idxentry}})}(\pvi{}) &= \textstyle\sup_{\vvi{} \in \kR} \pvi{}\vvi{} - \pospart{\conj{\pertfunc}(\vvi{}) - \reg} \\
&= \max\big\{
(\alpha + \reg/\bigM)\abs{\pvi{}}, \
(\alpha + \reg/\bigM)\abs{\pvi{}} + \icvx(\abs{\pvi{}} \leq \bigM)
\big\} \\
&= (\alpha + \reg/\bigM)\abs{\pvi{}} + \icvx(\abs{\pvi{}} \leq \bigM)
\end{align}

when $\idxentry \in \setnone$ in (<ref>).

Function $\pertfunc(\cdot)$ given by (<ref>).
With this parametrization, we obtain

\begin{align}
\conj{\pertfunc}(\vvi{})
&\triangleq \textstyle\sup_{\abs{\pvi{}} \leq \bigM} \vvi{}\pvi{} - \alpha\pvi{}^2 \\
&= -\textstyle\inf_{\abs{\pvi{}} \leq \bigM} \alpha\pvi{}^2 - \vvi{}\pvi{}
\label{eq:proof:relaxregfunc:square abs + interval:buf conjugate} \\
&= \begin{cases}
\tfrac{\vvi{}^2}{4\alpha} &\text{if} \ \abs{\vvi{}} \leq 2\alpha\bigM \\
\bigM\abs{\vvi{}} - \alpha\bigM^2 &\text{if} \ \abs{\vvi{}} > 2\alpha\bigM
\end{cases}
\end{align}

by noting that the scalar defined as the orthogonal projection   of \(\tfrac{\vvi{}}{2\alpha}\) onto the interval \([-\bigM, \bigM]\) satisfies the necessary and sufficient first order optimality condition (see e.g. Corollary 3.68 in [Beck, 2017]) associated to the convex minimization problem involved in (<ref>).
The function $\conj{\pertfunc}(\cdot)$ is continuous, monotone, minimized at $\vvi{} = 0$ and one has $\conj{\pertfunc}(\vvi{}) = \alpha\bigM^2$ at the threshold $\abs{\vvi{}} = 2\alpha\bigM$.
In this view, we deduce that

\begin{align}
\conj{\pertfunc}(\vvi{}) \leq \reg &\iff
\begin{cases}
\bigM\abs{\vvi{}} - \alpha\bigM^2 \leq \reg &\text{if} \ \alpha\bigM^2 \leq \reg  \\
\tfrac{\vvi{}^2}{4\alpha} \leq \reg &\text{if} \ \alpha\bigM^2 > \reg
\end{cases}
\\
\label{eq:proof:relaxregfunc:8}
&\iff \begin{cases}
\abs{\vvi{}} \leq \alpha\bigM + \reg/\bigM & \text{if} \ \bigM \leq \sqrt{\reg/\alpha} \\
\abs{\vvi{}} \leq 2\sqrt{\reg\alpha} &\text{if} \ \bigM > \sqrt{\reg/\alpha}
\end{cases}
\end{align}

By treating the two above cases separately, we next show that

\begin{align}
\biconj{(\node{\separable{\regfunc}{\idxentry}})}(\pvi{}) &= \textstyle\sup_{\vvi{} \in \kR} \pvi{}\vvi{} - \pospart{\conj{\pertfunc}(\vvi{}) - \reg} \\
\begin{cases}
\icvx(\abs{\pvi{}} \leq \bigM) + (\tfrac{\reg}{\bigM} + \alpha\bigM)\abs{\pvi{}} &\text{if} \ \bigM \leq \sqrt{\reg/\alpha} \\
\icvx(\abs{\pvi{}} \leq \bigM) + 2\reg B(\abs{\pvi{}}\sqrt{\alpha/\reg}) &\text{if} \ \bigM > \sqrt{\reg/\alpha}
\end{cases}
\end{align}

when $\idxentry \in \setnone$ in (<ref>).

Case $\bigM \leq \sqrt{\reg/\alpha}$.
We have

sup_ ∈  - () -     = max{
sup_ ≤α+ /         ,
sup_ ≥α+ /  - ^24α +  α2α - ^2 +       }

= max{
sup_ ≤α+ /         ,
sup_ ≥α+ /  - ^24α +  α(2α - )^2 +       }

= max{
(α+ /)        ,
sup_ ≥α+ /  -  + α^2 +       }

where equality (<ref>) holds since $\bigM \leq \sqrt{\reg/\alpha} \implies \alpha\bigM + \reg\bigM \geq 2\alpha\bigM$.
Further, we remark that the supremum of the right member in (<ref>) lies in $\kset{\vvi{} \in \kR \cup \{+\infty\}}{\vvi{} \geq 0}$ when $\pvi{} \geq 0$ and in $\kset{\vvi{} \in \kR \cup \{-\infty\}}{\vvi{} \leq 0}$ when $\pvi{} \leq 0$.
Therefore, we obtain

\begin{align}
\label{eq:proof:relaxregfunc:9:4}
\textstyle\sup_{\abs{\vvi{}} \geq \alpha\bigM + \reg/\bigM} \pvi{}\vvi{} - \bigM\abs{\vvi{}} + \alpha\bigM^2 + \reg &= \textstyle\sup_{\abs{\vvi{}} \geq \alpha\bigM + \reg/\bigM} \abs{\vvi{}}(\abs{\pvi{}} - \bigM) + \alpha\bigM^2 + \reg \\
\begin{cases}
+\infty &\text{if} \ \abs{\pvi{}} > \bigM \\
(\alpha\bigM + \reg/\bigM)\abs{\pvi{}} &\text{if} \ \abs{\pvi{}} \leq \bigM
\end{cases} \\
\label{eq:proof:relaxregfunc:9:5}
&= (\alpha\bigM + \reg/\bigM)\abs{\pvi{}} + \icvx(\abs{\pvi{}} \leq \bigM)
\end{align}

Overall, combining (<ref>)-(<ref>) with (<ref>)-(<ref>) gives

\begin{align}
\biconj{(\node{\separable{\regfunc}{\idxentry}})}(\pvi{}) &= \max\big\{
(\alpha\bigM + \reg/\bigM)\abs{\pvi{}}, \
(\alpha\bigM + \reg/\bigM)\abs{\pvi{}} + \icvx(\abs{\pvi{}} \leq \bigM)
\big\} \\
\label{eq:proof:relaxregfunc:9:6}
&= (\alpha\bigM + \reg/\bigM)\abs{\pvi{}} + \icvx(\abs{\pvi{}} \leq \bigM)
\end{align}

for the case $\bigM \leq \sqrt{\reg/\alpha}$.

Case $\bigM > \sqrt{\reg/\alpha}.$
We have

sup_ ∈  - () -     = max{
sup_ ≤2√(α)         ,
sup_ ≥2√(α)  - ^24α + α2α - ^2 +       }

= max{
2√(α)        ,
sup_ ≥2√(α)  - ^24α + α2α - ^2 +       }

Further, we note that since $\bigM > \sqrt{\reg/\alpha} \implies 2\alpha\bigM > 2\sqrt{\reg\alpha}$, the right member in equation (<ref>) splits into

sup_ ≥2√(α)  - ^24α + α2α - ^2 +

=    max{
sup_2√(α) ≤ ≤2α  - ^24α + α2α - ^2 +        ,
sup_ ≥2α  - ^24α + α2α - ^2 +       }

=    max{
sup_2√(α) ≤ ≤2α  - ^24α +        ,
sup_ ≥2α  - ^24α + α(2α - )^2 +       }

=    max{
sup_2√(α) ≤ ≤2α  - ^24α +        ,
sup_ ≥2α  -  + α^2 +       }

On the one hand, the left member in equation (<ref>) can be expressed in closed-form as

\begin{align}
\textstyle\sup_{2\sqrt{\reg\alpha} \leq \abs{\vvi{}} \leq 2\alpha\bigM} \pvi{}\vvi{} - \tfrac{\vvi{}^2}{4\alpha} + \reg
\reg
-\textstyle\inf_{2\sqrt{\reg\alpha} \leq \abs{\vvi{}} \leq 2\alpha\bigM} \tfrac{\vvi{}^2}{4\alpha} - \pvi{}\vvi{}
\label{eq:proof:relaxregfunc:square abs + interval:buf biconjugate}
\\
&= \begin{cases}
2\sqrt{\reg\alpha} \abs{\pvi{}} &\text{if} \ \abs{\pvi{}} \leq \sqrt{\reg/\alpha} \\
\alpha\pvi{}^2 + \reg &\text{if} \ \abs{\pvi{}} \in \ ]\sqrt{\reg/\alpha}, \bigM] \\
2\alpha\bigM\abs{\pvi{}} - \alpha\bigM^2 + \reg &\text{if} \ \abs{\pvi{}} > \bigM
\end{cases}
\label{eq:proof:relaxregfunc:12:1}
\end{align}

by noting that the scalar defined as the orthogonal projection of \(2\alpha\pvi{}\) onto the interval \([2\sqrt{\reg\alpha},2\alpha\bigM]\) satisfies the necessary and sufficient first order optimality condition (see e.g. Corollary 3.68 in [Beck, 2017]) associated to the convex minimization problem involved in the right-hand side of (<ref>).
On the other hand, we remark that the supremum of the right member in equation (<ref>) lies in $\kset{\vvi{} \in \kR \cup \{+\infty\}}{\vvi{} \geq 0}$ when $\pvi{} \geq 0$ and in $\kset{\vvi{} \in \kR \cup \{-\infty\}}{\vvi{} \leq 0}$ when $\pvi{} \leq 0$, which yields

\begin{align}
\label{eq:proof:relaxregfunc:12:2}
\textstyle\sup_{\abs{\vvi{}} \geq 2\alpha\bigM} \pvi{}\vvi{} - \bigM\abs{\vvi{}} + \alpha\bigM^2 + \reg &= \textstyle\sup_{\abs{\vvi{}} \geq 2\alpha\bigM} \abs{\vvi{}}(\abs{\pvi{}} - \bigM) + \alpha\bigM^2 + \reg \\
\begin{cases}
+\infty &\text{if} \ \abs{\pvi{}} > \bigM \\
2\alpha\bigM\abs{\pvi{}} - \alpha\bigM^2 + \reg &\text{if} \ \abs{\pvi{}} \leq \bigM
\end{cases} \\
\label{eq:proof:relaxregfunc:12:3}
&= 2\alpha\bigM\abs{\pvi{}} - \alpha\bigM^2 + \reg + \icvx(\abs{\pvi{}} \leq \bigM)
\end{align}

By plugging (<ref>) and (<ref>)-(<ref>) into each member of (<ref>)-(<ref>), we deduce that

sup_ ≥2√(α)  - ^24α + α2α - ^2 +

,                 2α - α^2 + }    if   ≤√(/α)

α^2 +              ,                 2α - α^2 + }    if   ∈ ]√(/α),]

2α - α^2 +  ,                 +∞                 }    if   >

2√(α)    if   ≤√(/α)

α^2 +    if   ∈ ]√(/α),]

+∞   if   >

where equality (<ref>) holds since $2\sqrt{\reg\alpha} \leq 2\alpha\bigM$ and $-\alpha\bigM^2 + \reg \leq 0$, reminding that we consider the case $\bigM > \sqrt{\reg/\alpha}$.
Overall, combining (<ref>)-(<ref>) with (<ref>)-(<ref>) gives

\begin{align}
\biconj{(\node{\separable{\regfunc}{\idxentry}})}(\pvi{}) &=
\begin{cases}
2\sqrt{\reg\alpha}\abs{\pvi{}} &\text{if} \ \abs{\pvi{}} \leq \sqrt{\reg/\alpha} \\
\alpha\pvi{}^2 + \reg &\text{if} \ \abs{\pvi{}} \in \ ]\sqrt{\reg/\alpha},\bigM] \\
+\infty &\text{if} \ \abs{\pvi{}} > \bigM \\
\end{cases} \\
\label{eq:proof:relaxregfunc:14}
&= 2\reg B(\abs{\pvi{}}\sqrt{\alpha/\reg}) + \icvx(\abs{\pvi{}} \leq \bigM)
\end{align}

for the case $\bigM > \sqrt{\reg/\alpha}$.

§.§ Hyperparameters Calibration

To calibrate the value of $\reg$ in (<ref>) and the hyperparameters of the function $\pertfunc(\cdot)$, we use the  package [Dedieu et al., 2021] that can approximately solve some specific instances of the problem.
More specifically, we call the  procedure to perform a grid search over the values of $\reg$ and the hyperparameters of the function $\pertfunc(\cdot)$.
For each point in the grid, an approximate solution $\hat{\pv}$ to problem (<ref>) is constructed and an associated cross-validation score is computed.
We select the value of $\reg$ and the hyperparameters in $\pertfunc(\cdot)$ leading to the best cross-validation score.
For the synthetic instances considered in <Ref>, we only consider the candidates $\hat{\pv}$ in the grid with the best F1-score for the support recovery of the ground truth $\groundtruth$.