





\section{Introduction}

This paper focuses on the resolution of the so-called ``\(\ell_0\)-regularized least-squares problem'' given by
\begin{align}
	\stepcounter{equation}
	\tag{\theequation-\(\mathcal{P}\)}
	\label{prob:prob-base} 
	\min_{\mathbf{\pvi} \in \kR^{n}} \ \mathrm{P}(\mathbf{\pvi}) \triangleq \tfrac{1}{2}\kvvbar{\mathbf{y} - \mathbf{A}\mathbf{\pvi}}_2^2 + \lambda \norm{\mathbf{\pvi}}{0}
\end{align}
where \(\mathbf{y} \in \kR^{m}\) and \(\mathbf{A} \in \kR^{m\timesn}\) are input data, \(\lambda>0\) is a regularization parameter and \(\kvvbar{\cdot}_0\) denotes the $\ell_0$-pseudonorm which counts the number of non-zero elements in its argument.  

Solving~\eqref{prob:prob-base} is of paramount interest in many scientific fields such as statistics, machine learning or inverse problems~\cite{candes2005decoding,donoho2005stable,tropp2010computational}. 
Unfortunately, this problem also turns out to be NP-hard~\cite[Th.~1]{Chen2012}. 
Hence, the last decades have seen a flurry of contributions proposing tractable procedures able to recover approximate solutions of~\eqref{prob:prob-base}. 
Canonical examples include greedy algorithms or methodologies based on convex relaxations, see~\cite[Sec.~3.2 and Ch.~4]{foucart2013invitation}. 
Although these procedures successfully recover the actual solutions of \eqref{prob:prob-base} in ``easy'' setups, they usually fall short 
for more challenging instances of the problem.  
This observation, combined with some recent advances in integer optimization and hardware performance,
has revived the interest in methods solving~\eqref{prob:prob-base} exactly, see~\textit{e.g.},~\cite{cplex2009v12,bourguignon2015exact,Bertsimas2016,bertsimas2021unified,hazimeh2021sparse,guyard2022node}.

In this paper, we propose a new strategy, dubbed \textit{``safe peeling''}, to accelerate the exact resolution of~\eqref{prob:prob-base} via a \gls{bnb} procedure. 
In a nutshell, our contribution is a computationally simple test to identify some intervals of \(\kR^n\) which cannot contain a solution of~\eqref{prob:prob-base}. 
This information is then embedded in the \gls{bnb} nodes, leading to tighter convex relaxations and more aggressive pruning of the nodes of the decision tree. 
Our numerical experiments show that the proposed method leads to a significant reduction of the solving time as compared to state-of-the-art concurrent methods. 
The name \textit{``safe peeling''} comes from the fact that the proposed method enables to reduce (or in more figurative terms, ``to peel'') the feasible set of the problem at each node of the decision tree while safely preserving the correctness of the \gls{bnb} procedure.


The rest of the paper is organized as follows.
\Cref{sec:bnb} describes the main ingredients of \gls{bnb} methods.
Our peeling strategy is presented in \Cref{sec:peeling} and its performance is illustrated in \Cref{sec:results}.
\ifthenelse{\boolean{techreport}}{
	Proofs of the results presented in the following are postponed to the appendix.
}{
	All the proofs are postponed to the technical report accompanying this paper~\cite{peeling2022technicalreport}.
} 

\section{Proof of the main results}


\subsection{Proof of \texorpdfstring{Lemma~\ref{lemma: consequence properties peeled interval}}{Lemma 1}}  \label{proof:lemma:consequence properties peeled interval}

We first have 
	\begin{align*}
		\node{p}(\bigL',\bigU')\geq \node{p}(\mathbf{\bigLi},\mathbf{\bigUi})
	\end{align*}
since $[\bigL',\bigU']\subseteq [\mathbf{\bigLi},\mathbf{\bigUi}]$ from \eqref{eq:peeling: constraint mnew <= m} and the left-hand side corresponds to the minimum value of a problem more constrained than the right-hand side.
The reverse implication in \eqref{eq:peeling: equivalence} is thus always verified. 

The direct implication results from the following observations.	
If
$\node{p}(\mathbf{\bigLi},\mathbf{\bigUi})\leq \UB{p}$, then we have 
	\begin{align*}
		\node{p}(\mathbf{\bigLi},\mathbf{\bigUi}) = \node{p}(\bigL',\bigU')
	\end{align*}
since \eqref{eq:peeling: constraint} ensures that the  minimizers of $\node{\mathrm{P}}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi})$ do not belong to $[\mathbf{\bigLi},\mathbf{\bigUi}]\setminus[\bigL',\bigU']$ and are therefore also minimizers of $\node{\mathrm{P}}(\mathbf{\pvi};\bigL',\bigU')$. 
We thus obtain the direct implication by contraposition.\\ 			



\subsection{Proof of~\texorpdfstring{\eqref{eq:lower bound on peeling metric}}{(16)}} \label{proof:lower bounds}

We first note that
\begin{subequations}
	\begin{align}
		\node{p}_\alpha(\mathbf{\bigLi},\mathbf{\bigUi}) &= \inf_{\mathbf{\pvi}\in\kR^n} \node{\mathrm{P}}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi}) + \eta(x_{\idxpeeled} > \alpha) \\
		&\geq \inf_{\mathbf{\pvi}\in\kR^n} \node{\mathrm{R}}_{\idxpeeled}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi}) + \eta(x_{\idxpeeled} > \alpha)\\
		&\geq \inf_{\mathbf{\pvi}\in\kR^n} \node{\mathrm{R}}_{\idxpeeled}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi}) + \eta(x_{\idxpeeled} \geq \alpha) \\
		&= \min_{\mathbf{\pvi}\in\kR^n} \node{\mathrm{R}}_{\idxpeeled}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi}) + \eta(x_{\idxpeeled} \geq \alpha)\label{eq:relaxed-perturbed-problem b}
	\end{align}
\end{subequations}
where 
\begin{multline} 
	\node{\mathrm{R}}_{\idxpeeled}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi})
	\triangleq 
	\tfrac{1}{2}\kvvbar{\mathbf{y} - \mathbf{A}\mathbf{\pvi}}_2^2 
	+ \lambda \sum_{\idxentry \in \setidx_\symbnone\backslash \idxpeeled} \Big(\tfrac{\pospart{x_{\idxentry}}}{u_\idxentry}
	- \tfrac{\pospart{-x_{\idxentry}}}{l_\idxentry}\Big)

	\\
	+ \lambda (\card{\setidx_\symbone}+1)
	+ \eta(\subzero{\mathbf{\pvi}} = \mathbf{0})
	+ \eta(\mathbf{\pvi} \in [\mathbf{\bigLi},\mathbf{\bigUi}]) 

\end{multline}
and the last equality follows from the fact that the objective function is lower-semi-continuous and its domain is closed and bounded. 

The lower bound in \eqref{eq:lower bound on peeling metric} then stems from the Fenchel-Rockafellar dual problem of \eqref{eq:relaxed-perturbed-problem b}.
More specifically,  we first notice that \eqref{eq:relaxed-perturbed-problem b} can be rewritten as
\begin{equation}\label{eq:primal formulation}
	\min_{\mathbf{\pvi}\in\kR^n} f(\mathbf{A}\mathbf{\pvi}) 
	+ g(\mathbf{\pvi})
\end{equation}
with
\begin{subequations}
	\begin{align}
		\label{eq:def-f-peel}
		f(\mathbf{z}) \,=\,& \tfrac{1}{2}\kvvbar{\mathbf{y} - \mathbf{z}}_2^2
		\\
		%
		\label{eq:def-g-peel}
		g(\mathbf{\pvi}) \,=\,& g_\idxpeeled(x_\idxpeeled) + \sum_{\idxentry\neq \idxpeeled} 
		g_\idxentry(x_\idxentry)
	\end{align}
\end{subequations}
and
\begin{align*}

	g_\idxpeeled(x)
	&= \eta(x \in [\alpha,u_{\idxpeeled}]) + \lambda\\

	g_\idxentry(x)
	&=
	\begin{cases}
		\eta(x=0) 	&\text{ if } \idxentry\in\setidx_\symbzero \\
		\eta(x \in [l_{\idxentry},u_{\idxentry}]) + \lambda &\text{ if } \idxentry\in\setidx_\symbone \\
		\eta(x \in [l_{\idxentry},u_{\idxentry}]) + \lambda(\tfrac{\pospart{x}}{u_\idxentry} - \tfrac{\pospart{-x}}{l_\idxentry}) & \text{ if } \idxentry\in\setidx_\symbnone\backslash \idxpeeled.
	\end{cases}
