We present a novel screening methodology to safely discard irrelevant nodes within a generic bnb algorithm solving the \(\ell_0\)-penalized least-squares problem.
Our contribution is a set of two simple tests to detect sets of feasible vectors that cannot yield optimal solutions.
This allows to prune nodes of the bnb exploration tree, thus reducing the overall solution time.
One cornerstone of our contribution is a nesting property between tests at different nodes that allows to implement screening at low computational cost.
Our work leverages the concept of safe screening, well known for sparsity-inducing convex problems, and some recent advances in this field for \(\ell_0\)-penalized regression problems.

Sparse approximation, Mixed-integer problems, Branch and bound, Safe screening.

§ INTRODUCTION

Finding a sparse representation is a fundamental problem in the field of statistics, machine learning and inverse problems.
It consists in decomposing some input vector \(\obs\in\kR^\ddim\) as a linear combination of a few columns (dubbed atoms) of a dictionary \(\dic\in\kR^{\ddim\times\pdim}\).
This task can be addressed by solving

\begin{equation}
\label{prob:prob-base}
\min_{\pv \in \kR^{\pdim}} \ \tfrac{1}{2}\norm{\obs-\dic\pv}{2}^2 + \reg
\norm{\pv}{0}
\end{equation}

where \(\norm{\pv}{0}\) counts the number of nonzero entries in \(\pv\) and \(\reg>0\) is a tuning parameter. Unfortunately, problem (<ref>) has proven to be NP-hard in the general case <cit.>. This has led researchers to develop sub-optimal procedures to approximate its solution and address large-scale problems. Among the most popular, one can mention greedy algorithms and methodologies based on convex relaxation, see e.g., <cit.>.

Sub-optimal procedures are unfortunately only guaranteed to solve (<ref>) under restrictive conditions that are rarely met when dealing with highly-correlated dictionaries. On the other hand, there has been recently a surge of interest for methods solving (<ref>) exactly, see [3, 4, 5, 6, 7] to name a few. Many approaches leverage the fact that finding

\begin{equation}
\label{prob:prob} \tag{\(\mpb\)}
\begin{split}
\optobj \ =
\min_{\pv \in \kR^{\pdim}}
&\tfrac{1}{2} \norm{\obs-\dic\pv}{2}^2 + \reg \norm{\pv}{0}\\
&\text{s.t.}\,\,\, \kvvbar{\pv}_{\infty} \leq \bigM
\end{split}
\end{equation}

is equivalent to solve (<ref>), provided that \(\bigM\) is chosen large enough.
Interestingly, (<ref>) can be reformulated as a mip by introducing binary variables encoding the nullity of the entries of \(\pv\), see e.g., [4]. Besides, (<ref>) has been shown to be solvable for moderate-size problems by both commercial solvers [6] and tailored bnb algorithms [7].

In a recent paper, Atamtürk and Gómez extended the notion of safe screening introduced by El Ghaoui et al. in [8] from sparsity-promoting convex problems to non-convex \(\ell_0\)-penalized problems. In particular, they introduced a new methodology that allows to detect (some of) the positions of zero and non-zero entries in the minimizers of a particular \(\ell_0\)-penalized problem [9].
Their methodology is used as a preprocessing step of any algorithmic procedure and allows
to reduce the problem dimensionality.

In this paper, we make one step forward in the development of numerical methods addressing large-scale \(\ell_0\)-penalized problems by proposing  rules that can be implemented within any bnb algorithm.
In contrast to [9], we emphasize the existence of a nesting property between screening tests at different nodes.
This enables to (potentially) prune multiple nodes at any step of the optimization process with marginal cost.

Our exposition is organized as follows.
<Ref> gathers the main notations used in the paper.
In <Ref>, we describe a bnb algorithm tailored to (<ref>).
<Ref> presents our new tests and explains how to implement them efficiently within the bnb process.
In <Ref>, we assess the performance of our method on synthetic data. All the proofs are postponed to Appendix <ref>.

§ NOTATIONS

We use the following notational conventions throughout the paper.
Boldface uppercase (e.g., \(\dic\)) and lowercase (e.g., \(\pv\)) letters respectively represent matrices and vectors.
\(\0\) denotes the all-zeros vector.
Since the dimension will usually be clear from the context, it is omitted in the notation.
The \(\idxentrynode\)th column of a matrix \(\dic\) is denoted \(\atom_\idxentrynode\).
Similarly, the \(\idxentrynode\)th entry of a vector \(\pv\) is denoted \(\pve_{\idxentrynode}\).
Calligraphic letters ( \(\setnodesymb\)) are used to denote sets and the notation \(\card{\cdot}\) refers to their cardinality.
If \(\setnodesymb\subseteq\{1,\dotsc,\pdim\}\) and \(\pv\in\kR^\pdim\), \(\pv_{\setnodesymb}\) denotes to the restriction of \(\pv\) to its elements indexed by \(\setnodesymb\).
Similarly, \(\dic_{\setnodesymb}\) corresponds to the restriction of \(\dic\in\kR^{\ddim\times\pdim}\) to its columns indexed by \(\setnodesymb\).