\end{align*}
Second, we have from standard results of convex optimization~\cite[Ch.~12]{Beck2017aa} that the Fenchel-Rockafellar dual problem of~\eqref{eq:primal formulation} reads
\begin{equation}
	\label{eq:dual formulation d'un certain probleme qui n'a pas de nom}
	\max_{\mathbf{\dvi}\in\kR^m} \ -\conj{f}(-\mathbf{\dvi}) - \conj{g}(\transp{\mathbf{A}}\mathbf{\dvi}),
\end{equation}
where $\conj{f}$ and $\conj{g}$ denote the convex conjugates of \(f\) and \(g\) (see~\cite[Def.~4.1]{Beck2017aa}).
In particular, we have \(\forall \mathbf{\pvi}\in\kR^n, \mathbf{\dvi}\in\kR^m\): 
\begin{align}\label{eq:duality}
	f(\mathbf{A}\mathbf{\pvi}) + g(\mathbf{\pvi})
	\geq
	-\conj{f}(-\mathbf{\dvi}) - \conj{g}(\transp{\mathbf{A}}\mathbf{\dvi})
\end{align}
as a consequence of weak duality applied to the pair~\eqref{eq:primal formulation}-\eqref{eq:dual formulation d'un certain probleme qui n'a pas de nom}. 
Our lower bound in \eqref{eq:lower bound on peeling metric}  corresponds to a particularization of the right-hand side of \eqref{eq:duality} to \eqref{eq:def-f-peel}-\eqref{eq:def-g-peel}. 
In particular, applying the definition of the convex conjugate to the function $f$, one easily obtains:
\begin{equation}\label{eq:conj f}
	\conj{f}(\mathbf{\dvi}) = -\tfrac{1}{2}\norm{\mathbf{y}}{2}^2 + \tfrac{1}{2}\norm{\mathbf{y} + \mathbf{\dvi}}{2}^2
	.
\end{equation}
Invoking \Cref{lemma:conj-setone} in \Cref{app:proof} with \((b,l,u)=(\lambda,\alpha,u_\idxpeeled)\), we obtain
\begin{equation}
	\label{eq:conj g2}
	\conj{g}_\idxpeeled(v)
	=
	u_{\idxpeeled}\pospart{v } - \alpha \pospart{-v } - \lambda.
\end{equation}
Applying again \Cref{lemma:conj-setone} with
\begin{equation*}
	(b,l,u) = 
	\begin{cases}
		 (0,0,0) &\text{if} \ \idxentry \in \setidx_\symbzero \\
		(\lambda,l_{\idxentry},u_{\idxentry}) &\text{if} \ \idxentry \in \setidx_\symbone,
	\end{cases}
\end{equation*} 
leads to 
\begin{equation}
	\label{eq:conj g1}
	\conj{g}_\idxentry(v)
	=
	\begin{cases}
		0 	&\text{ if } \idxentry\in\setidx_\symbzero \\
		\pospart{u_{\idxentry}v } + \pospart{l_{\idxentry}v} - \lambda &\text{ if } \idxentry\in\setidx_\symbone 
	\end{cases}
\end{equation}
where we have used the fact that for all \(\idxentry \in \setidx_\symbone\), \(u_{\idxentry}\pospart{v}=\pospart{u_{\idxentry}v}\) (resp. \(l_{\idxentry}\pospart{-v}=-\pospart{l_{\idxentry}v}\)) since \(u_{\idxentry}\geq0\) (resp. \(l_{\idxentry} \leq0\)). 	
Similarly, using \Cref{lemma:conj-setnone} in \Cref{app:proof} with $(a,l,u) = (\lambda,l_{\idxentry},u_{\idxentry}) $ 
we obtain: 
\begin{equation}
	\label{eq:conj g3}
	\conj{g}_\idxentry(v)
	=
	\pospart{u_{\idxentry}v  - \lambda } + \pospart{l_{\idxentry}v  -\lambda} \quad \text{ if } \idxentry\in\setidx_\symbnone\backslash \idxpeeled.
\end{equation}	
Finally, combining \eqref{eq:conj f}-\eqref{eq:conj g3} leads to the desired result. \\

\subsection{Proof of \texorpdfstring{Lemma~\ref{eq:lemma parallele application}}{Lemma~3}} \label{sec:proof lemma parallele application}

Let \([\mathbf{\bigLi}',\mathbf{\bigUi}']\) and \([\mathbf{\bigLi}'',\mathbf{\bigUi}'']\) be two intervals satisfying~\eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m}, and defines \(\mathcal{S}=[\mathbf{\bigLi}',\mathbf{\bigUi}']\cap[\mathbf{\bigLi}'',\mathbf{\bigUi}'']\).
We first note that \([\mathbf{\bigLi}',\mathbf{\bigUi}']\cap[\mathbf{\bigLi}'',\mathbf{\bigUi}'']\) can be rewritten as
\begin{equation*}
	\mathcal{S}
	=
	\bigtimes_{\idxentry=1}^n \kintervcc{
		\max(l'_\idxentry, l''_\idxentry)
	}{
		\min(u'_\idxentry, u''_\idxentry)
	}
\end{equation*}
where \(\bigtimes\) denotes the Cartesian product of sets, and therefore defines an interval.
It thus remains to show that \(\mathcal{S}\) fulfills~\eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m}.

First, the inclusion \(\mathcal{S}\subset[\mathbf{\bigLi},\mathbf{\bigUi}]\) follows from the fact that both \([\mathbf{\bigLi}',\mathbf{\bigUi}']\) and \([\mathbf{\bigLi}'',\mathbf{\bigUi}'']\) verify \eqref{eq:peeling: constraint mnew <= m}. 
Second, the fact that \(\mathcal{S}\) fulfills~\eqref{eq:peeling: constraint} is a consequence of the following set equality: 
\begin{equation*}
	[\mathbf{\bigLi},\mathbf{\bigUi}]\setminus\mathcal{S} = \kparen{[\mathbf{\bigLi},\mathbf{\bigUi}]\setminus[\mathbf{\bigLi}',\mathbf{\bigUi}']}\cup \kparen{[\mathbf{\bigLi},\mathbf{\bigUi}]\setminus[\mathbf{\bigLi}'',\mathbf{\bigUi}'']}
	.
\end{equation*}
Hence, if \(\mathbf{\pvi}\in[\mathbf{\bigLi},\mathbf{\bigUi}]\setminus\mathcal{S}\) then \(\mathbf{\pvi}\in\kparen{[\mathbf{\bigLi},\mathbf{\bigUi}]\setminus[\mathbf{\bigLi}',\mathbf{\bigUi}']}\) or \(\mathbf{\pvi}\in\kparen{[\mathbf{\bigLi},\mathbf{\bigUi}]\setminus[\mathbf{\bigLi}'',\mathbf{\bigUi}'']}\).
Assume without loss of generality that \(\mathbf{\pvi}\in\kparen{[\mathbf{\bigLi},\mathbf{\bigUi}]\setminus[\mathbf{\bigLi}',\mathbf{\bigUi}']}\).
Since \([\mathbf{\bigLi}',\mathbf{\bigUi}']\) fulfills~\eqref{eq:peeling: constraint}, we then have \(\node{\mathrm{P}}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi})>\UB{p}\).
One concludes the proof by noting that the latter rationale holds for all \(\mathbf{\pvi}\in[\mathbf{\bigLi},\mathbf{\bigUi}]\setminus\mathcal{S}\).


\subsection{Proof of {Lemma}~\texorpdfstring{\ref{corollary:heritage}}{4}} \label{sec:proof corollary:heritage}

If the peeling procedure defined in \Cref{prop:entry-wise peeling} is applied at node $\nu'$ with $[\bigL',\bigU']$  as input to generate a new interval $[\bigL'',\bigU'']$, we have by construction:
	\begin{subequations}
	\begin{align}
		\forall \mathbf{\pvi}\in[\bigL',\bigU']\setminus [\bigL'',\bigU'']&:\,
		\mathrm{P}^{\nu'}(\mathbf{\pvi};\bigL',\bigU')>\UB{p}
		\label{eq:peeling: constraint 2}
		\\
		[\bigL'',\bigU'']&\subseteq [\bigL',\bigU']
		.
		\label{eq:peeling: constraint mnew <= m 2}
	\end{align}
	\end{subequations}
We note that \eqref{eq:peeling: constraint 2} is equivalent to 
 	\begin{align}\label{eq:peeling: constraint 2b}
		\forall \mathbf{\pvi}\in[\bigL',\bigU']\setminus [\bigL'',\bigU'']:\,
		\mathrm{P}^{\nu'}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi})>\UB{p}
	\end{align}
since $\mathbf{\pvi}\in [\bigL',\bigU']$ and $[\bigL',\bigU']\subseteq[\mathbf{\bigLi},\mathbf{\bigUi}]$ by hypothesis. 

Moreover, we have 
	\begin{align*}
		\mathrm{P}^{\nu'}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi}) \geq \node{\mathrm{P}}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi})
		.
	\end{align*}
since $\mathcal{X}^{\nu'} \subset \node{\mathcal{X}}$ for any child node $\nu'$ of $\nu$. 
Hence, any interval $[\bigL',\bigU']$ verifying \eqref{eq:peeling: constraint} at node $\nu$ obviously also fulfills \eqref{eq:peeling: constraint} at $\nu'$, that is 
	\begin{align}\label{eq:peeling: constraint ho ho}
		\forall \mathbf{\pvi}\in[\mathbf{\bigLi},\mathbf{\bigUi}]\setminus [\bigL',\bigU']:\,
		\mathrm{P}^{\nu'}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi})>\UB{p}
		.
	\end{align}
Combining \eqref{eq:peeling: constraint mnew <= m 2}, \eqref{eq:peeling: constraint 2b} and \eqref{eq:peeling: constraint ho ho}, we finally obtain the result. 

\section{Technical lemmas} 
\label{app:proof}


In this appendix, we derive the expressions of two conjugate functions appearing in the derivation of \eqref{eq:lower bound on peeling metric}. The results are encapsulated in \Cref{lemma:conj-setone,lemma:conj-setnone} below. 

\begin{lemma}
	\label{lemma:conj-setone}
	Let \(\kfuncdef{g}{\kR}{\kR}\) be defined for all \(x\in\kR\) as
	\begin{equation*}
		g(x) = 
		b + \eta(x\in[l,u])		
	\end{equation*}
	for $b\in\kR$ and \(l \leq u\). 
	Then, 
	\begin{equation*}
		g^\star(v) =
		u \pospart{v} - l \pospart{-v} -b
		.
	\end{equation*}	
\end{lemma}

\begin{proof}
	By definition of the convex conjugate (see~\cite[Def.~4.1]{Beck2017aa}), one has \(\forall v\in\kR\):
	\begin{subequations}
	\begin{align*}
		\conj{g}(v) 
		&= \sup_{x\in\kR} \; \cvix - g(x)\\
		&= \max_{l\leq x \leq u} \; \cvix - b\\
		&= u \pospart{v} - l \pospart{-v} -b
	\end{align*}
	\end{subequations}
\end{proof}

{

\newcommand{\funcsymbAppBLemmaTwo}{g} 
\begin{lemma}
	\label{lemma:conj-setnone}
	Let \(\kfuncdef{\funcsymbAppBLemmaTwo}{\kR}{\kR}\) be defined for all \(x\in\kR\) as
	\begin{equation*}
		\funcsymbAppBLemmaTwo(x) = 
		\eta(x\in[l,u])		
		+\tfrac{a}{u} \pospart{x} 
		-\tfrac{a}{l}\pospart{-x} 
	
	\end{equation*}
	for $a \geq 0$

	 and \(l < 0 < u\). 
	Then, 
	\begin{equation*}
		\conj{\funcsymbAppBLemmaTwo}(v) 
		= 
		\pospart{\bigUiv-a} 
		+ \pospart{\bigLiv - a} 
	
		.
	\end{equation*}	
	This expression remains valid for \(l \leq 0 \leq u\) as long as we use the convention ``$0/0=0$''.
\end{lemma}
\begin{proof}

	\newcommand{\beta} %\overline{\cvi}{\beta}
	\newcommand{\mu} %\upsilon %\overline{\bigUi}{\mu}

		Prior to proving the lemma, let us first show that the following result holds true: for all scalars \(\beta} %\overline{\cvi\in\kR,\mu} %\upsilon %\overline{\bigUi\geq 0\)
		, we have
		\begin{equation}
			\label{eq:proof lemma conj setnone:intermediary result}
			\max_{0 \leq x \leq \mu} %\upsilon %\overline{\bigUi} \; \bufcvix - \tfrac{a}{\mu} %\upsilon %\overline{\bigUi}x = \pospart{\mu} %\upsilon %\overline{\bigUi\beta} %\overline{\cvi - a}.
		\end{equation}
		Note that the maximum in~\eqref{eq:proof lemma conj setnone:intermediary result} is attained since the optimization problem amounts to maximizing a linear function over a (nonempty) compact set.
		Now, if \(\mu} %\upsilon %\overline{\bigUi=0\) we have
		\begin{equation*}
			\max_{0 \leq x \leq \mu} %\upsilon %\overline{\bigUi} \; \bufcvix - \tfrac{a}{\mu} %\upsilon %\overline{\bigUi}x 
			= 0
		\end{equation*} 
		by using the convention ``$0/0=0$''. 
		Moreover, if \(\mu} %\upsilon %\overline{\bigUi>0\):	
		\begin{subequations}
			\begin{align*}
				\max_{0 \leq x \leq \mu} %\upsilon %\overline{\bigUi} \; \bufcvix - \tfrac{a}{\mu} %\upsilon %\overline{\bigUi}x
				\,=\,& \max_{0 \leq x' \leq 1} \; x'(\mu} %\upsilon %\overline{\bigUi\beta} %\overline{\cvi - a)
				\\
				\,=\,& \pospart{\mu} %\upsilon %\overline{\bigUi\beta} %\overline{\cvi - a}.
			\end{align*}
		\end{subequations}	 
		Overall, one can summarize the two cases as
		\begin{equation}
			\max_{0 \leq x \leq \mu} %\upsilon %\overline{\bigUi} \; \bufcvix - \tfrac{a}{\mu} %\upsilon %\overline{\bigUi}x
			=
			\pospart{\mu} %\upsilon %\overline{\bigUi\beta} %\overline{\cvi - a}
		\end{equation}
		since \(a\geq 0\).

	
	\vspace*{.2cm}

	We now turn to the proof of \Cref{lemma:conj-setnone}. 
	By definition of the convex conjugate (see~\cite[Def.~4.1]{Beck2017aa}), one has \(\forall v\in\kR\): 
	\begin{equation*}
		\conj{\funcsymbAppBLemmaTwo}(v) = \sup_{x\in\kR} \; \cvix - \funcsymbAppBLemmaTwo(x).
	\end{equation*}
	The above problem can be equivalently expressed as
	\begin{equation}
		\label{eq:conj-split}
		\conj{\funcsymbAppBLemmaTwo}(v) = \max\big\{\conj{\funcsymbAppBLemmaTwo}_1(v) ; \conj{\funcsymbAppBLemmaTwo}_2(v)\big\}
	\end{equation}
	where 
	\begin{subequations}
		\begin{align*}
			\conj{\funcsymbAppBLemmaTwo}_1(v) &\triangleq \sup_{x \geq 0} \; \cvix - \funcsymbAppBLemmaTwo(x) \\
			\conj{\funcsymbAppBLemmaTwo}_2(v) &\triangleq \sup_{x \leq 0} \; \cvix - \funcsymbAppBLemmaTwo(x).
		\end{align*}
	\end{subequations}
	On the one hand, \(\conj{\funcsymbAppBLemmaTwo}_1(v)\) can be rewritten as 
	\begin{subequations}
		\begin{align*}
			\conj{\funcsymbAppBLemmaTwo}_1(v) 
			\,=\,& \sup_{0 \leq x \leq u} \; \cvix - \tfrac{a}{u}x 
			.
		\end{align*}
	\end{subequations}
	By applying~\eqref{eq:proof lemma conj setnone:intermediary result} with \((\beta} %\overline{\cvi,\mu} %\upsilon %\overline{\bigUi) = (v,u)\), we then obtain:
	\begin{equation} \label{eq:expression g1*}
		\conj{\funcsymbAppBLemmaTwo}_1(v) = \pospart{\bigUiv - a}
		.
	\end{equation}
	On the other hand, we have
	\begin{subequations}
		\begin{align*}
			\conj{\funcsymbAppBLemmaTwo}_2(v) 
			\,=\,& \sup_{l \leq x \leq 0} \; \cvix - \tfrac{a}{l}x
			\nonumber \\
			\,=\,& \sup_{0 \leq x' \leq -l} \; -\cvix' - \tfrac{a}{-l}x' 
			.
		\end{align*}
	\end{subequations}
	Hence, using ~\eqref{eq:proof lemma conj setnone:intermediary result} with \((\beta} %\overline{\cvi, \mu} %\upsilon %\overline{\bigUi)=(v,-l)\) leads to 
	\begin{equation} \label{eq:expression g2*}
		\conj{\funcsymbAppBLemmaTwo}_2(v)
		=
		\pospart{\bigLiv - a}
		.
	\end{equation}
	Finally, plugging \eqref{eq:expression g1*}-\eqref{eq:expression g2*} into \eqref{eq:conj-split} yields
	\begin{subequations}
		\begin{align*}
			\conj{\funcsymbAppBLemmaTwo}(v) 
			&= \max\big\{
				\pospart{\bigUiv - a
				;
				\pospart{\bigLiv - a
				\big\} \\
			&= \pospart{\bigUiv - a} + \pospart{\bigLiv - a}
		\end{align*}
	\end{subequations}
	where the last equality follows from the fact that\footnote{If $\bigUiv - a>0$ then necessarily $u>0$ since $a\geq 0$. We thus have $v >\tfrac{a}{u}\geq 0$.
	This implies that $\bigLiv - a\leq 0$ since $l\leq 0$ and $a\geq 0$.} 
	\begin{align*}
		\bigUiv - a>0 &\implies \bigLiv - a\leq 0
		.
	\end{align*}
\end{proof}

}

\section{Macros}

\begin{itemize}
    \item Data : $k$, $m$, $n$, $\mathbf{\pvi}^{\dagger}$, $f$, $g$, $\mathbf{y}$, $\mathbf{A}$, $\mathbf{a}$, $\lambda$, $M$
    \item Sub/super scripts : $\opt{x}$, $\UB{x}$, $\LB{x}$, $\subzero{x}$, $\subone{x}$, $\subnone{x}$, $\node{x}$, $\peelleft{x}{i}$, $\peelright{x}{i}$, $\transp{x}$
    \item Problems : $\mathcal{P}$, $\mathcal{R}$, $\mathcal{D}$
    \item Objective values : $p$, $d$, $r$ 
    \item Objective functions : $\mathrm{P}$, $\mathrm{D}$, $\mathrm{R}$ 
    \item Variables : $x$, $w$, $v$, $\mathbf{\pvi}$, $\mathbf{\dvi}$, $\mathbf{\cvi}$, $u$, $l$, $\mathbf{\bigUi}$, $\mathbf{\bigLi}$
    \item Entries : $\idxentry$, $\idxentrynode$, $\idxpeeled$
    \item BnB : $\nu$, $\nodeSymbIter{0}$, $\setidx_\symbzero$, $\setidx_\symbone$, $\setidx_\symbnone$
    \item Peeling : $\pivot{\idxpeeled}(\cdot)$, $l'$, $m'$, $\epsilon$
    \item Experiments : $\boldsymbol{\Sigma}$, $\rho$, $\tau$, $\mathbf{n}} %\boldsymbol{\epsilon$, $\gamma$, \texttt{BnB}\xspace, \texttt{BnB+S}\xspace, \texttt{BnB+P}\xspace
    \item Operators : $\abs{\cdot}$, $\mathrm{argmin}$, $\mathrm{argmax}$, $\conj{f}$, $\card{\cdot}$, $\mathbf{I}$, $\norm{\cdot}{\cdot}$, $\pospart{\cdot}$, $\negpart{\cdot}$, $\sign{\cdot}$, $\partial$, $\intervint{1}{n}$
\end{itemize}
\section{Conclusion}

In this paper, we presented a tractable strategy (named ``\emph{peeling}'') to tighten the box constraints used in a \gls{bnb} procedure tailored to $\ell_0$-penalized least-squares problems. Unlike the standard approach which imposes \textit{one global} constraint to the problem, our strategy aims to locally refine the box constraints \textit{at each node} of the decision tree. 
This refinement enables to strengthen the convex relaxations used in the pruning decisions made by the \gls{bnb} procedure and can lead to significant improvements in terms of solving time, as emphasized by our simulation results.

\section{Principles of BnB Methods}
\label{sec:bnb}
In this section, we recall the main principles of \gls{bnb} procedures. 
Due to space limitation, we only review the elements of interest to introduce the proposed peeling method. 
We refer the reader to \cite[Ch.~7]{Wolsey:1998op} for an in-depth treatment of the subject. 

\subsection{Pruning}

The crux of \gls{bnb} methods consists in identifying and discarding some subsets of $\kR^{n}$ which do not contain a minimizer of \eqref{prob:prob-base}.  
To do so, one constructs a decision tree in which each node corresponds to a particular subset of $\kR^{n}$. 
In our context, a tree node is identified by two disjoint subsets of $\intervint{1}{n}$, say $\setidx_\symbzero$ and $\setidx_\symbone$.  
The goal at node $\nu \triangleq (\setidx_\symbzero,\setidx_\symbone)$ is to detect whether a solution of \eqref{prob:prob-base} can be attained within
\begin{align}
	\node{\mathcal{X}}
	\triangleq \kset{\mathbf{\pvi}\in\kR^n}{\subzero{\mathbf{\pvi}} = \mathbf{0}, \ \subone{\mathbf{\pvi}} \neq \mathbf{0}}
	,
\end{align}
where \(\mathbf{\pvi}_{\nu_k}\) denotes the restriction of $\mathbf{\pvi}$ to its elements in \(\nu_k\). 
In particular, let $\opt{\mathcal{X}}$ be the (non-empty) set of minimizers of \eqref{prob:prob-base}.
Then, if some upper bound $\UB{p}$ on $\opt{p}$ is known and if we let
\begin{equation}\label{eq:node_problem}
	\node{p} \triangleq \inf_{\mathbf{\pvi} \in \kR^{n}} \node{\mathrm{P}}(\mathbf{\pvi})
\end{equation}
with $\node{\mathrm{P}}(\mathbf{\pvi}) \triangleq 
	 \mathrm{P}(\mathbf{\pvi}) + \eta(\mathbf{\pvi} \in \node{\mathcal{X}})$, 
	we obtain the following implication:
\begin{align}\label{eq:pruning condition}
	\node{p} > \UB{p} \implies  \node{\mathcal{X}}\cap\opt{\mathcal{X}}=\emptyset
	.
\end{align}
In words, if the left-hand side of \eqref{eq:pruning condition} is satisfied, $\node{\mathcal{X}}$ contains \textit{no} solution of \eqref{prob:prob-base} and can therefore be discarded from the search space of the optimization problem. 
This operation is usually referred to as ``\textit{pruning}''. 



\subsection{Bounding and relaxing}\label{sec:sec:Bounding and relaxing}

Making a pruning decision at node $\nu$ 
 requires the knowledge of $\UB{p}$ and $\node{p}$. 
On the one hand, finding $\UB{p}$ is an easy task since the value of the objective function in \eqref{prob:prob-base} at any feasible point constitutes an upper bound on $\opt{p}$.
On the other hand, evaluating $\node{p}$ is NP-hard. 
This issue can nevertheless be circumvented by finding a tractable lower bound $\node{r}$ on $\node{p}$ and relaxing \eqref{eq:pruning condition} as
\begin{align}\label{eq:relaxed pruning condition}
	\node{r} > \UB{p} \implies  \node{\mathcal{X}}\cap\opt{\mathcal{X}}=\emptyset
	.
\end{align}
One ubiquitous approach in the literature to find such a lower bound consists in:
\begin{itemize}
	\item[\textit{i)}] 
	Adding an extra term ``$\eta(\mathbf{\pvi} \in [\mathbf{\bigLi},\mathbf{\bigUi}])$'' to the cost function of \eqref{prob:prob-base}, for some \textit{well-chosen}
	bounds $\mathbf{\bigLi}\in\kR-^n$, $\mathbf{\bigUi}\in\kR+^n$.\footnote{This additional constraint usually takes the form ``$\mathbf{\pvi}\in [-\1M, \1M]$'' with $M>0$ and is known as \textit{``Big-M''} constraint in the literature, see~\cite[Sec.~3]{bourguignon2015exact}}.	
	In particular, the new constraint ``$\mathbf{\pvi} \in [\mathbf{\bigLi},\mathbf{\bigUi}]$'' must lead to a problem fully equivalent to \eqref{prob:prob-base}, that is
	\begin{align}\label{eq: condition equivalence}
		\forall \opt{\mathbf{\pvi}}\in\opt{\mathcal{X}}:\ \opt{\mathbf{\pvi}}\in [\mathbf{\bigLi},\mathbf{\bigUi}]. 
	\end{align}
	\item[\textit{ii)}] Considering the convex relaxation of the function $\kvvbar{\cdot}_0$ on $\node{\mathcal{X}}\cap [\mathbf{\bigLi},\mathbf{\bigUi}]$, that is (with the convention ``$0/0=0$'') 
	\begin{align} \label{eq:convex relaxation l0 pseudo norm}
		\|\mathbf{\pvi}\|_0
		\geq 
		|\setidx_\symbone| + \sum_{\idxentry\in \setidx_\symbnone} 
		\frac{\pospart{x_{\idxentry}}}{u_{\idxentry}} - \frac{\pospart{-x_{\idxentry}}}{l_{\idxentry}}
		,
	\end{align}
	where $\setidx_\symbnone \triangleq \intervint{1}{n} \setminus(\setidx_\symbzero\cup\setidx_\symbone)$. 
\end{itemize}
  
\noindent
On the one hand, item \textit{i)} implies that pruning test \eqref{eq:pruning condition} involves the following quantity (rather than $\node{p}$):  
\begin{align}
	\label{eq:node problem constr}
	\stepcounter{equation}
	\tag{\theequation-\(\node{\mathcal{P}}\)}
	\node{p}(\mathbf{\bigLi},\mathbf{\bigUi}) 
		= \inf_{\mathbf{\pvi} \in \kR^{n}} \ \node{\mathrm{P}}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi})
\end{align}
where $\node{\mathrm{P}}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi}) \triangleq \node{\mathrm{P}}(\mathbf{\pvi}) + \eta(\mathbf{\pvi} \in [\mathbf{\bigLi},\mathbf{\bigUi}])$.
On the other hand, 
a lower bound $\node{r}(\mathbf{\bigLi},\mathbf{\bigUi})$ on $\node{p}(\mathbf{\bigLi},\mathbf{\bigUi})$ can be obtained by using \eqref{eq:convex relaxation l0 pseudo norm}  and solving
\begin{align}
	\label{eq:equivalent node relaxed problem:big M}
	\stepcounter{equation}
	\tag{\theequation-\(\node{\mathcal{R}}\)}
	\node{r}(\mathbf{\bigLi},\mathbf{\bigUi})  = \min_{\mathbf{\pvi}\in\kR^n} \node{\mathrm{R}}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi})