§ BRANCH-AND-BOUND PROCEDURES

Branch-and-bound procedures refer to an algorithmic solution to address mip, among others [10].
It consists in implicitly enumerating all feasible solutions and applying pruning rules to discard irrelevant candidates.
When particularized to problem (<ref>), it can be interpreted as a search tree where a new decision regarding the nullity of an entry of the variable \(\pv\) is taken at each node as illustrated in <Ref>.
Formally, we define a node as a triplet \(\nodeSymb=(\setzero,\setone,\setnone)\) where :
i) \(\setzero\) and \(\setone\) contain the indices of the entries of $\pv$ which are forced to be zero and non-zero, respectively;
ii) \(\setnone\) gathers all the entries of $\pv$ for which no decision has been taken at this stage of the decision tree.

The bnb algorithm reads as follows: starting from node \(\nodeSymb=(\emptyset,\emptyset,\{1,\dotsc,\pdim\})\), the method alternates between processing the current node (bounding step) and selecting a new node (branching step).
The algorithm identifies the global minimum in a finite number of steps with a worst-case complexity equal to an exhaustive search.
In practice, its efficiency depends on both the ability to process nodes quickly and the number of nodes processed.
We review below specific choices for these two steps tailored to problem (<ref>).

§.§ Bounding step

When processing node \(\nodeSymb=(\setzero,\setone,\setnone)\), we prospect if any $\pv$ with zeros on $\setzero$ and nonzero values on $\setone$ can yield a global minimizer of (<ref>).
To that end, we look for the smallest objective value achieved by feasible candidates with respect to these constraints, i.e., the value of

\begin{equation} \label{prob:prob-node} \tag{\(\node{\spb}\)}
\begin{split}
\node{\sobj} \ = \
\min_{\pv \in \kR^{\pdim}}\ & \tfrac{1}{2} \norm{\obs-\dic\pv}{2}^2 + \reg \norm{\subnone{\pv}}{0} + \reg \card{\setone} \\
&\mathrm{s.t. } \ \norm{\pv}{\infty} \leq \bigM \ \text{ and }\ \subzero{\pv} = \0.
\end{split}
\end{equation}

One verifies that \(\node{\sobj} \geq \optobj\) since all minimizers of (<ref>) are also feasible for (<ref>).
Moreover, equality holds whenever the constraints given by \(\setzero\) and \(\setone\) match the configuration of one of the minimizers of (<ref>).

If \(\node{\sobj} > \optobj\), \(\nodeSymb\) can be pruned from the tree since no vector $\pv$ verifying
the constraints defined at this node (and therefore at any sub-nodes) can yield an optimal solution.
Unfortunately, the latter pruning rule is of poor practical interest since \(\optobj\) is not available.
Moreover evaluating \(\node{\sobj}\) when \(\setnone\neq\emptyset\) is also  a combinatorial problem.
To circumvent these problems, a relaxed version of the rule is devised: let \(\node{\robj}\) and \(\UB{\pobj}\) be respectively lower and upper bounds on \(\node{\sobj}\) and \(\optobj\). Then, a sufficient condition to prune node $\nodeSymb$ reads \(\node{\robj}> \UB{\pobj}\).

To obtain \(\node{\robj}\), a standard approach is to consider the following relaxation of (<ref>):

\begin{equation} \label{prob:relax-node} \tag{\(\node{\rpb}\)}
\begin{split}
\node{\robj} \ = \
\min_{\pv \in \kR^{\pdim}} \ & \tfrac{1}{2} \norm{\obs-\dic\pv}{2}^2 + \tfrac{\reg}{\bigM} \norm{\subnone{\pv}}{1} + \reg \card{\setone} \\
&\mathrm{s.t. } \ \norm{\pv}{\infty} \leq \bigM \ \text{ and } \ \subzero{\pv} = \0
\end{split}
\end{equation}

where the \(\ell_0\)-penalty term has been replaced by an \(\ell_1\)-norm.
Interestingly, (<ref>) is a constrained LASSO problem [11] and can be solved (to machine precision) in polynomial time by using one of the many methods suitable for this class of problems [12, 13].

The upper bound $\UB{\pobj}$ can be obtained by keeping track of the best known objective value of (<ref>) during the bnb process.
More precisely, we construct a feasible solution of (<ref>) at each node to obtain \(\node{\hobj}\).
The best known upper bound is then updated as \(\UB{\pobj}\gets\min(\UB{\pobj},\node{\hobj})\).
During the bnb procedure, \(\UB{\pobj}\) converges toward \(\optobj\) and relaxation are strengthened as new variables are fixed, allowing to prune nodes more efficiently.

§.§ Branching step