\end{align}
where 
\begin{multline*} 
	\node{\mathrm{R}}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi}) 
	\triangleq 
	\tfrac{1}{2}\kvvbar{\mathbf{y} - \mathbf{A}\mathbf{\pvi}}_2^2 
	+ \lambda \sum_{\idxentry \in \setidx_\symbnone} \tfrac{\pospart{x_{\idxentry}}}{u_\idxentry}
	- \tfrac{\pospart{-x_{\idxentry}}}{l_\idxentry}

	\\
	+ \lambda \card{\setidx_\symbone} 
	+ \eta(\subzero{\mathbf{\pvi}} = \mathbf{0})
	+ \eta(\mathbf{\pvi} \in [\mathbf{\bigLi},\mathbf{\bigUi}]) 

	.
\end{multline*}
We note that \eqref{eq:equivalent node relaxed problem:big M} is a convex  problem and can be solved efficiently (to good accuracy) via numerous numerical procedures, see \emph{e.g.,}\xspace~\cite[Ch.~10]{Beck2017aa}.

In practice, the choice of $\mathbf{\bigLi}$ and $\mathbf{\bigUi}$ must respect two conflicting imperatives. 
First, the new constraint ``$\mathbf{\pvi} \in [\mathbf{\bigLi},\mathbf{\bigUi}]$'' should not modify the solution of our target problem~\eqref{prob:prob-base} and condition \eqref{eq: condition equivalence} must therefore be verified. 
Since $\opt{\mathcal{X}}$ is obviously not accessible beforehand, 
this suggests that the entries of $\mathbf{\bigLi}$ and $\mathbf{\bigUi}$ should be chosen  ``large-enough'' in absolute values.\footnote{Some heuristics are commonly used in the literature to select proper values of the bounds, see \cite[Sec.~V.B]{bourguignon2015exact}, \cite[Sec.~5.1]{guyard2022node} or \cite[Sec.~4]{mhenni2020sparse}.}
Second, the tightness of $\node{r}(\mathbf{\bigLi},\mathbf{\bigUi})$ with respect to $\node{p}(\mathbf{\bigLi},\mathbf{\bigUi})$ degrades with the spread of the set \([\mathbf{\bigLi},\mathbf{\bigUi}]\).\footnote{
	This impairment pertains to a large class of MIP and is pointed out since (at least) the 90s, see \textit{e.g.},~\cite{camm1990cutting}.
} 
In particular, the right-hand side of \eqref{eq:convex relaxation l0 pseudo norm} tends to zero  when $\mathbf{\bigLi} \ll \mathbf{\pvi}$ or $\mathbf{\pvi} \ll \mathbf{\bigUi}$.  
Therefore, setting the entries of $\mathbf{\bigLi}$ and $\mathbf{\bigUi}$ with too large absolute values is likely to degrade the effectiveness of the relaxed pruning decision~\eqref{eq:relaxed pruning condition}. 