If node \(\nodeSymb=(\setzero,\setone,\setnone)\) has not been pruned during the bounding step, the tree exploration goes on.
An index \(\idxentrynode \in \setnone\) is selected according to some branching rule and two direct sub-nodes are created below \(\nodeSymb\) by imposing either “\(\pve_{\idxentrynode} = 0\)” or “\(\pve_{\idxentrynode} \neq 0\)”.
Finally, when all nodes have been explored or pruned, the bnb algorithm stops and any candidate yielding the best upper bound \(\UB{\sobj}\) is a minimizer of (<ref>).

First nodes of the bnb tree.
The root node \(\nodeSymbIter{0}\) corresponds to problem (<ref>) and each sub-node corresponds to a sub-problem (<ref>) with different fixed variables.
Standard implementation of a bnb where all nodes are processed.
Impact of applying tests on the bnb search tree:
at node \(\nodeSymbIter{0}\), test (<ref>) is passed for entries \(\idxentrynode_1\) and \(\idxentrynode_2\);
one can thus directly switch to the sub-node including constraints “\(\pve_{\idxentrynode_1}=0\)” and “\(\pve_{\idxentrynode_2}=0\)”, namely \(\nodeSymbIter{3}\).
The shaded nodes need not be explored.

§ TESTS

In this section, we present our new tests.
They aim at identifying, with marginal cost, nodes of the search tree that cannot yield a global minimum of (<ref>).

§.§ Dual properties

The crux of our procedure is a connection between the Fenchel dual problem of (<ref>) at two consecutive nodes.
Prior to expose this result, let us define for all \(\dv \in \kR^{\ddim}\) and \(\idxentrynode \in \{1,\dots,\pdim\}\) three families of pivot values as

\begin{equation}
\begin{array}{rcl}
\pivot{\idxentrynode}{}(\dv) &\triangleq&
\bigM(|\ktranspose{\atom}_\idxentrynode\dv| - \tfrac{\reg}{\bigM}) \\
\pivot{\idxentrynode}{0}(\dv) &\triangleq&
\bigM\relu{|\ktranspose{\atom}_\idxentrynode\dv| - \tfrac{\reg}{\bigM}} \\
\pivot{\idxentrynode}{1}(\dv) &\triangleq&
\bigM\relu{\tfrac{\reg}{\bigM}-|\ktranspose{\atom}_\idxentrynode\dv|}
\end{array}
\end{equation}

where \(\relu{z} \triangleq \max(0, z)\) for all scalars \(z\). We now express the Fenchel dual of (<ref>) with respect to these quantities.

The Fenchel dual of (<ref>) is given by

\begin{equation}
\label{prob:dual-node}
\begin{split}
\node{\dobj} \ = \ \max_{\dv \in \kR^{\ddim}} \node{\dfunc}(\dv) \triangleq& \tfrac{1}{2} \norm{\obs}{2}^2 - \tfrac{1}{2} \norm{\obs-\dv}{2}^2 \\
&- \textstyle\sum_{\idxentrynode \in \setnone} \pivot{\idxentrynode}{0}(\dv) - \textstyle\sum_{\idxentrynode \in \setone} \pivot{\idxentrynode}{}(\dv)
\end{split}
\tag{$\node{\dpb}_l$}
\end{equation}

and strong duality holds for (<ref>)-(<ref>), i.e., \(\node{\robj} = \node{\dobj}\).
Moreover, if \((\opt{\pvrelax},\dvopt)\) is a couple of primal-dual solutions, one has

\begin{equation} \label{eq:optimality condition}
\dvopt = \obs - \dic\opt{\pvrelax}
\end{equation}

Hence, at a given node \(\nodeSymb\), the dual objective of (<ref>) is the sum of a term common to all nodes and some well-chosen pivot values.
Interestingly, the dual objective functions between two consecutive nodes only differ from one pivot value as shown in the following result:

Let \(\nodeSymb=(\setzero,\setone,\setnone)\) and \(\idxscreen \in \setnone\), then \(\forall\dv \in \kR^{\ddim}\),

\begin{align}
\dfunc^{\nodePlusZero{\nodeSymb}{\idxscreen}}(\dv) &= \node{\dfunc}(\dv) + \pivot{\idxscreen}{0}(\dv)
\label{eq:dual-link:nodeplus0}\\
\dfunc^{\nodePlusOne{\nodeSymb}{\idxscreen}}(\dv) &= \node{\dfunc}(\dv) + \pivot{\idxscreen}{1}(\dv)
\label{eq:dual-link:nodeplus1}
\end{align}

where \(\nodePlusZero{\nodeSymb}{\idxscreen}\) and \(\nodePlusOne{\nodeSymb}{\idxscreen}\) denote the two direct sub-nodes of \(\nodeSymb\) where index \(\idxscreen\) has been swapped from \(\setnone\) to \(\setzero\) or \(\setone\).

§.§ tests

We now expose our proposed screening strategy to identify nodes of the tree that provably cannot yield a global minimizer of (<ref>).
Let \(\nodeSymb\) be a node and \(\UB{\pobj}\) an upper bound on \(\opt{\pobj}\).
We have by strong duality between (<ref>)-(<ref>) that