In the next section, we propose a solution to address this problem by deriving a methodology which locally tightens the constraint $\mathbf{\pvi} \in [\mathbf{\bigLi},\mathbf{\bigUi}]$ at each node of the decision tree while preserving the correctness of the \gls{bnb} procedure. 

\section{Numerical results} 
\label{sec:results}

\subsection{Experimental setup}

We consider instances of problem~\eqref{prob:prob-base} with dimensions \((m,n)=(100,150)\).
For each trial, new realizations of \(\mathbf{A},\mathbf{y}\) and \(\lambda\) are generated as follows.
Each row of the dictionary \(\mathbf{A}\) is drawn from a multivariate normal distribution with zero mean and covariance matrix \({\boldsymbol{\Sigma}}\in\kR^{n\times n}\). The entry $(i,j)$ of $\boldsymbol{\Sigma}$ is defined as 
\ifthenelse{\boolean{techreport}}{\(\Sigma_{ij}=\rho^{\abs{i-j}}\) \(\forall i,j\in\intervint{1}{n}\) with $\rho=10^{-1}$} 
{\(\Sigma_{ij}=10^{-\abs{i-j}}\) \(\forall i,j\in\intervint{1}{n}\)}. 
Each realization of \(\mathbf{y}\) is generated in two steps.
We first create 
\ifthenelse{\boolean{techreport}}
{a $k$-sparse vector \(\mathbf{\pvi}^{\dagger} \in \kR^{n}\) with evenly-distributed non-zero components. We choose $k=5$.} 
{a \(5\)-sparse vector \(\mathbf{\pvi}^{\dagger} \in \kR^{n}\) with evenly-distributed non-zero components.} 
The non-zero entries are defined as \(x^{\dagger}_{\idxentry} = \sign{r_{\idxentry}} + r_{\idxentry}\) where  \(r_{\idxentry}\) is an independent realization of a zero-mean Gaussian with variance $\sigma^2$. 
We then set \(\mathbf{y} = \mathbf{A}\mathbf{\pvi}^{\dagger} + \mathbf{n}} %\boldsymbol{\epsilon\) for some  zero-mean white Gaussian noise \(\mathbf{n}} %\boldsymbol{\epsilon\). The variance of the noise is adjusted so that the $\mathrm{SNR} \triangleq 10 \log_{10}(\norm{\mathbf{A}\mathbf{\pvi}^{\dagger}}{2}^2/\norm{\mathbf{n}} %\boldsymbol{\epsilon}{2}^2)$
is equal to \(15\)dB.
The parameter \(\lambda\) is calibrated for each instance of~\eqref{prob:prob-base} using the cross-validation tools of the \textsc{L0Learn} package \cite{hazimeh2020fast} with the default parameters\ifthenelse{\boolean{techreport}}{. 
	More specifically, we call the
	\texttt{cv.fit} 
	procedure that takes \(\mathbf{y}\) and \(\mathbf{A}\) as inputs and returns couples 
	$(\mathbf{\pvi}_{\lambda_i},c_{\lambda_i})$ from a grid of values $\{\lambda_i\}_{i \in \kN}$ selected data-dependently by the package.
	The vector $\mathbf{\pvi}_{\lambda_i}$
	is an approximate solution of~\eqref{prob:prob-base} where
	the weight of the $\ell_0$-norm is set to $\lambda_i$ and \(c_{\lambda_i}\) is an associated cross-validation score 
	computed on 10 folds of $m/10$ randomly sampled rows in $\mathbf{A}$ and entries in $\mathbf{y}$.
	In our experiments, we then set
	\begin{equation}
		\lambda = \kargmin_{\lambda_i} \ c_{\lambda_i} 
		\ \text{s.t.} \ 
		\kvvbar{\mathbf{\pvi}_{\lambda_i}}_0 = \kvvbar{\mathbf{\pvi}^{\dagger}}_0
		.
	\end{equation}
}{, see \cite[Sec.~\ref{tecreport-sec:results}]{peeling2022technicalreport} for more details.} 

\subsection{Competing procedures}
We consider the following numerical solvers addressing~\eqref{prob:prob-base}: 
\textit{i)} \textsc{Cplex}~\cite{cplex2009v12}, a generic \gls{mip} solver; 
\textit{ii)} \textsc{L0bnb}~\cite{hazimeh2021sparse}, a standard \gls{bnb} procedure using a 
``breadth-first search'' exploration 
strategy, see \cite[Sec.~3.3]{hazimeh2021sparse};  
\textit{iii)} \textsc{Sbnb}, a standard \gls{bnb} procedure using a 
``depth-first search'' exploration strategy, see ~\cite[Sec.~2.2]{mhenni2020sparse};
\textit{iv)} \textsc{Sbnb-N}~\cite{guyard2022node}, corresponding to \textsc{Sbnb} enhanced with additional ``node-screening'' techniques;
\textit{v)} \textsc{Sbnb-P}, corresponding to \textsc{Sbnb} enhanced with the peeling strategy presented in this paper. 
\textsc{L0bnb}, \textsc{Sbnb}, \textsc{Sbnb-N} and \textsc{Sbnb-P} all use the same solving procedure for relaxed problem \eqref{eq:equivalent node relaxed problem:big M}, namely a coordinate descent method~\cite{wright2015coordinate}.  
We use the C++ implementation of \textsc{Cplex}\footnote{\protect{\url{https://github.com/jump-dev/CPLEX.jl}}}
and the Python implementation of \textsc{L0bnb}.\footnote{\protect{\url{https://github.com/hazimehh/L0Learn}}}
\textsc{Sbnb}, \textsc{Sbnb-N} and \textsc{Sbnb-P} are implemented in Julia. 
\footnote{\protect{\url{https://github.com/TheoGuyard/BnbPeeling.jl}}}

Regarding \textsc{Sbnb-P}, peeling is applied at each iteration of the numerical procedure solving the relaxed problem \eqref{eq:equivalent node relaxed problem:big M}.
We use the current iterate\addTG{,} say $\mathbf{\pvi}^{(k)}$\addTG{,} to define 
$\mathbf{\dvi} \triangleq \mathbf{y}-\mathbf{A} \mathbf{\pvi}^{(k)}$ 
and apply the peeling rules defined in \Cref{prop:entry-wise peeling} in parallel (\emph{i.e.,}\xspace simultaneously for all the components of $\mathbf{\pvi}$). The value of $\alpha$ satisfying \eqref{eq:entry-wise peeling treshold}, if any, is chosen as $\alpha=\UB{\alpha}+10^{-16}$. The peeled intervals are propagated through the decision tree as described in \Cref{sec:propagation}.
All the \gls{bnb} procedures 
are provided with the initial bounds $\mathbf{\bigLi}=-M\mathbf{1}$ and $\mathbf{\bigUi}=M\mathbf{1}$ for some proper value of $M$; this corresponds to the standard \textit{``Big-$M$''} constraint commonly considered in the literature \cite[Sec.~3]{bourguignon2015exact}. 
As far as our random simulation setup is concerned, it can be shown that \eqref{prob:prob-base} admits a unique minimizer $\opt{\mathbf{\pvi}}$ with probability one and  we thus choose $M= \gamma \norm{\opt{\mathbf{\pvi}}}{\infty}$ for some $\gamma\geq 1$ in our simulations.\footnote{This requires to solve \eqref{prob:prob-base} once beforehand to identify $\opt{\mathbf{\pvi}}$. This operation is here only done for the sake of comparing the sensibility of the solving methods to the choice of $\gamma$. } 
\ifthenelse{\boolean{techreport}}{}{
	More details are given in our companion paper~\cite[Sec.~\ref{tecreport-sec:results}]{peeling2022technicalreport}.
}

\ifthenelse{\boolean{techreport}}{
	In practice, we obtain $\opt{\mathbf{\pvi}}$ by solving a sequence of problems with an increasing value of $M$ in the \textit{Big-M} constraint.
	More specifically, letting $\opt{\mathbf{\pvi}}_{M}$ denote the solution of \eqref{prob:prob-base}  with the additional constraint ``$\|\mathbf{\pvi}\|_{\infty} \leq M$'', we are guaranteed that $\opt{\mathbf{\pvi}}=\opt{\mathbf{\pvi}}_{M}$ as soon as $\|\opt{\mathbf{\pvi}}_{M}\|_{\infty} < M$. 
	We thus compute the solution $\opt{\mathbf{\pvi}}_{M}$ for a 
	sequence of $M$ of the form 
	$\kset{(1.1)^iM_0}{i\in\kN}$ for some $M_0>0$ and stop as soon as $\|\opt{\mathbf{\pvi}}_{M}\|_{\infty} < M$.
}{} 

\subsection{Computational gains}

\begin{figure}[t]
    \input{img/sensibility}
	\vspace*{-1.5em}
    \caption{
		Left: Solving time as a function of $\gamma$  (top, $\sigma=1$) and $\sigma$ (bottom, $\gamma=5$). 
		Right: gain in terms of solving time (solid) and number of nodes explored (dashed) 
		with respect to 
		\textsc{Sbnb-N}.
	}
	\label{fig:sensibility}
\end{figure}

\Cref{fig:sensibility} presents the performance of the considered solving procedures. All results are averaged over 50 problem instances. Experiments were run on one Intel Xeon E5-2660 v3 CPU clocked at 2.60 GHz with 16 GB of RAM.
The left column in \Cref{fig:sensibility} represents the average solving time of each procedure as a function of $\gamma$ (top) and $\sigma$ (bottom);
the right column illustrates the gain allowed by the proposed method in terms of solving time (solid) and number of nodes explored (dashed) as compared to its best competitor, that is \textsc{Sbnb-N}.    

We note that \textsc{Sbnb-P} 
leads to the smallest running time in  all the considered setups. 
Since 
the latter 
corresponds to \textsc{Sbnb} where peeling has been added, the spacing between the red and green curves materializes the gain provided by peeling. 
As far as our simulation setup is concerned, we see that the proposed method enables an acceleration of almost one order of magnitude with respect to \textsc{Sbnb}. 
It is noticeable that this acceleration occurs even if $\gamma=1$, that is the \textit{Big-$M$} constraint is perfectly tuned to the problem at hand. This is due to the fact that peeling can refine \textit{individually} each component of the initial bounds $\mathbf{\bigLi}$ and $\mathbf{\bigUi}$ at \textit{each node} of the \gls{bnb} decision tree to fit the local geometry of the problem. 

We also notice that \textsc{Sbnb-P} improves over \textsc{Sbnb-N}, which can be seen as another acceleration of \textsc{Sbnb}. 
In particular, \textsc{Sbnb-P} performs always as well as \textsc{Sbnb-N} as emphasized by the gains in the right-hand side of Fig.~\ref{fig:sensibility}. 
We note in particular the gain provided by peeling in terms of number of nodes processed by the \gls{bnb} procedure: as expected, peeling allows for more aggressive pruning and thus reduces the number of nodes to be explored. 

\section{Notations}

We use the following notations.
\(\mathbf{0}\) and \(\mathbf{1}\) denote the all-zero and all-one vectors.
The \(\idxentry\)-th column of a matrix \(\mathbf{A}\) is denoted \(\mathbf{a}_\idxentry\) and the \(\idxentry\)-th entry of a vector \(\mathbf{\pvi}\) is denoted \(x_{\idxentry}\). Any vectorial relation has to be understood component-wise, \emph{e.g.,}\xspace $\mathbf{\pvi} \in [\mathbf{\bigLi},\mathbf{\bigUi}]$ means $x_\idxentry \in [l_\idxentry,u_\idxentry]$ $\forall \idxentry$.  
$\eta(\cdot)$ denotes the indicator function which equals to $0$ if the condition in argument is fulfilled and to $+\infty$ otherwise.
\(\pospart{x}=\max(x,0)\) refers to the positive-part function;  
\(\card{\cdot}\) denotes the cardinality of a set. $\intervint{1}{n}$ with \(n\in\kN*\) is a short-hand notation for the set $\{1,\ldots,n\}$. 

\section{Peeling}
\label{sec:peeling}


In this section, we introduce our proposed peeling procedure.
As an initial assumption, suppose that we know some interval $[\mathbf{\bigLi},\mathbf{\bigUi}]$ verifying condition \eqref{eq: condition equivalence}. This assumption will be relaxed later on in \Cref{sec:propagation}. 

Our goal is to find a new interval $[\bigL',\bigU']$ such that  
\begin{subequations}
	\begin{align}
		\forall \mathbf{\pvi}\in [\mathbf{\bigLi},\mathbf{\bigUi}]\setminus[\bigL',\bigU']&: \ \node{\mathrm{P}}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi})>\UB{p} 
		\label{eq:peeling: constraint}
		\\
		[\bigL',\bigU']&\subseteq [\mathbf{\bigLi},\mathbf{\bigUi}]
		\label{eq:peeling: constraint mnew <= m}
		.
	\end{align}	
\end{subequations}
These requirements imply that the pruning decision \eqref{eq:pruning condition} made at node $\nu$ remains unchanged when replacing  constraint ``$\mathbf{\pvi}\in [\mathbf{\bigLi},\mathbf{\bigUi}]$'' by ``$\mathbf{\pvi}\in [\bigL',\bigU']$'' in \eqref{eq:node problem constr}. More specifically, the following result holds: 
\begin{lemma}\label{lemma: consequence properties peeled interval}
 	Assume $[\mathbf{\bigLi},\mathbf{\bigUi}]$ and $[\bigL',\bigU']$ verify \eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m}, then 
 	\begin{align}\label{eq:peeling: equivalence}
 		\node{p}(\bigL',\bigU')> \UB{p}	
		&\iff
		\node{p}(\mathbf{\bigLi},\mathbf{\bigUi}) > \UB{p}.
 	\end{align}
\end{lemma} 
A proof of this result is available in \ifthenelse{\boolean{techreport}}{\Cref{proof:lemma:consequence properties peeled interval}}{\cite[App.~\ref{tecreport-proof:lemma:consequence properties peeled interval}]{peeling2022technicalreport}}.  
A consequence of preserving the pruning decision is that taking the new constraint ``$\mathbf{\pvi}\in [\bigL',\bigU']$'' into account at node $\nu$ does not alter the output of the \gls{bnb} procedure. 
In particular, it still correctly identifies the solutions of \eqref{prob:prob-base}.
The second requirement \eqref{eq:peeling: constraint mnew <= m} implies that $\node{r}(\bigL',\bigU')$ can possibly be larger than $\node{r}(\mathbf{\bigLi},\mathbf{\bigUi})$ since the lower bound in \eqref{eq:convex relaxation l0 pseudo norm} is tightened by considering larger absolute values for $\mathbf{\bigLi}$ and $\mathbf{\bigUi}$.
Overall, any choice of $[\bigL',\bigU']$ verifying \eqref{eq:peeling: constraint}\addTG{-}\eqref{eq:peeling: constraint mnew <= m} thus keeps unchanged the output of the \gls{bnb} procedure while allowing for potentially more aggressive (relaxed) pruning decisions. 

In the rest of this section, we describe a strategy to find some interval $[\bigL',\bigU']$ satisfying~\eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m}.  
Because of the symmetry of the problem at stake, we only focus on the construction of the upper bound $\bigU'$.
Nevertheless, the identification of a lower bound $\bigL'$ can be done along the same lines.  



\subsection{Target peeling strategy}

Given some index $\idxpeeled \in \setidx_\symbnone$ and $\alpha>0$, we consider the following perturbed versions of \eqref{eq:node problem constr}:
		\begin{align}
			\label{eq:peeling condition b}
			\node{p}_\alpha(\mathbf{\bigLi},\mathbf{\bigUi}) \triangleq \inf_{\mathbf{\pvi}\in\kR^n} 
			\node{\mathrm{P}}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi}) 
			+ \eta(x_{\idxpeeled} > \alpha).	
		\end{align}	
Problem \eqref{eq:peeling condition b} corresponds to~\eqref{eq:node problem constr} where $x_{\idxpeeled}$ is additionally constrained to be strictly greater than $\alpha$. 
The following lemma then trivially follows from the definition of $\node{p}_\alpha(\mathbf{\bigLi},\mathbf{\bigUi})$: 
	\begin{lemma}\label{lemma:ideal peeling test}
		If $\alpha \in [0,u_{\idxpeeled}[$ and 
		\begin{align}\label{eq:ideal peeling test}
			\node{p}_\alpha(\mathbf{\bigLi},\mathbf{\bigUi}) > \UB{p},
		\end{align}
		then \eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m} holds with 
		\begin{align} \label{eq:considered bigU'}
		u_\idxentry' 
		\,=\,
		\begin{cases}
			\alpha 			   & \mbox{if $\idxentry=\idxpeeled$}\\
			u_{\idxentry}   & \mbox{otherwise}.
		\end{cases}
		\end{align}
	\end{lemma}
In words, \Cref{lemma:ideal peeling test} states that 
 any $\alpha \in [0,u_{\idxpeeled}[$ verifying \eqref{eq:ideal peeling test} enables to  construct some $\bigU'$ automatically fulfilling~\eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m}. 
Unfortunately, evaluating~\eqref{eq:ideal peeling test} involves the same computational burden as solving~\eqref{eq:node problem constr}.
 This problem can nevertheless be circumvented by finding some proper lower bound on $\node{p}_\alpha(\mathbf{\bigLi},\mathbf{\bigUi})$ as described in the next section.


\subsection{Tractable implementation}

\ifthenelse{\boolean{techreport}}{In \Cref{proof:lower bounds} we show}{It can be shown (see~\cite[\Cref{tecreport-proof:lower bounds}]{peeling2022technicalreport})} that for any \(\mathbf{\dvi}\in\kR^m\), the following lower bound on $\node{p}_\alpha(\mathbf{\bigLi},\mathbf{\bigUi})$ holds: 

\begin{align}\label{eq:lower bound on peeling metric}
	\node{p}_\alpha(\mathbf{\bigLi},\mathbf{\bigUi})
	\geq \node{\mathrm{D}}(\mathbf{\dvi};\mathbf{\bigLi},\mathbf{\bigUi}) + \psi_\idxpeeled(\mathbf{\dvi};\mathbf{\bigLi},\mathbf{\bigUi})+ \alpha \pospart{-\ktranspose{\mathbf{a}}_{\idxpeeled}\mathbf{\dvi} }	
\end{align}
where 
\begin{align*}
	\node{\mathrm{D}}(\mathbf{\dvi};\mathbf{\bigLi},\mathbf{\bigUi}) 
	&\triangleq
	\tfrac{1}{2}\norm{\mathbf{y}}{2}^2 - \tfrac{1}{2}\norm{\mathbf{y} - \mathbf{\dvi}}{2}^2+ |\setidx_\symbone| \lambda \\
	&\qquad- \sum_{\idxentry\in\setidx_\symbone} 
	\pivot{0,\idxentry}(\ktranspose{\mathbf{a}}_{\idxentry}\mathbf{\dvi})
	- \sum_{\idxentry\in\setidx_\symbnone} \pivot{\lambda,\idxentry}(\ktranspose{\mathbf{a}}_{\idxentry}\mathbf{\dvi})
	\\
	\psi_\idxpeeled(\mathbf{\dvi};\mathbf{\bigLi},\mathbf{\bigUi})
	&\triangleq 
	\pivot{\lambda,\idxpeeled}(\ktranspose{\mathbf{a}}_{\idxpeeled}\mathbf{\dvi})
	- u_{\idxpeeled}\pospart{\ktranspose{\mathbf{a}}_{\idxpeeled}\mathbf{\dvi}}+ \lambda
\end{align*}
and $\pivot{\rho,\idxentry}(v)\triangleq \pospart{u_{\idxentry} v  - \rho } + \pospart{l_{\idxentry}v  -\rho}$.	



Using this result, condition \eqref{eq:ideal peeling test} can obviously be relaxed as
	\begin{align}\label{eq: alpha peeling condition}
		\node{\mathrm{D}}(\mathbf{\dvi};\mathbf{\bigLi},\mathbf{\bigUi}) + \psi_\idxpeeled(\mathbf{\dvi};\mathbf{\bigLi},\mathbf{\bigUi})+ \alpha \pospart{-\ktranspose{\mathbf{a}}_{\idxpeeled}\mathbf{\dvi} }	>\UB{p}
		.
	\end{align}
Hence, choosing any \(\alpha\in[0,u_\idxpeeled[\) verifying~\eqref{eq: alpha peeling condition} for some \(\mathbf{\dvi}\in\kR^m\) defines a new valid constraint (via \eqref{eq:considered bigU'}) in the sense of \eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m}. 
Interestingly, the left-hand side of~\eqref{eq: alpha peeling condition} depends linearly on \(\alpha\), thus allowing to precisely characterize the range of possible values for $\alpha$. 
This leads us to the main result of this section:
\begin{proposition}[Entry-wise peeling] \label{prop:entry-wise peeling}
	Let \(\mathbf{\dvi}\in\kR^m\).
	If \(\transp{\mathbf{a}}_{\idxpeeled}\mathbf{\dvi}\geq 0\) and 
	\begin{align}\label{eq:peeling: cas positif}
		\node{\mathrm{D}}(\mathbf{\dvi};\mathbf{\bigLi},\mathbf{\bigUi}) + \psi_\idxpeeled(\mathbf{\dvi};\mathbf{\bigLi},\mathbf{\bigUi})>\UB{p}
	\end{align}
	then $\bigU'$ defined as in~\eqref{eq:considered bigU'} with $\alpha=0$ fulfills~\eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m}.

	If \(\transp{\mathbf{a}}_{\idxpeeled}\mathbf{\dvi}<0\), then 
	$\bigU'$ defined as in~\eqref{eq:considered bigU'} with any 
	\(\alpha\in[0,u_{\idxpeeled}[\) verifying
	\begin{equation} \label{eq:entry-wise peeling treshold}
		\alpha
		>
		\UB{\alpha}\triangleq 
		\frac{
			\UB{p} - \node{\mathrm{D}}(\mathbf{\dvi};\mathbf{\bigLi},\mathbf{\bigUi}) - \psi_\idxpeeled(\mathbf{\dvi};\mathbf{\bigLi},\mathbf{\bigUi}) 
		}{
			\pospart{-\ktranspose{\mathbf{a}}_\idxpeeled\mathbf{\dvi}} 
		}
	\end{equation}
	fulfills~\eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m}.
\end{proposition}
Our next result shows that \Cref{prop:entry-wise peeling} can be applied to all indices \(\idxpeeled\in\intervint{1}{n}\) either sequentially or in parallel, while preserving the correctness of the \gls{bnb} procedure:
\begin{lemma} \label{eq:lemma parallele application}
	Let \([\mathbf{\bigLi}',\mathbf{\bigUi}']\) and \([\mathbf{\bigLi}'',\mathbf{\bigUi}'']\) be two intervals satisfying~\eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m}.
	Then, the interval \([\mathbf{\bigLi}',\mathbf{\bigUi}']\cap[\mathbf{\bigLi}'',\mathbf{\bigUi}'']\) also fulfills~\eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m}.
\end{lemma}
A proof is available in \ifthenelse{\boolean{techreport}}{\Cref{sec:proof lemma parallele application}}{\cite[App.~\ref{tecreport-sec:proof lemma parallele application}]{peeling2022technicalreport}}. 
We note that in terms of complexity the parallel 
application of \Cref{prop:entry-wise peeling} to all indices \(\idxpeeled\in\intervint{1}{n}\) requires the computation of the inner products $\{\transp{\mathbf{a}}_{\idxentry}\mathbf{\dvi}\}_{\idxentry=1}^{n}$ and \textit{one single} evaluation of \(\node{\mathrm{D}}(\mathbf{\dvi};\mathbf{\bigLi},\mathbf{\bigUi})\).
Interestingly\addCE{,} these inner products are already computed in most numerical procedures solving~\eqref{eq:equivalent node relaxed problem:big M} and are thus usually available at no additional cost, see \emph{e.g.,}\xspace~\cite[Sec.~4.3]{guyard2022node} for a discussion. 
The overhead complexity of our proposed peeling strategy thus scales as \(\mathcal{O}(n+m)\). 

\subsection{Propagating peeling down the tree}
\label{sec:propagation}

In this section, we emphasize that any interval $[\bigL',\bigU']$ verifying \eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m} at node $\nu$ can be used as a starting point to apply our peeling procedure at the child nodes of $\nu$. 
 More specifically, the following result holds: 

\begin{lemma}\label{corollary:heritage}
	Let $[\bigL',\bigU']$ be some interval verifying \eqref{eq:peeling: constraint}-\eqref{eq:peeling: constraint mnew <= m} at node $\nu$ and let $\nu'$ be some child node of $\nu$. 
	Assume that the peeling procedure defined in \Cref{prop:entry-wise peeling} is applied at node $\nu'$ with $[\bigL',\bigU']$  as input (rather than $[\mathbf{\bigLi},\mathbf{\bigUi}]$) to generate a new interval $[\bigL'',\bigU'']$.
	Then we have 
	\begin{subequations}
	\begin{align}
		\forall \mathbf{\pvi}\in[\mathbf{\bigLi},\mathbf{\bigUi}]\setminus [\bigL'',\bigU'']&:\,
		\mathrm{P}^{\nu'}(\mathbf{\pvi};\mathbf{\bigLi},\mathbf{\bigUi})>\UB{p}
		\label{eq:peeling: constraint 3}
		\\
		[\bigL'',\bigU'']&\subseteq [\mathbf{\bigLi},\mathbf{\bigUi}]
		\label{eq:peeling: constraint mnew <= m 3}
		.
	\end{align}
	\end{subequations}
\end{lemma}
A proof of this result is available in \ifthenelse{\boolean{techreport}}{\Cref{sec:proof corollary:heritage}}{\cite[App.~\ref{tecreport-sec:proof corollary:heritage}]{peeling2022technicalreport}}. 
In other words, \Cref{corollary:heritage} states that any peeled interval $[\bigL',\bigU']$ computed at node $\nu$ can be used as a starting point to apply a new peeling step at any child node $\nu'$.  This allows to propagate the peeled interval $[\bigL',\bigU']$ down the decision tree to hopefully improve sequentially the tightness of the convex relaxation \eqref{eq:equivalent node relaxed problem:big M}. 


\section{Additional numerical results}

\begin{figure*}
    \input{img/additional.tex}
    \caption{
        Three first columns: Solving time as a function of $\gamma$  (top, $\sigma=1$) and $\sigma$ (bottom, $\gamma=3$). 
        Last column: gain in terms of solving time (solid) and number of nodes explored (dashed) 
        with respect to 
        \textsc{Sbnb-N} for the three different types of hardness.
    }
    \label{fig:additional}
\end{figure*}

In \Cref{sec:results}, we only vary the parameters that may directly influence the \textit{Big-M} constraint.
To give a broader overview of the performance allowed by our peeling methodology, we have selected three different types of instances of \eqref{prob:prob-base} with increasing hardness.
We still fix $m$, $n$ and $\mathrm{SNR}$ but we vary $k$ and $\rho$ 
to control the instance inherent complexity.
On the one hand, larger $k$ leads to problems with a more complex combinatorial nature.
On the other hand, increasing $\rho$ leads to more correlation between the columns of $\mathbf{A}$, which brings the local minima of \eqref{prob:prob-base} closer to the global ones.
Both of these effects make the problem harder to solve.
Our three setups are constructed with the following combination of the parameters:
\begin{center}
    \begin{tabular}{cccccc}
        \toprule
        & $k$ & $m$ & $n$ & $\rho$ & SNR \\
        \midrule
        Easy & $5$ & $100$ & $150$ & $0.1$ & 15 \\
        Medium & $7$ & $100$ & $150$ & $0.1$ & 15 \\
        Hard & $7$ & $100$ & $150$ & $0.8$ & 15 \\
        \bottomrule
    \end{tabular}
\end{center}
We note that the maximum solution time allowed for the solver is 30 hours, that is just above $10^5$ seconds.
For particularly hard instances, some solvers hit this time limit.
Since it only concerns a minority of the trials, we have decided not to remove them, but this may have slightly lower their mean solution time.

\Cref{fig:additional} shows the solution time of the different methods for the three different types of instances when varying $\gamma$ and $\sigma$.
We roughly observe the same behaviors as in \Cref{fig:sensibility}. 
A notable remark is that the more difficult the instance, the larger the gains permitted by our peeling methodology.
Our empirical analysis is that on harder instances, the \gls{bnb} tree has to be explored more extensively.
Tightening the bounds at a given node will therefore allow a greater number of nodes to benefit from the strengthening of the relaxations.
The solving time will therefore be further reduced.
When varying $\gamma$, the gains of the method implementing our peeling strategy (\textsc{Sbnb-P}) again its best concurrent (\textsc{Sbnb-N}) almost reach a factor $2$ on \textsc{Easy} instances. 
On \textsc{Medium} instances, it reaches a factor $4$ and on \textsc{Hard} instances, it reaches a factor $5$. 