\begin{equation} \label{eq:obj-ordering}
\forall\dv\in\kR^\ddim,
\quad
\node{\dfunc}(\dv) \leq \node{\dobj} = \node{\robj} \leq \node{\sobj}
\end{equation}

Hence, combining (<ref>) with <Ref> leads to the next result:

Let \(\idxscreen \in \setnone\) be some index.
Then, \(\forall\dv \in \kR^{\ddim}\),

() + 0()    >              ^ >

() + 1()    >            ^ >

Stated otherwise, <Ref> describes a simple procedure to identify some sub-nodes of \(\nodeSymb\) that cannot yield an optimal solution of (<ref>). In particular, if both (<ref>) and (<ref>) pass for a given index, no feasible vector with respect to the constraints defined at \(\nodeSymb\) is a minimizer of (<ref>) (\(\nodeSymb\) can therefore be pruned).

Our proposed procedure differs from the pruning methodology described in <Ref> in several aspects.
Performing a test only requires the evaluation of one single inner product.
They can therefore be implemented at marginal cost compared to the overall cost of the bounding step.
Very importantly, they also inherit from the following nesting property:

Let \(\nodeSymb=(\setzero,\setone,\setnone)\) be a node. Then, if test (<ref>) or (<ref>) passes for index \(\idxscreen\in\setnone\), then it also passes for any sub-node \(\nodeSymb'=(\setzero',\setone',\setnone')\subset\nodeSymb\) such that \(\idxscreen\in\setnone'\).

In particular, <Ref> has the following consequence.
Assume that two distinct indexes \(\idxscreen\) and \(\idxscreen'\) pass test (<ref>) at node \(\nodeSymb\).
Then index \(\idxscreen'\) will also pass (<ref>) at node \(\nodePlusOne{\nodeSymb}{\idxscreen}\).
The same consequence holds for test (<ref>) and node \(\nodePlusZero{\nodeSymb}{\idxscreen}\).
In other words, if several tests pass at a given node, one can simultaneously prune several sub-nodes of \(\nodeSymb\), as illustrated in <Ref>. This is in contrast with the bounding procedure described in <Ref> which has to process each sub-node individually before applying any potential pruning operation.

§.§ Implementation considerations

Implementing tests described by <Ref> requires the knowledge of an upper bound \(\UB{\sobj}\) of \(\optobj\) and a suitable \(\dv\in\kR^\ddim\), as close as possible to the maximizer of (<ref>).
The value of \(\UB{\sobj}\) can be obtained at no additional cost since the standard implementation of bnb already requires storing such valid upper bounds on \(\opt{\pobj}\).
To obtain a relevant candidate for \(\dv\), we suggest the following procedure: assuming the method used to solve (<ref>) generates a sequence of iterates \(\{\pv^{(t)}\}_{t\in\kN}\) that converges to a global minimizer, we set

\begin{equation}
\label{eq:dual scaling-ish}
\kforall t\in\kN,\quad
\dv^{(t)} = \obs - \dic\pv^{(t)}
\end{equation}

This choice enjoys two desirable properties.
First, the sequence \(\{\dv^{(t)}\}_{t\in\kN}\) converges toward a maximizer of (<ref>) as a consequence of the optimality condition (<ref>).
Second, both \(\dv^{(t)}\) and \(\ktranspose{\dic}\dv^{(t)}\) are already evaluated by most solvers as they correspond to the residual error and the (negative) gradient of the least-squares term of the objective function, respectively.
Thus, the computational cost of evaluating the screening tests at all undecided entries is marginal compared to the cost needed to address (<ref>).
The latter choice for \(\dv^{(t)}\) suggests performing tests during the bounding step.
Hence, when some entries passes the test at node $\nodeSymb$, the bnb algorithm immediately switches to some sub-node, without solving (<ref>) at high precision.

§ NUMERICAL RESULTS

We finally report simulation results illustrating the relevance of the methodology on two different setups.

§.§ Experimental setups

For each trial, we generate a new realization of \(\dic \in \kR^{\ddim\times\pdim}\) and \(\obs\in\kR^{\ddim}\) as follows.
In the Gaussian setup, we set \((\ddim,\pdim)=(500,1000)\) and the entries of \(\dic\) are i.i.d. realizations of a normal distribution.
In the Tœplitz setup, we set \((\ddim,\pdim)=(500,300)\) and the first column of \(\dic\) contains a sampled sinc function.
The other columns are obtained by shifting the entries in a row-wise fashion so that \(\dic\) inherits from a Tœplitz structure.
In both setups, the columns of $\dic$ are normalized to one.
To obtain \(\obs\), we first sample a \(\sparsitylevel\)-sparse vector $\sparsevector\in\kR^\pdim$ with uniformly distributed non-zero entries.
Each non-zero entry is set to \(s(1 + \kvbar{a})\), where \(s\in\{-1,+1\}\) is a random sign and \(a\) is a realization of the normal distribution.
In our experiment, we chose \(\sparsitylevel \in \{5,7,9\}\).
Finally, the observation is constructed as \(\obs = \dic\sparsevector + \boldsymbol{\epsilon}\)
where the entries of \(\boldsymbol{\epsilon}\) are i.i.d. realizations of a centered Gaussian distribution with standard deviation \(\sigma = \norm{\dic\sparsevector}{2} / \sqrt{10\ddim}\).
Such a design leads to a SNR of \(10\)dB in average [14].
We tune \(\reg\) statistically as in [15], i.e., by setting \(\reg = 2\sigma^2\log(\pdim/\sparsitylevel - 1)\).
Finally, we empirically set \(\bigM = 1.5\norm{\ktranspose{\dic}\obs}{\infty}\) as advised in [7].

We compare three methods that address (<ref>):
i) , that uses CPLEX [6], a state-of-the-art commercial mip solver;
ii) , the algorithm presented in [7], which is (up to our knowledge) the fastest bnb algorithm tailored for (<ref>);
iii)  which corresponds to  enhanced with the methodology presented in <Ref>.
Note that both  and  leverage an efficient implementation of the Active-Set algorithm <cit.> to solve (<ref>) and the branching rule described in <cit.>.

Experiments are run on a MacOS with an i7 CPU, clocked at 2.2 GHz and with 8Go of RAM.
We restrict calculations on a single core to avoid bias due to parallelization capabilities.
The Academic Version 20.1 of CPLEX is used and both  and  are implemented in Julia v1.5 [17].
Results are averaged over 100 instances of problem (<ref>).

§.§ Method comparison

<Ref> compares the average number of nodes treated (the number of relaxations (<ref>) solved at machine precision) and the solution time for the three methods and all simulation setups.
One observes that  outperforms the two other methods on all scenarii and all figures of merit.
We also note that, compared to , the reduction in the solution time is more significant than the reduction in the number of processed nodes.
A thorough examination of our results indicates that the bounding step is performed all the faster as many variables are fixed to zero in (<ref>). Our node-screening methodology allows to reach quickly nodes where the bounding step is performed with a lower computational cost.

As a final remark, we mention that  relies on an efficient C++ implementation while  and  are implemented in Julia.
As C++ usually runs faster than Julia, there is still room for improvements in the comparison of  and  with .

§ CONCLUSION

In this paper, we presented a novel methodology aiming at accelerating the resolution of the \(\ell_0\)-penalized least-squares problem with a bnb solver. Our contribution leverages a nesting property between the node-screening tests at two consecutive nodes. Our method leads to significant improvements in terms of number of nodes explored and resolution time on two simulated datasets.

§ PROOFS

This appendix gathers the proofs of the results derived in <Ref>. We refer to [7] for proofs regarding <Ref>.

§.§ Proof of <Ref>

Our proof of <Ref> leverages the Fenchel conjugate of a function \(f\) and denoted \(\conj{f}\). In particular, we first state the following technical lemma whose proof is postponed to the end of the section

Define for all \(\z\in\kR^\ddim\) and \(\pv\in\kR^\pdim\)

\begin{align}
f_1(\z) &\triangleq \tfrac{1}{2}\norm{\obs-\z}{2}^2 + \varepsilon
\label{eq:lemma:conjugate f_1} \\
f_2(\pv) &\triangleq
\tfrac{\reg}{\bigM}\norm{\pv_{\setnodesymb}}{1}
+ \Icvx_{\kset{\pv'}{\pv'_{\setnodesymb'}=\0}}(\pv)
+ \Icvx_{\mathcal{B}_\infty(M)}(\pv)
\label{eq:lemma:conjugate f_2}
\end{align}

Then, for all \(\dv\in\kR^\ddim\) and \(\fenchdualvar\in\kR^\pdim\)

\begin{align}
\conj{f_1}(\dv) &=
\tfrac{1}{2}(\norm{\obs+\dv}{2}^2 - \norm{\obs}{2}^2)
-  \varepsilon \\
\conj{f_2}(\fenchdualvar) &=
\sum_{\idxentrynode\in\setnodesymb}
\bigM\relu{\kvbar{\fenchdualvarel_\idxentrynode}-\tfrac{\reg}{\bigM}}
\bigM\norm{\fenchdualvar_{\setnodesymb''}}{1}
\end{align}

In the latter result, \(\varepsilon\) is a scalar, \((\setnodesymb,\setnodesymb', \setnodesymb'')\) is a partition of \(\{1,\dotsc,\pdim\}\), \(\mathcal{B}_\infty(M)\) denotes the \(\ell_\infty\)-ball of \(\kR^\pdim\) centered at \(\0\) with radius \(\bigM\), \(\Icvx_{\mathcal{A}}\) denotes the indicator function of set \(\mathcal{A}\subseteq\kR^\pdim\) defined for all \(\pv\in\kR^\pdim\) by \(\Icvx_{\mathcal{A}}(\pv)=0\) if \(\pv\in\mathcal{A}\) and \(+\infty\) otherwise.

We now use <Ref> to prove : a) the formulation of the dual problem (<ref>), b) the strong-duality relation linking (<ref>)-(<ref>) and c) the relation (<ref>). In the sequel, we let \(\nodeSymb=(\setzero,\setone,\setnone)\) be a given node.

a) Dual problem (<ref>).
Using the functions defined in <Ref> with \(\varepsilon=\reg\card{\setone}\), \(\setnodesymb=\setnone\), \(\setnodesymb'=\setzero\) and \(\setnodesymb'' = \setone\), the relaxed problem (<ref>) can be rewritten as

\begin{equation}
\label{eq:proof:rewritting relaxed}
\min_{\pv\in\kR^{\pdim}} \
f_1(\dic\pv) + f_2(\pv)
\end{equation}

Then, the Fenchel dual of (<ref>) is given by (see <cit.>[Note that, unlike <cit.>, we use the formulation of the dual problem that involves a maximization problem.])

\begin{equation}
\max_{\dv\in\kR^\ddim} \
- \conj{f_1}(-\dv)
- \conj{f_2}(\ktranspose{\dic}\dv)
\end{equation}

One concludes the proof by remarking that for all \(\dv\in\kR^\ddim\),

\begin{align}
\conj{f_2}(\ktranspose{\dic}\dv) - \varepsilon
\;=\;&
\sum_{\idxentrynode\in\setnone} \pivot{\idxentrynode}{0}(\dv)
+ \bigM\kvvbar{\ktranspose{\dic}_{\setone}\dv}_1 - \tfrac{\bigM}{\bigM} \lambda\card{\setone}
\nonumber \\
\;=\;&
\sum_{\idxentrynode\in\setnone} \pivot{\idxentrynode}{0}(\dv)
\sum_{\idxentrynode\in\setone} \pivot{\idxentrynode}{}(\dv)
\nonumber.
\end{align}

b) Strong duality.
Strong duality holds as a consequence of <cit.>.
More particularly, one easily verifies that the condition in item vii) is fulfilled since the functions \(f_1\) and \(f_2\) are proper and continuous.

c) Relation (<ref>).
Let \((\node{\pvrelax},\node{\dv})\) be a couple of primal-dual solution of (<ref>)-(<ref>).
Since strong duality holds, we have by item ii) of <cit.> that[Again, we note that (<ref>) involves the opposite of \(\node{\dv}\) since the authors consider the reformulation of the dual as a minimization problem.]

\begin{equation}
\label{eq:proof:dual link opt condition}
-\node{\dv}\in\partial f_1(\dic\node{\pvrelax})
\end{equation}

where \(\partial f_1(\dic\node{\pvrelax})\) denotes the sub-differential of \(f_1\) evaluated at \(\dic\node{\pvrelax}\).
One finally obtains (<ref>) by noticing that \(f_1\) is differentiable at \(\dic\node{\pvrelax}\) so that the sub-differential reduces to the singleton \(\{\nabla f_1(\dic\node{\pvrelax})\}\).

§.§ Proof of <Ref>.

The Fenchel conjugate of \(f_1\) evaluated at \(\dv\in\kR^{\ddim}\) is defined as

\begin{equation}
\label{eq:proof:fenchel dual f_1}
\conj{f_1}(\dv) = \max_{\z\in\kR^{\ddim}}\
\ktranspose{\dv}\z - f_1(\z)
\end{equation}

One easily sees that (<ref>) is an unconstrained concave optimization problem whose objective function is differentiable.
We then obtain (<ref>) by noticing that the maximum is attained at
\(\opt{\z}=\obs + \dv\).

Let \(\fenchdualvar\in\kR^\pdim\).
By definition of \(\conj{f}_2\) and using the fact that \(f_2(\pv)=+\infty\) for all \(\pv\in\kR^\pdim\) such that \(\pv_{\setnodesymb'}\neq\0\), We have

\begin{align}
\conj{f}_2(\fenchdualvar)
\;=\;&
\max_{ \kset{\pv\in\kR^\pdim}{\kvvbar{\pv}_{\infty} \leq \bigM} }
\ktranspose{\fenchdualvar}\pv
- \tfrac{\lambda}{\bigM} \kvvbar{\pv_{\setnodesymb}}_1
\nonumber \\
\;=\;&
\sum_{\idxentrynode\in\setnodesymb}
\max_{ \kset{\pve\in\kR}{\kvbar{\pve} \leq \bigM} }
\fenchdualvarel_{\idxentrynode}\pve_{\idxentrynode}
- \tfrac{\lambda}{\bigM} \kvbar{\pve_{\idxentrynode}}
\nonumber \\
&\quad +\max_{ \kset{\pv_{\setnodesymb''}\in\kR^{\card{\setnodesymb''}}}{\kvvbar{\pv_{\setnodesymb''}}_{\infty} \leq \bigM} }
\ktranspose{\fenchdualvar}_{\setnodesymb''}\pv_{\setnodesymb''}
\nonumber
\end{align}

where the optimization problem can be split since it is linear and \(\setnodesymb,\setnodesymb''\) have non empty intersection by definition.
Let \(\idxentrynode\in\setnodesymb\) and consider the problem

\begin{equation}
\label{eq:proof:subproblem 1}
\opt{\pve_{\idxentrynode}} \in \kargmax_{ \kset{\pve\in\kR}{\kvbar{\pve} \leq \bigM} } \
\fenchdualvarel_{\idxentrynode}\pve_{\idxentrynode}
- \tfrac{\lambda}{\bigM} \kvbar{\pve_{\idxentrynode}}
\end{equation}

We let the reader check that \(\opt{\pve}_\idxentrynode\) defined by

\begin{equation}
\opt{\pve}_\idxentrynode =
\begin{cases}
0 &\text{ if } \kvbar{\fenchdualvarel_{\idxentrynode}} - \tfrac{\lambda}{\bigM} \leq 0 \\
\sign{\fenchdualvarel_{\idxentrynode}} \bigM &\text{ otherwise}
\end{cases}
\end{equation}

is the maximizer of (<ref>). One finally expresses the maximum as a sum of pivot values by injecting the value of \(\opt{\pve}_\idxentrynode\) in the objective function.

Second, see that

\begin{align}
&\max_{ \kset{\pv_{\setnodesymb''}\in\kR^{\card{\setnodesymb''}}}{\kvvbar{\pv_{\setnodesymb''}}_{\infty} \leq \bigM} }
\ktranspose{\fenchdualvar}_{\setnodesymb''}\pv_{\setnodesymb''}
\nonumber \\
&\qquad \;=\;
\bigM\times \max_{ \kset{\pv_{\setnodesymb''}\in\kR^{\card{\setnodesymb''}}}{\kvvbar{\pv_{\setnodesymb''}}_{\infty} \leq 1} }
\ktranspose{\fenchdualvar}_{\setnodesymb''}\pv_{\setnodesymb''}
\nonumber \\
&\qquad \;=\;
\bigM\kvvbar{\fenchdualvar_{\setnodesymb''}}_1
\nonumber
\end{align}

where the last equality holds since one recognizes the Fenchel dual of the indicator function of the unit-ball with respect to the \(\ell_\infty\)-norm <cit.> in the penultimate line.

§.§ Proof of <Ref>

Let \(\nodeSymb=(\setzero,\setone,\setnone)\) be a node and \(\idxentrynode\) be some element of \(\setnone\) such that \(\nodePlusZero{\nodeSymb}{\idxentrynode}\) defines a child node.
Using (<ref>), we have for all \(\dv\in\kR^\ddim\)

\begin{align}
\dfunc^{\nodePlusZero{\nodeSymb}{\idxentrynode}}(\dv)
\dfunc^{\nodeSymb}(\dv)
\;=\;&
\sum_{\idxentrynode' \in \setnone} \pivot{\idxentrynode'}{0}(\dv)
\sum_{\idxentrynode' \in \setnone\setminus\{\idxentrynode\}} \pivot{\idxentrynode'}{0}(\dv)
\nonumber \\
\;=\;&
\pivot{\idxentrynode}{0}(\dv).
\nonumber
\end{align}

Similar calculations lead for all \(\dv\in\kR^\ddim\) to

\begin{align}
\dfunc^{\nodePlusOne{\nodeSymb}{\idxentrynode}}(\dv)
\dfunc^{\nodeSymb}(\dv)
\;=\;&
\pivot{\idxentrynode}{0}(\dv)
\pivot{\idxentrynode}{}(\dv)
\nonumber \\
\;=\,&
\pivot{\idxentrynode}{1}(\dv)
\nonumber
\end{align}

where the last line uses the relation \(\relu{\dve} - \relu{-\dve} = \dve\) that hold for all scalar \(\dve\).

§.§ Proof of <Ref>

Let \(\nodeSymb=(\setzero,\setone,\setnone)\) be a node.
We first state the following lemma which can easily be proved by induction using (<ref>) and (<ref>).

Let \(\nodeSymb'=(\setzero',\setone',\setnone')\) be a child node of \(\nodeSymb\).
Then, for all \(\dv\in\kR^\ddim\)

\begin{equation}
\dfunc^{\nodeSymb'}(\dv)
\dfunc^{\nodeSymb}(\dv)
+ \sum_{\idxentrynode\in\setzero'\setminus\setzero} \pivot{\idxentrynode}{0}(\dv)
+ \sum_{\idxentrynode\in\setone'\setminus\setzero} \pivot{\idxentrynode}{1}(\dv)
\end{equation}

We now prove <Ref>.
Let \(\idxscreen\in\setnone\) be some index such that test (<ref>) or (<ref>) passes at node \(\nodeSymb\).
To easier our exposition, we only prove the result for test (<ref>) as the same rationale can be used for test (<ref>).

Hence, assume that test (<ref>) passes for index \(\idxscreen\) at node \(\nodeSymb\), that

\begin{equation}
\node{\dfunc}(\dv) + \pivot{\idxscreen}{1}(\dv) > \UB{\pobj}
\end{equation}

Let \(\nodeSymb'=(\setzero',\setone',\setnone')\) be a child node of \(\nodeSymb\) such that \(\idxscreen\notin\setnone'\).
Using <Ref>, we have

\begin{equation}
\dfunc^{\nodeSymb'}(\dv)
\dfunc^{\nodeSymb}(\dv)
+ \sum_{\idxentrynode\in\setzero'\setminus\setzero} \pivot{\idxentrynode}{0}(\dv)
+ \sum_{\idxentrynode\in\setone'\setminus\setzero} \pivot{\idxentrynode}{1}(\dv)
\end{equation}

Since the pivots values \(\{\pivot{\idxentrynode}{0}(\dv)\}_{\idxentrynode=1}^\pdim\)
\(\{\pivot{\idxentrynode}{1}(\dv)\}_{\idxentrynode=1}^\pdim\) are all nonnegative by construction, we necessarily have \(\dfunc^{\nodeSymb'}(\dv) \geq \dfunc^{\nodeSymb}(\dv)\).
Hence the test (<ref>) also passes for index \(\idxscreen\) at node \(\nodeSymb'\).

[1]
Xiaojun Chen, Dongdong Ge, Zizhuo Wang, and Yinyu Ye,
“Complexity of unconstrained $\ell_2-\ell_p$ minimization,”
Mathematical Programming, vol. 143, no. 1-2, pp. 371–383,
November 2014.

[2]
Simon Foucart and Holger Rauhut,
A Mathematical Introduction to Compressive Sensing,
Springer New York, 2013.

[3]
Ryuhei Miyashiro and Yuichi Takano,
“Subset selection by Mallows' ${C}_p$: A mixed integer programming
Expert Systems with Applications, vol. 42, no. 1, pp. 325–331,
Jan. 2015.

[4]
Sébastien Bourguignon, Jordan Ninin, Hervé Carfantan, and Marcel
“Exact sparse approximation problems via mixed-integer programming:
Formulations and computational performance,”
IEEE Transactions on Signal Processing, vol. 64, no. 6, pp.
1405–1419, 2015.

[5]
Dimitris Bertsimas, Angela King, and Rahul Mazumder,
“Best subset selection via a modern optimization lens,”
The Annals of Statistics, vol. 44, no. 2, Apr. 2016.

[6]
CPLEX User's Manual,
“Ibm ilog cplex optimization studio,”
Version, vol. 12, pp. 1987–2018, 1987.

[7]
Ramzi Ben Mhenni, Sébastien Bourguignon, Marcel Mongeau, Jordan Ninin, and
Hervé Carfantan,
“Sparse branch and bound for exact optimization of \(\ell_0\)-norm
penalized least squares,”
in ICASSP. IEEE, 2020, pp. 5735–5739.

[8]
Laurent El Ghaoui, Vivian Viallon, and Tarek Rabbani,
“Safe feature elimination for the lasso and sparse supervised
learning problems,”

[9]
Alper Atamturk and Andrés Gómez,
“Safe screening rules for l0-regression from perspective
in International conference on machine learning. PMLR, 2020,
pp. 421–430.

[10]
Eugene L Lawler and David E Wood,
“Branch-and-bound methods: A survey,”
Operations research, vol. 14, no. 4, pp. 699–719, 1966.

[11]
Brian R Gaines, Juhyun Kim, and Hua Zhou,
“Algorithms for fitting the constrained lasso,”
Journal of Computational and Graphical Statistics, vol. 27, no.
4, pp. 861–871, 2018.

[12]
Neal Parikh and Stephen Boyd,
“Proximal algorithms,”
Foundations and Trends in optimization, vol. 1, no. 3, pp.
127–239, 2014.

[13]
Honglak Lee, Alexis Battle, Rajat Raina, and Andrew Y Ng,
“Efficient sparse coding algorithms,”
in Advances in neural information processing systems, 2007, pp.

[14]
Don H Johnson,
“Signal-to-noise ratio,”
Scholarpedia, vol. 1, no. 12, pp. 2088, 2006.

[15]
Charles Soussen, Jérôme Idier, David Brie, and Junbo Duan,
“From Bernoulli-Gaussian deconvolution to sparse signal
IEEE Transactions on Signal Processing, vol. 59, no. 10, pp.
4572–4584, 2011.

[16]
Stephen Wright, Jorge Nocedal, et al.,
“Numerical optimization,”
Springer Science, vol. 35, no. 67-68, pp. 7, 1999.

[17]
Jeff Bezanson, Stefan Karpinski, Viral B Shah, and Alan Edelman,
“Julia: A fast dynamic language for technical computing,”

[18]
Heinz H. Bauschke and Patrick L. Combettes,
Convex Analysis and Monotone Operator Theory in Hilbert Spaces,
Springer International Publishing, 2017.