<figure><img src="image_1.png" /></figure>

HAL
open science

# Straight-Through meets Sparse Recovery: the Support Exploration Algorithm

Mimoun Mohamed, François Malgouyres, Valentin Emiya, Caroline Chaux

► To cite this version:

Mimoun Mohamed, François Malgouyres, Valentin Emiya, Caroline Chaux. Straight-Through meets Sparse Recovery: the Support Exploration Algorithm. ICML 2024, the 41st International Conference on Machine Learning, Jul 2024, Vienna, Austria. hal-03964976v3

HAL Id: hal-03964976

https://hal.science/hal-03964976v3

Submitted on 24 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Straight-Through meets Sparse Recovery: the Support Exploration Algorithm

Mimoun Mohamed<fnref n="1" />,<fnref n="2" />, François Malgouyres<fnref n="4" />, Valentin Emiya<fnref n="1" />, and Caroline Chaux<fnref n="3" />

## Abstract

The *straight-through estimator* (STE) is commonly used to optimize quantized neural networks, yet its contexts of effective performance are still unclear despite empirical successes. To make a step forward in this comprehension, we apply STE to a well-understood problem: *sparse support recovery*. We introduce the *Support Exploration Algorithm* (SEA), a novel algorithm promoting sparsity, and we analyze its performance in support recovery (a.k.a. model selection) problems. SEA explores more supports than the state-of-the-art, leading to superior performance in experiments, especially when the columns of A are strongly coherent. The theoretical analysis considers recovery guarantees when the linear measurements matrix A satisfies the *Restricted Isometry Property* (RIP). The sufficient conditions of recovery are comparable but more stringent than those of the state-of-the-art in sparse support recovery. Their significance lies mainly in their applicability to an instance of the STE.

# 1 Introduction

**Straight-through estimator.** The use of quantized neural networks spares memory, energy, and computing resources during inference, making them essential for embedding neural networks [47, 40]. An effective strategy is to learn the quantized weights. Seminal works [15, 30] rely on full-precision weights $w$ that evolve in the parameter space, while quantized weights $w_q = H_q(w)$ are obtained by applying a piecewise-constant quantization operator $H_q$.

Denoting by $F$ the computational chain from $w_q$ to the loss, the learning procedure Minimize $F(H_q(w))$ relies on the computational graph:

<figure><img src="image_1.png" /></figure>

Given a step-size $\eta$, the update of $w$ is performed as $w \leftarrow w - \eta \frac{\partial F}{\partial w}|_{w_q}$. The motivation, as explained in [29, 5], is that since $\frac{\partial H_q}{\partial w}|_w$ is either undefined or 0, we cannot backpropagate using the chain rule $\frac{\partial F \circ H_q}{\partial w}|_w = \frac{\partial F}{\partial w_q}|_{w_q} \frac{\partial H_q}{\partial w}|_w$. The STE makes the coarse approximation $\frac{\partial F \circ H_q}{\partial w}|_w \approx \frac{\partial F}{\partial w_q}|_{w_q}$ to backpropagate the gradient through the piecewise-constant operator $H_q$. Many subsequent works improve these methods in various aspects [47, 40].

Although STE achieves state-of-the-art performance in training quantized weights for neural networks, it is poorly understood and has not been investigated beyond the context of quantization. We introduce an STE principle for sparsification, leading to a novel algorithm named the Support Exploration Algorithm (SEA) and present experimental evidence of its benefits in challenging, coherent settings such as spike deconvolution, as well as in systematic experiments like the phase transition diagram, see Figure 1. Additionally, we establish theoretical guarantees for the STE-based algorithm.

**Sparse support recovery.** For a sparsity $k \in \mathbb{N}$, we assume $x^* \in \mathbb{R}^n$ is an unknown sparse vector of unknown support $S^* = \text{supp}(x^*)$, of sparsity $|S^*| \le k$, $A \in \mathbb{R}^{m \times n}$ is a known matrix, and

$$y = Ax^* + e \in \mathbb{R}^m$$

<footnote n="1">¹Aix Marseille Univ, CNRS, LIS, Marseille, France</footnote>

<footnote n="2">²Aix Marseille Univ, CNRS, I2M, Marseille, France</footnote>

<footnote n="3">³CNRS, IPAL, Singapour</footnote>

<footnote n="4">⁴Institut de Mathématiques de Toulouse ; UMR5219, Université de Toulouse ; CNRS, UPS IMT F-31062 Toulouse Cedex 9, France</footnote>
# Straight-Through meets Sparse Recovery: the Support Exploration Algorithm

Mimoun Mohamed<fnref n="1" />,<fnref n="2" />, François Malgouyres<fnref n="4" />, Valentin Emiya<fnref n="1" />, and Caroline Chaux<fnref n="3" />

## Abstract

The *straight-through estimator* (STE) is commonly used to optimize quantized neural networks, yet its contexts of effective performance are still unclear despite empirical successes. To make a step forward in this comprehension, we apply STE to a well-understood problem: *sparse support recovery*. We introduce the *Support Exploration Algorithm* (SEA), a novel algorithm promoting sparsity, and we analyze its performance in support recovery (a.k.a. model selection) problems. SEA explores more supports than the state-of-the-art, leading to superior performance in experiments, especially when the columns of A are strongly coherent. The theoretical analysis considers recovery guarantees when the linear measurements matrix A satisfies the *Restricted Isometry Property* (RIP). The sufficient conditions of recovery are comparable but more stringent than those of the state-of-the-art in sparse support recovery. Their significance lies mainly in their applicability to an instance of the STE.

# 1 Introduction

**Straight-through estimator.** The use of quantized neural networks spares memory, energy, and computing resources during inference, making them essential for embedding neural networks [47, 40]. An effective strategy is to learn the quantized weights. Seminal works [15, 30] rely on full-precision weights $w$ that evolve in the parameter space, while quantized weights $w_q = H_q(w)$ are obtained by applying a piecewise-constant quantization operator $H_q$.

Denoting by $F$ the computational chain from $w_q$ to the loss, the learning procedure Minimize $F(H_q(w))$ relies on the computational graph:

<figure><img src="image_2.png" /></figure>

Given a step-size $\eta$, the update of $w$ is performed as $w \leftarrow w - \eta \frac{\partial F}{\partial w}|_{w_q}$. The motivation, as explained in [29, 5], is that since $\frac{\partial H_q}{\partial w}|_w$ is either undefined or 0, we cannot backpropagate using the chain rule $\frac{\partial F \circ H_q}{\partial w}|_w = \frac{\partial F}{\partial w_q}|_{w_q} \frac{\partial H_q}{\partial w}|_w$. The STE makes the coarse approximation $\frac{\partial F \circ H_q}{\partial w}|_w \approx \frac{\partial F}{\partial w_q}|_{w_q}$ to backpropagate the gradient through the piecewise-constant operator $H_q$. Many subsequent works improve these methods in various aspects [47, 40].

Although STE achieves state-of-the-art performance in training quantized weights for neural networks, it is poorly understood and has not been investigated beyond the context of quantization. We introduce an STE principle for sparsification, leading to a novel algorithm named the Support Exploration Algorithm (SEA) and present experimental evidence of its benefits in challenging, coherent settings such as spike deconvolution, as well as in systematic experiments like the phase transition diagram, see Figure 1. Additionally, we establish theoretical guarantees for the STE-based algorithm.

**Sparse support recovery.** For a sparsity $k \in \mathbb{N}$, we assume $x^* \in \mathbb{R}^n$ is an unknown sparse vector of unknown support $S^* = \text{supp}(x^*)$, of sparsity $|S^*| \le k$, $A \in \mathbb{R}^{m \times n}$ is a known matrix, and

$$y = Ax^* + e \in \mathbb{R}^m$$

<footnote n="1">¹Aix Marseille Univ, CNRS, LIS, Marseille, France</footnote>

<footnote n="2">²Aix Marseille Univ, CNRS, I2M, Marseille, France</footnote>

<footnote n="3">³CNRS, IPAL, Singapour</footnote>

<footnote n="4">⁴Institut de Mathématiques de Toulouse ; UMR5219, Université de Toulouse ; CNRS, UPS IMT F-31062 Toulouse Cedex 9, France</footnote>
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Figure 1: Overview of the main results. Left: phase transition diagram showing the recovery limits in dimension $n = 500$ while sparsity $k$ and number of observations $m$ varies (the higher, the better, see details in Section 5.1). Right: spike deconvolution in dimension $m = n = 500$ - Average distance between the supports of the solution $x^*$ and the estimations obtained from various algorithms, plotted against the sparsity level $k$ (the lower, the better, see details in Section 5.2).</figcaption></figure>

is a linear observation of $x^*$, contaminated with an additive error/noise $e \in \mathbb{R}^m$.

The support recovery objective<fnref n="1" />, also coined variable or model selection, searches for a support $S$ with cardinality at most $k$ such that $S^* \subseteq S$. We say that an algorithm recovers $S^*$ if it finds such an $S$. Using a least-square criterion

$$F(x) = \frac{1}{2} \|Ax - y\|_2^2,$$

a famous model for support recovery is the optimization problem with sparsity constraint:

$$\text{Minimize } F(x) \text{ s.t. } \|x\|_0 \le k, \qquad (1)$$

where $\|x\|_0$ is the $l^0$ pseudo-norm of $x$. Problem (1) is known to be NP-hard [17] and to recover the correct support under mild conditions [24, Chapter 5.2.2].

**Proposed STE-based approach for sparse recovery.** Let us define an unconstrained optimization problem which is equivalent to problem (1) and whose structure is compatible with the STE. We set

$$\text{Minimize } F(H(\mathcal{X})) \qquad (2)$$

where $H$ is the sparsification operator<fnref n="2" />

$$H(\mathcal{X}) \in \underset{\substack{x \in \mathbb{R}^n \\ \text{supp}(x) \subseteq \text{largest}_k(\mathcal{X})}}{\text{argmin}} \frac{1}{2} \|Ax - y\|_2^2. \qquad (3)$$

The equivalence between (1) and (2) is established in Appendix A.

Operator $H$ is piece-wise constant and finds the non-zero values in $x$, the sparse support being induced by dense vector $\mathcal{X}$. Formulation (2) has a similar structure as in the quantization case presented above. It uses the suitable loss $F$ and the sparsification operator $H$ in place of the quantification operator $H_q$. Here, vector $\mathcal{X}$ is dense and $x = H(\mathcal{X})$ is $k$-sparse. Applying the STE to the new formulation (2), we obtain the update $\mathcal{X} \leftarrow \mathcal{X} - \eta \frac{\partial F}{\partial x}|_{H(\mathcal{X})}$, for a step-size $\eta > 0$, where we have approximated $\frac{\partial F \circ H}{\partial \mathcal{X}}|_{\mathcal{X}} = \frac{\partial F}{\partial x}|_{H(\mathcal{X})} \frac{\partial H}{\partial \mathcal{X}}|_{\mathcal{X}} \approx \frac{\partial F}{\partial x}|_{H(\mathcal{X})}$. This leads to an original algorithm for sparse recovery, based on the STE, which we analyze in this article.

**Contributions.** The first contribution of the article is to adapt the STE to the sparse support recovery problem (as explained above). Doing so, we obtain a new sparsity-inducing algorithm that we call Support Exploration Algorithm (SEA). It uses the full gradient history over iterations as a heuristic in order to select the next support to optimize over. SEA is supported by support recovery guarantees. In Theorem 4.1 and Corollary 4.2, the sufficient hypotheses guaranteeing the support recovery are on the Restricted Isometry Property (RIP) constants of $A$ and $x^*$. These conditions are comparable to those in the state-of-the-art,

<footnote n="1">¹The adaptation of the article to “signed support recovery” is possible and is straightforward. We chose to simplify the presentation and not discuss sign recovery.</footnote>

<footnote n="2">²In scenarios of interest, the minimization problem (3) has a unique and easy to compute solution.</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Figure 1: Overview of the main results. Left: phase transition diagram showing the recovery limits in dimension $n = 500$ while sparsity $k$ and number of observations $m$ varies (the higher, the better, see details in Section 5.1). Right: spike deconvolution in dimension $m = n = 500$ - Average distance between the supports of the solution $x^*$ and the estimations obtained from various algorithms, plotted against the sparsity level $k$ (the lower, the better, see details in Section 5.2).</figcaption></figure>

is a linear observation of $x^*$, contaminated with an additive error/noise $e \in \mathbb{R}^m$.

The support recovery objective<fnref n="1" />, also coined variable or model selection, searches for a support $S$ with cardinality at most $k$ such that $S^* \subseteq S$. We say that an algorithm recovers $S^*$ if it finds such an $S$. Using a least-square criterion

$$F(x) = \frac{1}{2} \|Ax - y\|_2^2,$$

a famous model for support recovery is the optimization problem with sparsity constraint:

$$\text{Minimize } F(x) \text{ s.t. } \|x\|_0 \le k, \qquad (1)$$

where $\|x\|_0$ is the $l^0$ pseudo-norm of $x$. Problem (1) is known to be NP-hard [17] and to recover the correct support under mild conditions [24, Chapter 5.2.2].

**Proposed STE-based approach for sparse recovery.** Let us define an unconstrained optimization problem which is equivalent to problem (1) and whose structure is compatible with the STE. We set

$$\text{Minimize } F(H(\mathcal{X})) \qquad (2)$$

where $H$ is the sparsification operator<fnref n="2" />

$$H(\mathcal{X}) \in \underset{\substack{x \in \mathbb{R}^n \\ \text{supp}(x) \subseteq \text{largest}_k(\mathcal{X})}}{\text{argmin}} \frac{1}{2} \|Ax - y\|_2^2. \qquad (3)$$

The equivalence between (1) and (2) is established in Appendix A.

Operator $H$ is piece-wise constant and finds the non-zero values in $x$, the sparse support being induced by dense vector $\mathcal{X}$. Formulation (2) has a similar structure as in the quantization case presented above. It uses the suitable loss $F$ and the sparsification operator $H$ in place of the quantification operator $H_q$. Here, vector $\mathcal{X}$ is dense and $x = H(\mathcal{X})$ is $k$-sparse. Applying the STE to the new formulation (2), we obtain the update $\mathcal{X} \leftarrow \mathcal{X} - \eta \frac{\partial F}{\partial \mathcal{X}}|_{H(\mathcal{X})}$, for a step-size $\eta > 0$, where we have approximated $\frac{\partial F \circ H}{\partial \mathcal{X}}|_{\mathcal{X}} = \frac{\partial F}{\partial x}|_{H(\mathcal{X})} \frac{\partial H}{\partial \mathcal{X}}|_{\mathcal{X}} \approx \frac{\partial F}{\partial x}|_{H(\mathcal{X})}$. This leads to an original algorithm for sparse recovery, based on the STE, which we analyze in this article.

**Contributions.** The first contribution of the article is to adapt the STE to the sparse support recovery problem (as explained above). Doing so, we obtain a new sparsity-inducing algorithm that we call *Support Exploration Algorithm (SEA)*. It uses the full gradient history over iterations as a heuristic in order to select the next support to optimize over. SEA is supported by support recovery guarantees. In Theorem 4.1 and Corollary 4.2, the sufficient hypotheses guaranteeing the support recovery are on the Restricted Isometry Property (RIP) constants of $A$ and $x^*$. These conditions are comparable to those in the state-of-the-art,

<footnote n="1">¹The adaptation of the article to “signed support recovery” is possible and is straightforward. We chose to simplify the presentation and not discuss sign recovery.</footnote>

<footnote n="2">²In scenarios of interest, the minimization problem (3) has a unique and easy to compute solution.</footnote>
albeit slightly more stringent. Their interest mainly lies in the fact that they apply to an instance of
STE, for which very few guarantees of convergence exist. However, the successes of SEA observed in the
experiments extend to coherent problems where the RIP hypothesis is no longer satisfied. Additional
support recovery statements are in Theorem C.4 and in Corollary C.7 of Appendix C. The proofs are based
on the interpretation of SEA as a noisy version of an ‘Oracle algorithm’ which is analyzed in Appendix C.1.

The performances of SEA are compared to those of state-of-the-art algorithms on: 1/ phase-transition synthetic experiments for Gaussian matrices; 2/ spike deconvolution problems; 3/ classification and regression problems for real datasets. An important feature of SEA is that it can be used as a post-processing to improve the results of existing algorithms, as shown in the experiments. Also, because SEA has the ability to explore more supports, it performs remarkably well when the matrix $A$ is coherent. The code is available in the git repository of the project.<fnref n="3" />

**Organization of the article.** Related works are detailed in Section 2. Then, SEA is described in Section 3. The theoretical analysis is in Section 4. The experiments are in Section 5. Conclusions and perspectives are in Section 6.

In Appendices, a thorough comparison between SEA and the most similar algorithms of the state-of-the-art is detailed in Appendix B. Appendix B also details an efficient implementation of SEA. The proofs of the theoretical statements as well as complementary support-recovery statements are in Appendix C. Complementary experimental results are in Appendices D, E, and F.

## 2 Related Works

**On the STE.** Although STE achieves performances defining the state-of-the-art, it is poorly understood. In [34], the authors show that STE behaves well on convex problems and that a stochastic variant of STE does not on non-convex ones. For a two-layer linear neural network with a quantized activation function, a well-chosen STE converges to a critical point of the population risk [46] or reproduces a teacher network [35]. These are the only known formal guarantees for STE. Other works study the large dimension geometry of binary weights [1]. In [28] the authors interpret $w$ as an inertia variable and design a new (related) algorithm. In [14], the authors view the STE as a projected Wasserstein gradient flow.

**On sparse prior and support recovery.** Sparse representations and sparsity-inducing algorithms are widely used in statistics and machine learning [27], as well as in signal processing [24]. For instance, in machine learning, sparse representations are used to select relevant variables. They are also sought to interpret trained models. In signal processing, linear inverse problems have a wide array of applications. The sparsity assumption is ubiquitous since most real signals can be exactly or approximately represented as sparse signals in some domains, e.g., communication signals in Fourier space, natural images in wavelet space. While sparse models are appealing, they are hard to estimate due to the underlying combinatorial difficulty of identifying the correct sparse support.

Our algorithm has been designed in a support recovery context. In the noisy case $e \neq 0$, support recovery is a stronger guarantee than the one in the most standard compressed sensing setting, initiated in [11, 20], when the goal is to upper-bound $\|x - x^*\|_2$, for a well-chosen $x$. The first particularity of support recovery is to assume $x^*$ is truly $k$-sparse – not just compressible. Also, support recovery guarantees always involve a hypothesis on $\min_{i \in S^*} |x_i^*|$, in addition of the incoherence hypothesis on $A$ [45, 37, 49, 10, 48]. We cannot indeed expect to recover an element $i \in S^*$ if $|x_i^*|$ is negligible when compared to all the other quantities involved in the problem [45].

**Support recovery models and algorithms.** Beyond (1), various algorithms were investigated. There are three main families of algorithms: relaxation, combinatorial approaches, and greedy algorithms.

The most famous relaxed model uses the $ℓ^1$ norm and is known as the LASSO [42] or Basis Pursuit Algorithm [13]. Combinatorial approaches like Branch and Bound algorithms [4], find the global minimum of (1) but lack scalability. Greedy algorithms can be divided into two categories. Greedy Pursuits like Matching Pursuit (MP) [36] and Orthogonal Matching Pursuit (OMP) [39] are algorithms that start from an empty support and build up an estimate of $x^*$ by iteratively adding components to the current support and optimizing the components. As for thresholding algorithms like Iterative Hard Thresholding (IHT) [8], Normalized Iterative Hard Thresholding (NIHT) [9], Hard Thresholding Pursuit (HTP) [25], Compressive Sampling Matching Pursuit (CoSaMP) [38], OMP with Replacement (OMPR) [32], Exhaustive Local

<footnote n="3">³https://gitlab.lis-lab.fr/valentin.emiya/sea-icml-2024</footnote>
albeit slightly more stringent. Their interest mainly lies in the fact that they apply to an instance of
STE, for which very few guarantees of convergence exist. However, the successes of SEA observed in the
experiments extend to coherent problems where the RIP hypothesis is no longer satisfied. Additional
support recovery statements are in Theorem C.4 and in Corollary C.7 of Appendix C. The proofs are based
on the interpretation of SEA as a noisy version of an ‘Oracle algorithm’ which is analyzed in Appendix C.1.

The performances of SEA are compared to those of state-of-the-art algorithms on: 1/ phase-transition synthetic experiments for Gaussian matrices; 2/ spike deconvolution problems; 3/ classification and regression problems for real datasets. An important feature of SEA is that it can be used as a post-processing to improve the results of existing algorithms, as shown in the experiments. Also, because SEA has the ability to explore more supports, it performs remarkably well when the matrix $A$ is coherent. The code is available in the git repository of the project.<fnref n="3" />

**Organization of the article.** Related works are detailed in Section 2. Then, SEA is described in Section 3. The theoretical analysis is in Section 4. The experiments are in Section 5. Conclusions and perspectives are in Section 6.

In Appendices, a thorough comparison between SEA and the most similar algorithms of the state-of-the-art is detailed in Appendix B. Appendix B also details an efficient implementation of SEA. The proofs of the theoretical statements as well as complementary support-recovery statements are in Appendix C. Complementary experimental results are in Appendices D, E, and F.

## 2 Related Works

**On the STE.** Although STE achieves performances defining the state-of-the-art, it is poorly understood. In [34], the authors show that STE behaves well on convex problems and that a stochastic variant of STE does not on non-convex ones. For a two-layer linear neural network with a quantized activation function, a well-chosen STE converges to a critical point of the population risk [46] or reproduces a teacher network [35]. These are the only known formal guarantees for STE. Other works study the large dimension geometry of binary weights [1]. In [28] the authors interpret $w$ as an inertia variable and design a new (related) algorithm. In [14], the authors view the STE as a projected Wasserstein gradient flow.

**On sparse prior and support recovery.** Sparse representations and sparsity-inducing algorithms are widely used in statistics and machine learning [27], as well as in signal processing [24]. For instance, in machine learning, sparse representations are used to select relevant variables. They are also sought to interpret trained models. In signal processing, linear inverse problems have a wide array of applications. The sparsity assumption is ubiquitous since most real signals can be exactly or approximately represented as sparse signals in some domains, e.g., communication signals in Fourier space, natural images in wavelet space. While sparse models are appealing, they are hard to estimate due to the underlying combinatorial difficulty of identifying the correct sparse support.

Our algorithm has been designed in a support recovery context. In the noisy case $e \neq 0$, support recovery is a stronger guarantee than the one in the most standard compressed sensing setting, initiated in [11, 20], when the goal is to upper-bound $\|x - x^*\|_2$, for a well-chosen $x$. The first particularity of support recovery is to assume $x^*$ is truly $k$-sparse – not just compressible. Also, support recovery guarantees always involve a hypothesis on $\min_{i \in S^*} |x_i^*|$, in addition of the incoherence hypothesis on $A$ [45, 37, 49, 10, 48]. We cannot indeed expect to recover an element $i \in S^*$ if $|x_i^*|$ is negligible when compared to all the other quantities involved in the problem [45].

**Support recovery models and algorithms.** Beyond (1), various algorithms were investigated. There are three main families of algorithms: relaxation, combinatorial approaches, and greedy algorithms.

The most famous relaxed model uses the $ℓ^1$ norm and is known as the LASSO [42] or Basis Pursuit Algorithm [13]. Combinatorial approaches like Branch and Bound algorithms [4], find the global minimum of (1) but lack scalability. Greedy algorithms can be divided into two categories. Greedy Pursuits like Matching Pursuit (MP) [36] and Orthogonal Matching Pursuit (OMP) [39] are algorithms that start from an empty support and build up an estimate of $x^*$ by iteratively adding components to the current support and optimizing the components. As for thresholding algorithms like Iterative Hard Thresholding (IHT) [8], Normalized Iterative Hard Thresholding (NIHT) [9], Hard Thresholding Pursuit (HTP) [25], Compressive Sampling Matching Pursuit (CoSaMP) [38], OMP with Replacement (OMPR) [32], Exhaustive Local

<footnote n="3">³https://gitlab.lis-lab.fr/valentin.emiya/sea-icml-2024</footnote>
Search (ELS) [2] (a.k.a. Fully Corrective Forward Greedy Selection with Replacement [41]) and Subspace Pursuit (SP) [16], they start from any vector and add a replacement step in the iterative process. It allows them to explore various supports before stopping at a local optimum.

**Position of the article.** In this work, we take a different approach and apply the STE to a well-understood problem to compare its behavior, empirical performances, and theoretical guarantees to those of the well-established state-of-the-art.

Compared to other sparse support recovery algorithms, the algorithm introduced in this article may belong to the family of greedy algorithms. A clear difference is the introduction of a non-sparse vector $\mathcal{X}^t \in \mathbb{R}^n$, which evolves during the iterative process and indicates which support should be tested at iteration $t$. We call $\mathcal{X}^t$ the support exploration variable. It is the analog of the full-precision weights – used by BinaryConnect that also instantiates the STE – to optimize binary weights of neural networks [15, 31]. We exhibit that the adaptation of STE for sparsification enables a different exploration/exploitation trade-off compared to the state-of-the-art. It explores more. This permits to obtain better performances than the state-of-the-art on very difficult –coherent– problems. We establish that it is possible, in the sparse support recovery context, to obtain theoretical guarantees for the STE.

# 3 Method

After clarifying the notations in Section 3.1, SEA is described in detail in Section 3.2 and its computational complexity is discussed in Section 3.3.

## 3.1 Notations

For any $a, b \in \mathbb{R}$ ($a$ and $b$ can be real numbers), the set of integers between $a$ and $b$ is denoted by $\llbracket a, b \rrbracket$ and $\lceil a \rceil$ denotes the floor of $a$. For any set $S \subseteq \llbracket 1, n \rrbracket$, we denote the cardinality of $S$ by $|S|$. The complement of $S$ in $\llbracket 1, n \rrbracket$ is denoted by $\bar{S}$.

Given $x \in \mathbb{R}^n$ and $i \in \llbracket 1, n \rrbracket$, the $i^{th}$ entry of $x$ is denoted by $x_i$. The support of $x$ is denoted by $\text{supp}(x) = \{i : x_i \neq 0\}$. The $\ell^0$ pseudo-norm of $x$ is defined by $\|x\|_0 = |\text{supp}(x)|$. The set containing the indices of the $k$ largest absolute entries of $x$ is denoted by $\text{largest}_k(x)$. When ties lead to multiple possible choices for $\text{largest}_k(x)$, we select the solution with the highest indices. For instance, $\text{largest}_k(0) = \llbracket n - k + 1, n \rrbracket$.

For any $S \subseteq \llbracket 1, n \rrbracket$, $A \in \mathbb{R}^{m \times n}$, and $x \in \mathbb{R}^n$, we define $x|_S \in \mathbb{R}^{|S|}$, the restriction of the vector $x$ to the indices in $S$, and $A_S \in \mathbb{R}^{m \times |S|}$, the restriction of the matrix $A$ to the set $S$ as the matrix composed of the columns of $A$ whose indexes are in $S$. The transpose of $A$ is denoted by $A^T \in \mathbb{R}^{n \times m}$. The pseudoinverse of $A$ is denoted by $A^\dagger \in \mathbb{R}^{n \times m}$ and the pseudoinverse of $A_S$ by $A_S^\dagger = (A_S)^\dagger \in \mathbb{R}^{|S| \times m}$. For any $d \in \mathbb{N}$, the identity matrix of size $d$ is denoted by $I_d$. The symbol $\odot$ denotes the Hadamard product.

## 3.2 The Support Exploration Algorithm

The proposed Support Exploration Algorithm (SEA) is given in Algorithm 1. In terms of pseudocode, SEA resembles many state-of-the-art algorithms and is close to HTP and IHT (see comparison between SEA, HTP, and IHT in Appendix B). However, it stands out from the others for its exploratory behavior, which stems from the STE principle behind it.

Each iteration begins by the forward pass given by (3), in which the current sparse solution $x^t$ is computed by applying the sparsification operator $H$ to a dense vector $\mathcal{X}^t$: after selecting the support $S^t$ from $\mathcal{X}^t$ at line 6, the sparse solution $x^t = H(\mathcal{X}^t)$ is computed at line 7. The backward pass uses the STE principle

$$ \frac{\partial F \circ H}{\partial \mathcal{X}}|_{\mathcal{X}^t} \approx \frac{\partial F}{\partial x}|_{\mathcal{X}^t} = A^T(Ax^t - y), $$

at line 8. To illustrate with the analogy with BinaryConnect [15], the non-sparse vector $\mathcal{X}$ is the analog of the full-precision weights and $H(\mathcal{X})$ is the analog of the quantized weights. To the best of our knowledge, this is the first use of the STE to solve a sparse linear inverse problem.

The key idea is that support $S^t$ is designated at line 6 by a non-sparse variable $\mathcal{X}^t$ called the support exploration variable. It offers an original mechanism to explore supports in a more diverse way than existing algorithms. Variable $\mathcal{X}^{t+1} = \mathcal{X}^0 - \eta \sum_{t'=0}^t A^T(Ax^{t'} - y)$ is actually an accumulation of gradients taken in the sparse iterates and is used to designate the support of the next sparse iterate $x^{t+1}$. Consequently, unlike other descent-based algorithms, $\mathcal{X}^t$ is not confined to the neighborhood of $k$-sparse vectors. Its
Search (ELS) [2] (a.k.a. Fully Corrective Forward Greedy Selection with Replacement [41]) and Subspace Pursuit (SP) [16], they start from any vector and add a replacement step in the iterative process. It allows them to explore various supports before stopping at a local optimum.

**Position of the article.** In this work, we take a different approach and apply the STE to a well-understood problem to compare its behavior, empirical performances, and theoretical guarantees to those of the well-established state-of-the-art.

Compared to other sparse support recovery algorithms, the algorithm introduced in this article may belong to the family of greedy algorithms. A clear difference is the introduction of a non-sparse vector $\mathcal{X}^t \in \mathbb{R}^n$, which evolves during the iterative process and indicates which support should be tested at iteration $t$. We call $\mathcal{X}^t$ the support exploration variable. It is the analog of the full-precision weights – used by BinaryConnect that also instantiates the STE – to optimize binary weights of neural networks [15, 31]. We exhibit that the adaptation of STE for sparsification enables a different exploration/exploitation trade-off compared to the state-of-the-art. It explores more. This permits to obtain better performances than the state-of-the-art on very difficult –coherent– problems. We establish that it is possible, in the sparse support recovery context, to obtain theoretical guarantees for the STE.

# 3 Method

After clarifying the notations in Section 3.1, SEA is described in detail in Section 3.2 and its computational complexity is discussed in Section 3.3.

## 3.1 Notations

For any $a, b \in \mathbb{R}$ ($a$ and $b$ can be real numbers), the set of integers between $a$ and $b$ is denoted by $\llbracket a, b \rrbracket$ and $\lceil a \rceil$ denotes the floor of $a$. For any set $S \subseteq \llbracket 1, n \rrbracket$, we denote the cardinality of $S$ by $|S|$. The complement of $S$ in $\llbracket 1, n \rrbracket$ is denoted by $\bar{S}$.

Given $x \in \mathbb{R}^n$ and $i \in \llbracket 1, n \rrbracket$, the $i^{th}$ entry of $x$ is denoted by $x_i$. The support of $x$ is denoted by $\text{supp}(x) = \{i : x_i \neq 0\}$. The $\ell^0$ pseudo-norm of $x$ is defined by $\|x\|_0 = |\text{supp}(x)|$. The set containing the indices of the $k$ largest absolute entries of $x$ is denoted by $\text{largest}_k(x)$. When ties lead to multiple possible choices for $\text{largest}_k(x)$, we select the solution with the highest indices. For instance, $\text{largest}_k(0) = \llbracket n - k + 1, n \rrbracket$.

For any $S \subseteq \llbracket 1, n \rrbracket$, $A \in \mathbb{R}^{m \times n}$, and $x \in \mathbb{R}^n$, we define $x|_S \in \mathbb{R}^{|S|}$, the restriction of the vector $x$ to the indices in $S$, and $A_S \in \mathbb{R}^{m \times |S|}$, the restriction of the matrix $A$ to the set $S$ as the matrix composed of the columns of $A$ whose indexes are in $S$. The transpose of $A$ is denoted by $A^T \in \mathbb{R}^{n \times m}$. The pseudoinverse of $A$ is denoted by $A^\dagger \in \mathbb{R}^{n \times m}$ and the pseudoinverse of $A_S$ by $A_S^\dagger = (A_S)^\dagger \in \mathbb{R}^{|S| \times m}$. For any $d \in \mathbb{N}$, the identity matrix of size $d$ is denoted by $I_d$. The symbol $\odot$ denotes the Hadamard product.

## 3.2 The Support Exploration Algorithm

The proposed Support Exploration Algorithm (SEA) is given in Algorithm 1. In terms of pseudocode, SEA resembles many state-of-the-art algorithms and is close to HTP and IHT (see comparison between SEA, HTP, and IHT in Appendix B). However, it stands out from the others for its exploratory behavior, which stems from the STE principle behind it.

Each iteration begins by the forward pass given by (3), in which the current sparse solution $x^t$ is computed by applying the sparsification operator $H$ to a dense vector $\mathcal{X}^t$: after selecting the support $S^t$ from $\mathcal{X}^t$ at line 6, the sparse solution $x^t = H(\mathcal{X}^t)$ is computed at line 7. The backward pass uses the STE principle

$$ \frac{\partial F \circ H}{\partial \mathcal{X}}|_{\mathcal{X}^t} \approx \frac{\partial F}{\partial x}|_{\mathcal{X}^t} = A^T(Ax^t - y), $$

at line 8. To illustrate with the analogy with BinaryConnect [15], the non-sparse vector $\mathcal{X}$ is the analog of the full-precision weights and $H(\mathcal{X})$ is the analog of the quantized weights. To the best of our knowledge, this is the first use of the STE to solve a sparse linear inverse problem.

The key idea is that support $S^t$ is designated at line 6 by a non-sparse variable $\mathcal{X}^t$ called the support exploration variable. It offers an original mechanism to explore supports in a more diverse way than existing algorithms. Variable $\mathcal{X}^{t+1} = \mathcal{X}^0 - \eta \sum_{t'=0}^t A^T(Ax^{t'} - y)$ is actually an accumulation of gradients taken in the sparse iterates and is used to designate the support of the next sparse iterate $x^{t+1}$. Consequently, unlike other descent-based algorithms, $\mathcal{X}^t$ is not confined to the neighborhood of $k$-sparse vectors. Its
**Algorithm 1** Support Exploration Algorithm

1: **Input:** noisy observation $y$, sampling matrix $A$, sparsity $k$, step size $\eta$
2: **Output:** sparse vector $x$
3: Initialize $\mathcal{X}^0$
4: $t \leftarrow 0$
5: **repeat**
6:   $S^t \leftarrow \text{largest}_k(\mathcal{X}^t)$
7:   $\left\{ \begin{array}{ll} x_i^t & \leftarrow 0 \text{ for } i \in \overline{S^t} \\ x_{S^t}^t & \leftarrow A_{S^t}^\dagger y \end{array} \right.$
8:   $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - \eta A^T (Ax^t - y)$
9:   $t \leftarrow t + 1$
10:  **until** halting criterion is true
11:  $t_{\text{BEST}} \leftarrow \underset{t' \in [0, t]}{\operatorname{argmin}} \|Ax^{t'} - y\|_2^2$
12:  **return** $x^{t_{\text{BEST}}}$

evolution is not intended to make the objective function decrease at each iteration. In this regard, since the algorithm explores supports in a way that allows the functional to sometimes increase, the retained solution is the best one encountered along the iterations (line 11). Illustrations of this phenomenon are given in Appendix E.2 where one can see that the behavior of the loss along the iterations shows important variations when a new support is explored. This is an important difference with the aforementioned state-of-the-art algorithms, resulting in increased exploration.

An important feature of SEA is that it can be used as a post-processing of the solution $\hat{x}$ of another algorithm. This is simply done by initializing $\mathcal{X}^0 = \hat{x}$. In this case $S^0 = \text{supp}(\hat{x})$ (line 6) and $x^0$ improves or is equal to $\hat{x}$ (line 7). In the experiments, we have investigated the initialization with the result of OMP [39], ELS [2, 41] and the initialization $\mathcal{X}^0 = 0$. We observe that the initialization with ELS is generally preferable except for difficult problems, when columns of $A$ are very coherent (see Section 5.2).

Finally, as often, there are many possible strategies to design the halting criterion of the 'repeat' loop of Algorithm 1. It is clear that a more permissive criterion allows for more exploration and better results, at the expense of computation time. We have not investigated this aspect in the experiments and leave this study for the future. We preferred to focus our experiments on the illustration of the potential benefits of SEA and, as a consequence, we always used a large fixed number of passes in the 'repeat' loop of Algorithm 1.

Similarly, since we have $\mathcal{X}^t = \mathcal{X}^0 - \eta \sum_{t'=0}^{t-1} A^T (Ax^{t'} - y)$ for all $t \ge 1$, $\eta$ has no impact on $S^t$ and $x^t$ when $\mathcal{X}^0 = 0$ and therefore on the output $x^{t_{\text{BEST}}}$ of Algorithm 1 as largest$_k(\mathcal{X}^t)$ does not depend on $\eta$. In this case, indeed, the whole trajectory $(\mathcal{X}^t)_{t \in \mathbb{N}}$ is dilated by $\eta > 0$ and the dilation has no effect on the selected supports $S^t$. When $\mathcal{X}^0 \ne 0$, the initial support exploration variable is forgotten as the iterations progress. It is forgotten more rapidly when $\eta$ increases. We have not studied the tuning of this parameter in depth, leaving it for future research.

## 3.3 Computational Complexity

An efficient implementation of SEA is described in Appendix B.2. The analysis of the computational complexity of SEA is based on two facts. First, if the support $S^t$ obtained at line 6 has already been explored, then the sparse vector $x^t$ and the gradient $\eta A^T (Ax^t - y)$ have already been computed. So, if these quantities have been memorized (as in Algorithm 5, Appendix B.2), the cost of the iteration is negligible. The overall cost thus depends on the number of explored supports rather than on the number of iterations. Second, each time a new support is extracted, the cost of the iteration is dominated by solving the (unconstrained) linear system $A_{S^t}^T A_{S^t} x_{S^t} = A_{S^t}^T y$. While the pseudo-inverse is a convenient notation at line 7, the solution may be obtained more efficiently, e.g. in $\mathcal{O}(k^2n)$ to compute $A_{S^t}^T A_{S^t}$ and $A_{S^t}^T y$, and apply the conjugate gradient algorithm. The overall complexity is thus in $\mathcal{O}(n_{\text{supp}}k^2n)$ where $n_{\text{supp}}$ is the number of supports actually explored.

The complexity of HTP, OMP, OMPR and ELS is also dominated by the number of times $A_S^T A_S x_S = A_S^T y$ is solved for $S$ such that $|S| = k$. As for SEA, efficient implementations of HTP, OMPR and ELS can save computations by storing all the explored support and related iterates. The HTP, OMP and OMPR then depend on the number of explored supports in a similar way as SEA. The OMP solves $k$ instances of them which results in less exploration and less computational cost. The ELS is much more
**Algorithm 1** Support Exploration Algorithm

1: **Input:** noisy observation $y$, sampling matrix $A$, sparsity $k$, step size $\eta$
2: **Output:** sparse vector $x$
3: Initialize $\mathcal{X}^0$
4: $t \leftarrow 0$
5: **repeat**
6: &nbsp;&nbsp;&nbsp;&nbsp;$S^t \leftarrow \text{largest}_k(\mathcal{X}^t)$
7: &nbsp;&nbsp;&nbsp;&nbsp;$$ \begin{cases} x_i^t & \leftarrow 0 \text{ for } i \in \bar{S}^t \\ x_{S^t}^t & \leftarrow A_{S^t}^\dagger y \end{cases} $$
8: &nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - \eta A^T (Ax^t - y)$
9: &nbsp;&nbsp;&nbsp;&nbsp;$t \leftarrow t + 1$
10: **until** halting criterion is true
11: $t_{\text{BEST}} \leftarrow \underset{t' \in [0, t]}{\operatorname{argmin}} \|Ax^{t'} - y\|_2^2$
12: **return** $x^{t_{\text{BEST}}}$

evolution is not intended to make the objective function decrease at each iteration. In this regard, since the algorithm explores supports in a way that allows the functional to sometimes increase, the retained solution is the best one encountered along the iterations (line 11). Illustrations of this phenomenon are given in Appendix E.2 where one can see that the behavior of the loss along the iterations shows important variations when a new support is explored. This is an important difference with the aforementioned state-of-the-art algorithms, resulting in increased exploration.

An important feature of SEA is that it can be used as a post-processing of the solution $\hat{x}$ of another algorithm. This is simply done by initializing $\mathcal{X}^0 = \hat{x}$. In this case $S^0 = \text{supp}(\hat{x})$ (line 6) and $x^0$ improves or is equal to $\hat{x}$ (line 7). In the experiments, we have investigated the initialization with the result of OMP [39], ELS [2, 41] and the initialization $\mathcal{X}^0 = 0$. We observe that the initialization with ELS is generally preferable except for difficult problems, when columns of $A$ are very coherent (see Section 5.2).

Finally, as often, there are many possible strategies to design the halting criterion of the 'repeat' loop of Algorithm 1. It is clear that a more permissive criterion allows for more exploration and better results, at the expense of computation time. We have not investigated this aspect in the experiments and leave this study for the future. We preferred to focus our experiments on the illustration of the potential benefits of SEA and, as a consequence, we always used a large fixed number of passes in the 'repeat' loop of Algorithm 1.

Similarly, since we have $\mathcal{X}^t = \mathcal{X}^0 - \eta \sum_{t'=0}^{t-1} A^T (Ax^{t'} - y)$ for all $t \ge 1$, $\eta$ has no impact on $S^t$ and $x^t$ when $\mathcal{X}^0 = 0$ and therefore on the output $x^{t_{\text{BEST}}}$ of Algorithm 1 as largest$_k(\mathcal{X}^t)$ does not depend on $\eta$. In this case, indeed, the whole trajectory $(\mathcal{X}^t)_{t \in \mathbb{N}}$ is dilated by $\eta > 0$ and the dilation has no effect on the selected supports $S^t$. When $\mathcal{X}^0 \ne 0$, the initial support exploration variable is forgotten as the iterations progress. It is forgotten more rapidly when $\eta$ increases. We have not studied the tuning of this parameter in depth, leaving it for future research.

## 3.3 Computational Complexity

An efficient implementation of SEA is described in Appendix B.2. The analysis of the computational complexity of SEA is based on two facts. First, if the support $S^t$ obtained at line 6 has already been explored, then the sparse vector $x^t$ and the gradient $\eta A^T (Ax^t - y)$ have already been computed. So, if these quantities have been memorized (as in Algorithm 5, Appendix B.2), the cost of the iteration is negligible. The overall cost thus depends on the number of explored supports rather than on the number of iterations. Second, each time a new support is extracted, the cost of the iteration is dominated by solving the (unconstrained) linear system $A_{S^t}^T A_{S^t} x_{S^t} = A_{S^t}^T y$. While the pseudo-inverse is a convenient notation at line 7, the solution may be obtained more efficiently, e.g. in $\mathcal{O}(k^2n)$ to compute $A_{S^t}^T A_{S^t}$ and $A_{S^t}^T y$, and apply the conjugate gradient algorithm. The overall complexity is thus in $\mathcal{O}(n_{\text{supp}}k^2n)$ where $n_{\text{supp}}$ is the number of supports actually explored.

The complexity of HTP, OMP, OMPR and ELS is also dominated by the number of times $A_S^T A_S x_S = A_S^T y$ is solved for $S$ such that $|S| = k$. As for SEA, efficient implementations of HTP, OMPR and ELS can save computations by storing all the explored support and related iterates. The HTP, OMP and OMPR then depend on the number of explored supports in a similar way as SEA. The OMP solves $k$ instances of them which results in less exploration and less computational cost. The ELS is much more
demanding since it explores $(n-k)$ supports at each iteration, many of which are irrelevant. IHT has a lower complexity than SEA since it never inverses the system $A_S^T A_S x_S = A_S^T y$.

As we will see in the deconvolution experiments in Section 5.2, SEA outperforms ELS (see Figure 4) while exploring two times less supports (see Appendices E.2 and E.3). The possibility of performing a random search with the same computational cost as SEA has been studied in Appendix E.11.

# 4 Theoretical Analysis

In this section, we provide the theorem stating that SEA recovers the correct support for some$^4 x^*$ when the matrix $A$ satisfies a RIP constraint. Then, we compare the conditions with existing support recovery conditions for state-of-the-art algorithms. In addition to the statements in this section, recovery statements are given in Appendices C. The interest of the theorems lies mainly in the fact that they apply to an instance of STE, for which guarantees are rare. From the practitioner's point of view, the theoretical analysis is not useful since SEA mostly shows promises in coherent scenarios in which the RIP hypothesis is not satisfied.

In this section, we assume that columns of $A$ are normalized: for any $i \in [1, n]$, $\|A_i\|_2 = 1$. As has been standard practice since Candès and Tao first proposed it in [12], we define for all $l \in [1, n]$ the $l$th Restricted Isometry Constant of $A$ as the smallest non-negative number $\delta_l$ such that for any $x \in \mathbb{R}^n$, such that $\|x\|_0 \le l$,

$$ (1 - \delta_l) \|x\|_2^2 \le \|Ax\|_2^2 \le (1 + \delta_l) \|x\|_2^2. \quad (4) $$

If $\delta_l < 1$, $A$ is said to satisfy the Restricted Isometry Property of order $l$ or the $l$-RIP.

In this section, we assume that $A$ satisfies the $(2k+1)$-RIP. In the scenarios of interest, $\delta_{2k+1}$ is small. We define

$$ \alpha_k^{\text{RIP}} = \delta_{2k+1} \left( \frac{\delta_{2k}}{1 - \delta_k} + 1 \right) \in \mathbb{R}_+^* \quad \text{and} \quad \gamma_k^{\text{RIP}} = \delta_{2k+1} \frac{\sqrt{1 + \delta_k}}{1 - \delta_k} + 1 \in \mathbb{R}_+^*. \quad (5) $$

As soon as $\delta_k$ is far from 1 (for example $\delta_k \le \frac{1}{2}$), $\alpha_k^{\text{RIP}}$ has the order of magnitude of $\delta_{2k+1}$ (in the example $\delta_{2k+1} \le \alpha_k^{\text{RIP}} \le 3\delta_{2k+1}$) and $\gamma_k^{\text{RIP}}$ has the order of magnitude of $1 + \delta_{2k+1}$ (in the example $(1 + \delta_{2k+1}) \le \gamma_k^{\text{RIP}} \le \sqrt{6}(1 + \delta_{2k+1}))$.

As is typical of support recovery statements, the next theorem includes a condition on $x^*$. We call this condition the Recovery Condition for the RIP case (RC$_{RIP}$). It is defined by

$$ \gamma_k^{\text{RIP}} \|e\|_2 < \frac{\min_{i \in S^*} |x_i^*|}{2k} - \alpha_k^{\text{RIP}} \|x^*\|_2. \quad (RC_{\text{RIP}}) $$

**Theorem 4.1 (Recovery - RIP case).** Assume $A$ satisfies the $(2k+1)$-RIP and$^5$ for all $i \in [1, n]$, $\|A_i\|_2 = 1$. Assume moreover that $x^*$ satisfies (RC$_{RIP}$).

Then for all initializations $\mathcal{X}^0$ and all $\eta > 0$, there exists $t_s \le T_{\text{RIP}}$ such that $S^* \subseteq S^{t_s}$, where

$$ T_{\text{RIP}} = \frac{2k \frac{\|\mathcal{X}^0\|_{\infty}}{\eta} + (k+1) \min_{i \in S^*} |x_i^*|}{\min_{i \in S^*} |x_i^*| - 2k (\alpha_k^{\text{RIP}} \|x^*\|_2 + \gamma_k^{\text{RIP}} \|e\|_2)}. \quad (6) $$

If moreover, $x^*$ is such that

$$ \min_{i \in S^*} |x_i^*| > \frac{2}{\sqrt{1 - \delta_{2k}}} \|e\|_2 \quad (7) $$

and SEA performs more than $T_{\text{RIP}}$ iterations, then $S^* \subseteq S^{t_{\text{BEST}}}$ and $\|x^{t_{\text{BEST}}}-x^*\|_2 \le \frac{2}{\sqrt{1-\delta_k}} \|e\|_2$.

The proof is in Appendix C. To introduce the proof and provide the main intuition, we first detail in Appendix C.1 a theorem and its proof stating that a variant of the support exploration algorithm, see Algorithm 6, using the oracle update rule defined by $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - u^t$ where for all $i \in [1, n]$

$$ u_i^t = \begin{cases} -\eta x_i^* & i \in S^* \cap \overline{S}^t \\ 0 & i \in \overline{S}^* \cup S^t, \end{cases} $$

$$ (1 - \delta_l) \|\Delta^{-1} x\|_2^2 \le \|A \Delta^{-1} x\|_2^2 \le (1 + \delta_l) \|\Delta^{-1} x\|_2^2. $$

Using $1 - \delta_1 \le \|A_i\|_2^2 \le 1 + \delta_1$, we can derive $l$-RIP constants for the normalized matrix $A\Delta^{-1}$.

<footnote>$^4$In particular it is necessary that $\min_{i \in S^*} |x_i^*|$ is sufficiently large.</footnote>

<footnote>$^5$The normalization aims at simplifying formulas by guaranteeing that $\delta_1 = 0$. It is done at no expense since, if $A$ is not normalized but satisfies (4) for $l > 1$, its normalization only has a small impact on $\delta_l$. Indeed, considering $\Delta \in \mathbb{R}^{n \times n}$ diagonal such that $\Delta_{i,i} = \|A_i\|_2$, $A\Delta^{-1}$ is normalized and for all $l$-sparse vector $x$</footnote>
demanding since it explores $(n-k)$ supports at each iteration, many of which are irrelevant. IHT has a lower complexity than SEA since it never inverses the system $A_S^T A_S x_S = A_S^T y$.

As we will see in the deconvolution experiments in Section 5.2, SEA outperforms ELS (see Figure 4) while exploring two times less supports (see Appendices E.2 and E.3). The possibility of performing a random search with the same computational cost as SEA has been studied in Appendix E.11.

# 4 Theoretical Analysis

In this section, we provide the theorem stating that SEA recovers the correct support for some<fnref n="4" /> $x^*$ when the matrix $A$ satisfies a RIP constraint. Then, we compare the conditions with existing support recovery conditions for state-of-the-art algorithms. In addition to the statements in this section, recovery statements are given in Appendices C. The interest of the theorems lies mainly in the fact that they apply to an instance of STE, for which guarantees are rare. From the practitioner's point of view, the theoretical analysis is not useful since SEA mostly shows promises in coherent scenarios in which the RIP hypothesis is not satisfied.

In this section, we assume that columns of $A$ are normalized: for any $i \in [1, n]$, $\|A_i\|_2 = 1$. As has been standard practice since Candès and Tao first proposed it in [12], we define for all $l \in [1, n]$ the $l$th Restricted Isometry Constant of $A$ as the smallest non-negative number $\delta_l$ such that for any $x \in \mathbb{R}^n$, such that $\|x\|_0 \le l$,

$$ (1 - \delta_l) \|x\|_2^2 \le \|Ax\|_2^2 \le (1 + \delta_l) \|x\|_2^2. \quad (4) $$

If $\delta_l < 1$, $A$ is said to satisfy the Restricted Isometry Property of order $l$ or the $l$-RIP.

In this section, we assume that $A$ satisfies the $(2k+1)$-RIP. In the scenarios of interest, $\delta_{2k+1}$ is small. We define

$$ \alpha_k^{\text{RIP}} = \delta_{2k+1} \left( \frac{\delta_{2k}}{1 - \delta_k} + 1 \right) \in \mathbb{R}_+^* \quad \text{and} \quad \gamma_k^{\text{RIP}} = \delta_{2k+1} \frac{\sqrt{1 + \delta_k}}{1 - \delta_k} + 1 \in \mathbb{R}_+^*. \quad (5) $$

As soon $\delta_k$ is far from 1 (for example $\delta_k \le \frac{1}{2}$), $\alpha_k^{\text{RIP}}$ has the order of magnitude of $\delta_{2k+1}$ (in the example $\delta_{2k+1} \le \alpha_k^{\text{RIP}} \le 3\delta_{2k+1}$) and $\gamma_k^{\text{RIP}}$ has the order of magnitude of $1 + \delta_{2k+1}$ (in the example $(1 + \delta_{2k+1}) \le \gamma_k^{\text{RIP}} \le \sqrt{6}(1 + \delta_{2k+1}))$.

As is typical of support recovery statements, the next theorem includes a condition on $x^*$. We call this condition the Recovery Condition for the RIP case (RC$_{RIP}$). It is defined by

$$ \gamma_k^{\text{RIP}} \|e\|_2 < \frac{\min_{i \in S^*} |x_i^*|}{2k} - \alpha_k^{\text{RIP}} \|x^*\|_2. \quad (\text{RC}_{\text{RIP}}) $$

**Theorem 4.1 (Recovery - RIP case).** Assume $A$ satisfies the $(2k+1)$-RIP and<fnref n="5" /> for all $i \in [1, n]$, $\|A_i\|_2 = 1$. Assume moreover that $x^*$ satisfies (RC$_{RIP}$).

Then for all initializations $\mathcal{X}^0$ and all $\eta > 0$, there exists $t_s \le T_{\text{RIP}}$ such that $S^* \subseteq S^{t_s}$, where

$$ T_{\text{RIP}} = \frac{2k \frac{\|\mathcal{X}^0\|_{\infty}}{\eta} + (k+1) \min_{i \in S^*} |x_i^*|}{\min_{i \in S^*} |x_i^*| - 2k (\alpha_k^{\text{RIP}} \|x^*\|_2 + \gamma_k^{\text{RIP}} \|e\|_2)}. \quad (6) $$

If moreover, $x^*$ is such that

$$ \min_{i \in S^*} |x_i^*| > \frac{2}{\sqrt{1 - \delta_{2k}}} \|e\|_2 \quad (7) $$

and SEA performs more than $T_{\text{RIP}}$ iterations, then $S^* \subseteq S^{t_{\text{BEST}}}$ and $\|x^{t_{\text{BEST}}}-x^*\|_2 \le \frac{2}{\sqrt{1-\delta_k}} \|e\|_2$.

The proof is in Appendix C. To introduce the proof and provide the main intuition, we first detail in Appendix C.1 a theorem and its proof stating that a variant of the support exploration algorithm, see Algorithm 6, using the oracle update rule defined by $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - u^t$ where for all $i \in [1, n]$

$$ u_i^t = \begin{cases} -\eta x_i^* & i \in S^* \cap \overline{S}^t \\ 0 & i \in \overline{S}^* \cup S^t, \end{cases} $$

$$ (1 - \delta_l) \|\Delta^{-1} x\|_2^2 \le \|A \Delta^{-1} x\|_2^2 \le (1 + \delta_l) \|\Delta^{-1} x\|_2^2. $$

Using $1 - \delta_1 \le \|A_i\|_2^2 \le 1 + \delta_1$, we can derive $l$-RIP constants for the normalized matrix $A\Delta^{-1}$.

<footnote n="4">⁴In particular it is necessary that $\min_{i \in S^*} |x_i^*|$ is sufficiently large.</footnote>

<footnote n="5">⁵The normalization aims at simplifying formulas by guaranteeing that $\delta_1 = 0$. It is done at no expense since, if $A$ is not normalized but satisfies (4) for $l > 1$, its normalization only has a small impact on $\delta_l$. Indeed, considering $\Delta \in \mathbb{R}^{n \times n}$ diagonal such that $\Delta_{i,i} = \|A_i\|_2$, $A\Delta^{-1}$ is normalized and for all $l$-sparse vector $x$</footnote>
always recovers the true support $S^*$. This update rule depends on $x^*$ and has no practical application. Its interest lies in providing an 'ideal' update rule that allows for the fast recovery of the true support, as seen in Theorem C.1.

The intuition behind the success of the oracle update rule is that the non-zero entries of $u_i^t$ are for indices $i$ from the true support $S^*$ but for which $|\mathcal{X}_i^t|$ is too small to be selected in $S^t$ at line 6. Whatever the initial content of $\mathcal{X}^0$, the oracle update rule always adds the same increment to $\mathcal{X}_i^t$, for $i \in S^* \cap \bar{S}^t$, and those for $i \in \bar{S}^*$ never change. This guarantees that, at some subsequent iteration $t' \ge t$, the true support $S^*$ is recovered among the $k$ largest absolute entries in $\mathcal{X}^{t'}$, i.e., $S^* \subseteq S^{t'}$.

The proof of Theorem 4.1 relies on measuring the deviation of the trajectory $(\mathcal{X}^t)_{t \in \mathbb{N}}$ defined by Algorithm 1 from the trajectory defined with the oracle update rule. More precisely, in Theorem C.4 of Appendix C.2, we provide a sufficient condition on the discrepancy between the oracle update and the STE-update guaranteeing that SEA visits the true support $S^*$. Then, in Appendix C.4, we establish that the hypotheses of Theorem 4.1 ensure that the discrepancy is sufficiently small to satisfy the hypothesis of Theorem C.4.

We also establish in Corollary C.7 of Appendix C.3 that when the columns of $A$ are orthonormal and $e = 0$, SEA recovers $x^*$ in less than $k+1$ iterations. Despite being a sanity check with no interest in applications, this result provides a meaningful case where the oracle update rule and the STE update rule coincide.

We emphasize that none of the proofs rely on the fact that a function decays. In particular, as will be illustrated in the experiments, $(F(H(\mathcal{X}^t)))_{t \in \mathbb{N}}$ generally exhibits erratic behavior. This is because, by construction in Algorithm 1, the next support, designated by $\mathcal{X}^{t+1}$, is not restricted to supports for which $F \circ H$ decays. SEA explores more supports than algorithms with this restriction and, for instance, does not get trapped in local minima.

Let us now discuss $T_{RIP}$.

When $(RC_{RIP})$ holds, $T_{RIP}$ increases as $\min_{i \in S^*} |x_i^*| - 2k (\alpha_k^{RIP} \|x^*\|_2 + \gamma_k^{RIP} \|e\|_2)$ decreases. In particular, the number of iterations required by the algorithm to provide the correct solution increases when the information on some of the columns of $S^*$ diminishes, i.e. when $\min_{i \in S^*} |x_i^*|$ decreases. Also, the initializations $\mathcal{X}^0 \neq 0$ have an apparent negative impact on the number of iterations required in the worst case. This is because in the worst-case $\mathcal{X}^0$ is poorly chosen and SEA needs iterations to correct this poor choice.

When the conditions of Theorem 4.1 are met, any $\mathcal{X}^0$ and $\eta$ permit the recovery of $S^*$. $\mathcal{X}^0$ and $\eta$ only influence $T_{RIP}$. In this regard, since the larger $\eta$, the faster SEA overrides the initialization $\mathcal{X}^0$, the choice of $\eta$ is very much related to the question of the quality of the initialization. The latter is often beneficial in practice.

To illustrate $(RC_{RIP})$, we provide below a simplified condition which is shown in Corollary 4.2 to be stronger than $(RC_{RIP})$ in the noiseless scenario. We say $x^*$ satisfies the Simplified Recovery Condition in the RIP case if there exists $\Lambda \in (0, 1)$ such that

$$2k \alpha_k^{RIP} \frac{\|x^*\|_2}{\min_{i \in S^*} |x_i^*|} \le \Lambda. \qquad (RC_{SRIP})$$

**Corollary 4.2 (Noiseless recovery - simplified RIP case).** Assume $\|e\|_2 = 0$, $A$ satisfies the $(2k+1)$-RIP and for all $i \in [1, n]$, $\|A_i\|_2 = 1$.

If moreover $x^*$ satisfies $(RC_{SRIP})$, then $x^*$ satisfies $(RC_{RIP})$. As a consequence, for $\mathcal{X}^0 = 0$ and for all $\eta > 0$, if SEA performs more than $T_{SRIP} = \frac{k+1}{1-\Lambda}$ iterations, we have $S^* \subseteq S^{\text{BEST}}$ and $x^{\text{BEST}} = x^*$.

The proof is in Appendix C.5.

Compared to the support recovery guarantees for the LASSO [45, 37, 49], the OMP [10], the HTP [25, 48] and the ARHT [2] the recovery conditions provided in Theorem 4.1 and Corollary 4.2 for SEA are stronger. All conditions involve a condition on the incoherence of $A$ and a condition similar to (7). In the case of the LASSO algorithm, the latter is not very explicit. However, none of the support recovery conditions involve a condition like $(RC_{RIP})$ and $(RC_{SRIP})$. Let elaborate on these two conditions.

One notable limitation of $(RC_{RIP})$ and $(RC_{SRIP})$ is that if $\alpha_k^{RIP} \neq 0$, there is no guarantee of recovering the support of an $x^*$ such that $\max_{i \in S^*} |x_i^*| \gg \min_{i \in S^*} |x_i^*|$. In fact, by increasing $x_j^*$ for $j \in S^*$ such that $|x_j^*| \neq \min_{i \in S^*} |x_i^*|$, we can maintain $\min_{i \in S^*} |x_i^*|$ unchanged while sufficiently increasing $\|x^*\|_2$, causing $(RC_{RIP})$ and $(RC_{SRIP})$ to fail. We do not have this problem with conditions similar to (7). We experimented on the influence of $\frac{\|x^*\|_2}{\min_{i \in S^*} |x_i^*|}$ on the recovery performances in Appendix E.9.

Another critical question arises: What (necessary and sufficient) condition on $\alpha_k^{RIP}$, $\gamma_k^{RIP}$ and $\|e\|_2$ ensures the existence of $x^*$ satisfying $(RC_{RIP})$ or $(RC_{SRIP})$?
always recovers the true support $S^*$. This update rule depends on $x^*$ and has no practical application. Its interest lies in providing an 'ideal' update rule that allows for the fast recovery of the true support, as seen in Theorem C.1.

The intuition behind the success of the oracle update rule is that the non-zero entries of $u_i^t$ are for indices $i$ from the true support $S^*$ but for which $|\mathcal{X}_i^t|$ is too small to be selected in $S^t$ at line 6. Whatever the initial content of $\mathcal{X}^0$, the oracle update rule always adds the same increment to $\mathcal{X}_i^t$, for $i \in S^* \cap \bar{S}^t$, and those for $i \in \bar{S}^*$ never change. This guarantees that, at some subsequent iteration $t' \ge t$, the true support $S^*$ is recovered among the $k$ largest absolute entries in $\mathcal{X}^{t'}$, i.e., $S^* \subseteq S^{t'}$.

The proof of Theorem 4.1 relies on measuring the deviation of the trajectory $(\mathcal{X}^t)_{t \in \mathbb{N}}$ defined by Algorithm 1 from the trajectory defined with the oracle update rule. More precisely, in Theorem C.4 of Appendix C.2, we provide a sufficient condition on the discrepancy between the oracle update and the STE-update guaranteeing that SEA visits the true support $S^*$. Then, in Appendix C.4, we establish that the hypotheses of Theorem 4.1 ensure that the discrepancy is sufficiently small to satisfy the hypothesis of Theorem C.4.

We also establish in Corollary C.7 of Appendix C.3 that when the columns of $A$ are orthonormal and $e = 0$, SEA recovers $x^*$ in less than $k + 1$ iterations. Despite being a sanity check with no interest in applications, this result provides a meaningful case where the oracle update rule and the STE update rule coincide.

We emphasize that none of the proofs rely on the fact that a function decays. In particular, as will be illustrated in the experiments, $(F(H(\mathcal{X}^t)))_{t \in \mathbb{N}}$ generally exhibits erratic behavior. This is because, by construction in Algorithm 1, the next support, designated by $\mathcal{X}^{t+1}$, is not restricted to supports for which $F \circ H$ decays. SEA explores more supports than algorithms with this restriction and, for instance, does not get trapped in local minima.

Let us now discuss $T_{RIP}$.

When $(RC_{RIP})$ holds, $T_{RIP}$ increases as $\min_{i \in S^*} |x_i^*| - 2k (\alpha_k^{RIP} \|x^*\|_2 + \gamma_k^{RIP} \|e\|_2)$ decreases. In particular, the number of iterations required by the algorithm to provide the correct solution increases when the information on some of the columns of $S^*$ diminishes, i.e. when $\min_{i \in S^*} |x_i^*|$ decreases. Also, the initializations $\mathcal{X}^0 \neq 0$ have an apparent negative impact on the number of iterations required in the worst case. This is because in the worst-case $\mathcal{X}^0$ is poorly chosen and SEA needs iterations to correct this poor choice.

When the conditions of Theorem 4.1 are met, any $\mathcal{X}^0$ and $\eta$ permit the recovery of $S^*$. $\mathcal{X}^0$ and $\eta$ only influence $T_{RIP}$. In this regard, since the larger $\eta$, the faster SEA overrides the initialization $\mathcal{X}^0$, the choice of $\eta$ is very much related to the question of the quality of the initialization. The latter is often beneficial in practice.

To illustrate $(RC_{RIP})$, we provide below a simplified condition which is shown in Corollary 4.2 to be stronger than $(RC_{RIP})$ in the noiseless scenario. We say $x^*$ satisfies the Simplified Recovery Condition in the RIP case if there exists $\Lambda \in (0, 1)$ such that

$$2k \alpha_k^{RIP} \frac{\|x^*\|_2}{\min_{i \in S^*} |x_i^*|} \le \Lambda. \qquad (RC_{SRIP})$$

**Corollary 4.2 (Noiseless recovery - simplified RIP case).** Assume $\|e\|_2 = 0$, $A$ satisfies the $(2k+1)$-RIP and for all $i \in [1, n]$, $\|A_i\|_2 = 1$.

If moreover $x^*$ satisfies $(RC_{SRIP})$, then $x^*$ satisfies $(RC_{RIP})$. As a consequence, for $\mathcal{X}^0 = 0$ and for all $\eta > 0$, if SEA performs more than $T_{SRIP} = \frac{k+1}{1-\Lambda}$ iterations, we have $S^* \subseteq S^{\text{BEST}}$ and $x^{\text{BEST}} = x^*$.

The proof is in Appendix C.5.

Compared to the support recovery guarantees for the LASSO [45, 37, 49], the OMP [10], the HTP [25, 48] and the ARHT [2] the recovery conditions provided in Theorem 4.1 and Corollary 4.2 for SEA are stronger. All conditions involve a condition on the incoherence of $A$ and a condition similar to (7). In the case of the LASSO algorithm, the latter is not very explicit. However, none of the support recovery conditions involve a condition like $(RC_{RIP})$ and $(RC_{SRIP})$. Let elaborate on these two conditions.

One notable limitation of $(RC_{RIP})$ and $(RC_{SRIP})$ is that if $\alpha_k^{RIP} \neq 0$, there is no guarantee of recovering the support of an $x^*$ such that $\max_{i \in S^*} |x_i^*| \gg \min_{i \in S^*} |x_i^*|$. In fact, by increasing $x_j^*$ for $j \in S^*$ such that $|x_j^*| \neq \min_{i \in S^*} |x_i^*|$, we can maintain $\min_{i \in S^*} |x_i^*|$ unchanged while sufficiently increasing $\|x^*\|_2$, causing $(RC_{RIP})$ and $(RC_{SRIP})$ to fail. We do not have this problem with conditions similar to (7). We experimented on the influence of $\frac{\|x^*\|_2}{\min_{i \in S^*} |x_i^*|}$ on the recovery performances in Appendix E.9.

Another critical question arises: What (necessary and sufficient) condition on $\alpha_k^{RIP}$, $\gamma_k^{RIP}$ and $\|e\|_2$ ensures the existence of $x^*$ satisfying $(RC_{RIP})$ or $(RC_{SRIP})$?
To answer this question for the condition ($RC_{RIP}$), we first remark that there exists $x^*$ satisfying ($RC_{RIP}$) if and only if there exists some constant $c \in \mathbb{R}$ and $x^*$, such that for all $i \in S^*$, $x_i^* = c$, satisfying ($RC_{RIP}$). Indeed, if $x^*$ satisfies ($RC_{RIP}$), we take $c = \min_{i \in S^*} |x_i^*|$ and the vector $\tilde{x}^*$ defined by $\tilde{x}_i^* = \min_{i \in S^*} |x_i^*|$ for all $i \in S^*$, and $\tilde{x}_i^* = 0$ for all $i \notin S^*$ also satisfies ($RC_{RIP}$), since we have $\min_{i \in S^*} |\tilde{x}_i^*| = \min_{i \in S^*} |x_i^*|$ and $\|\tilde{x}^*\|_2 \le \|x^*\|_2$.

We can rewrite ($RC_{RIP}$) when for all $i \in S^*$, $x_i^* = c$ for some constant $c \in \mathbb{R}$ and obtain

$$ \gamma_k^{RIP} \|e\|_2 < |c| \left( \frac{1 - 2k\alpha_k^{RIP} |S^*|^{\frac{1}{2}}}{2k} \right). $$

The existence of $c \in \mathbb{R}$ such that this condition holds only depends on the sign of $1 - 2k\alpha_k^{RIP} |S^*|^{\frac{1}{2}}$. If $1 - 2k\alpha_k^{RIP} |S^*|^{\frac{1}{2}} \le 0$, there does not exist any $c$ satisfying the condition; if $1 - 2k\alpha_k^{RIP} |S^*|^{\frac{1}{2}} > 0$, any $c \in \mathbb{R}$ satisfying

$$ |c| \ge \left( \frac{2k\gamma_k^{RIP}}{1 - 2k\alpha_k^{RIP} |S^*|^{\frac{1}{2}}} \right) \|e\|_2 $$

leads to an $x^*$ that satisfies ($RC_{RIP}$). Therefore, when $|S^*| = k$, the condition $\alpha_k^{RIP} < \frac{1}{2}k^{-\frac{3}{2}}$ is necessary and sufficient to guarantee the existence of an $x^*$ satisfying ($RC_{RIP}$). Similar developments concerning ($RC_{SRIIP}$) are provided in Appendix C.6.

We remind that $\alpha_k^{RIP}$ has the order of magnitude of $\delta_{2k+1}$. Therefore, the condition is $\alpha_k^{RIP} < \frac{1}{2}k^{-\frac{3}{2}}$ is more stringent than the equivalent conditions for other state-of-the-art methods. For instance, the conditions described in [2] are of the form $\delta_{2k} < C$, for a universal constant $C < 1$.

Initializing SEA with the solution of an algorithm enjoying better conditions of recovery is a simple and effective way for SEA to inherit its support recovery guarantee as soon as (7) holds. This can formally be proved using the same proof as in Appendix C.4.3.

As will be seen later in the experiments of Section 5, SEA performs well even when $A$ is coherent. This is not explained by Theorem 4.1 and Corollary 4.2 which use the RIP assumption. The main interest of the above theoretical results lies in the fact that they apply to an instance of the STE. Another theory needs to be developed to explain the good behavior of SEA for sparse support recovery when $A$ is coherent.

# 5 Experimental Analysis

We compare SEA to state-of-the-art algorithms on two tasks in the noisy setting: phase transition diagrams (Section 5.1 and Appendix D) and spike deconvolution problems for signal processing (Section 5.2 and Appendix E). For completeness, additional comparisons between SEA and state-of-the-art algorithms for linear and logistic regression tasks in supervised learning settings are provided in Appendix F.

The tested algorithms are Iterative Hard Thresholding (IHT) [8], Hard Thresholding Pursuit (HTP) [25], Orthogonal Matching Pursuit [36, 39], OMP with Replacement (OMPR) [32] and Exhaustive Local Search (ELS) [2]. OMPR and ELS are initialized with the solution of OMP. Three versions of SEA are studied: the cold-start version $SEA_0$, where SEA is initialized with the null vector, and the warm-start versions $SEA_{ELS}$ and $SEA_{OMP}$, where SEA is initialized with the solutions of ELS and OMP, respectively. We have also studied HTP and IHT initialized with OMP and ELS. They are called $HTP_{OMP}$, $HTP_{ELS}$, $IHT_{OMP}$ and $IHT_{ELS}$.

For all algorithms, each least-square projection for a fixed support, as in Line 7 of Algorithm 1, is solved using the conjugate gradient descent of SciPy [44]. For all algorithms, 256k iterations are performed. The results of HTP and to a lesser extent IHT and SEA depend on the choice of the step size. For the sake of fairness of the comparison with OMP, OMPR and ELS, we did not optimize the choice of the step size. The step size of SEA, HTP, and IHT is arbitrarily<fnref n="6" /> fixed to $\eta = \frac{1.8}{L}$, where $L$ is the spectral radius of $A$. The columns of $A$ are normalized before solving the problem. The sparse vector $x^* \in \mathbb{R}^n$ is random. Indexes of the support are randomly picked, uniformly without replacement. The non-zero entries of $x^*$ are drawn uniformly in $[-2, -1] \cup [1, 2]$ as in [24]. The noise $e$ is drawn uniformly using the same method as described in [7]. Their detailed descriptions are in the next two sections. For each experiment, the metrics used for performance evaluation are defined in the corresponding subsection. The code is implemented in Python 3 and is available in the git repository of the project <fnref n="7" />. As explained

<footnote n="6">⁶We do not report further experiments for $\eta \in \{\frac{2l}{L} | \text{ for } l = [-3, +3]\}$ that do not significantly alter the results in terms of running time, stability, performance, and do not impact our conclusions on phase transition diagrams. Variation of the step size for IHT and HTP in the deconvolution experiment is reported in Appendix E.10 and does not impact our conclusions.</footnote>

<footnote n="7">⁷https://gitlab.lis-lab.fr/valentin.emiya/sea-icml-2024</footnote>
To answer this question for the condition ($RC_{RIP}$), we first remark that there exists $x^*$ satisfying ($RC_{RIP}$) if and only if there exists some constant $c \in \mathbb{R}$ and $x^*$, such that for all $i \in S^*$, $x_i^* = c$, satisfying ($RC_{RIP}$). Indeed, if $x^*$ satisfies ($RC_{RIP}$), we take $c = \min_{i \in S^*} |x_i^*|$ and the vector $\tilde{x}^*$ defined by $\tilde{x}_i^* = \min_{i \in S^*} |x_i^*|$ for all $i \in S^*$, and $\tilde{x}_i^* = 0$ for all $i \notin S^*$ also satisfies ($RC_{RIP}$), since we have $\min_{i \in S^*} |\tilde{x}_i^*| = \min_{i \in S^*} |x_i^*|$ and $\|\tilde{x}^*\|_2 \le \|x^*\|_2$.

We can rewrite ($RC_{RIP}$) when for all $i \in S^*$, $x_i^* = c$ for some constant $c \in \mathbb{R}$ and obtain

$$ \gamma_k^{RIP} \|e\|_2 < |c| \left( \frac{1 - 2k\alpha_k^{RIP} |S^*|^{\frac{1}{2}}}{2k} \right). $$

The existence of $c \in \mathbb{R}$ such that this condition holds only depends on the sign of $1 - 2k\alpha_k^{RIP} |S^*|^{\frac{1}{2}}$. If $1 - 2k\alpha_k^{RIP} |S^*|^{\frac{1}{2}} \le 0$, there does not exist any $c$ satisfying the condition; if $1 - 2k\alpha_k^{RIP} |S^*|^{\frac{1}{2}} > 0$, any $c \in \mathbb{R}$ satisfying

$$ |c| \ge \left( \frac{2k\gamma_k^{RIP}}{1 - 2k\alpha_k^{RIP} |S^*|^{\frac{1}{2}}} \right) \|e\|_2 $$

leads to an $x^*$ that satisfies ($RC_{RIP}$). Therefore, when $|S^*| = k$, the condition $\alpha_k^{RIP} < \frac{1}{2}k^{-\frac{3}{2}}$ is necessary and sufficient to guarantee the existence of an $x^*$ satisfying ($RC_{RIP}$). Similar developments concerning ($RC_{S RIP}$) are provided in Appendix C.6.

We remind that $\alpha_k^{RIP}$ has the order of magnitude of $\delta_{2k+1}$. Therefore, the condition is $\alpha_k^{RIP} < \frac{1}{2}k^{-\frac{3}{2}}$ is more stringent than the equivalent conditions for other state-of-the-art methods. For instance, the conditions described in [2] are of the form $\delta_{2k} < C$, for a universal constant $C < 1$.

Initializing SEA with the solution of an algorithm enjoying better conditions of recovery is a simple and effective way for SEA to inherit its support recovery guarantee as soon as (7) holds. This can formally be proved using the same proof as in Appendix C.4.3.

As will be seen later in the experiments of Section 5, SEA performs well even when $A$ is coherent. This is not explained by Theorem 4.1 and Corollary 4.2 which use the RIP assumption. The main interest of the above theoretical results lies in the fact that they apply to an instance of the STE. Another theory needs to be developed to explain the good behavior of SEA for sparse support recovery when $A$ is coherent.

# 5 Experimental Analysis

We compare SEA to state-of-the-art algorithms on two tasks in the noisy setting: phase transition diagrams (Section 5.1 and Appendix D) and spike deconvolution problems for signal processing (Section 5.2 and Appendix E). For completeness, additional comparisons between SEA and state-of-the-art algorithms for linear and logistic regression tasks in supervised learning settings are provided in Appendix F.

The tested algorithms are Iterative Hard Thresholding (IHT) [8], Hard Thresholding Pursuit (HTP) [25], Orthogonal Matching Pursuit [36, 39], OMP with Replacement (OMPR) [32] and Exhaustive Local Search (ELS) [2]. OMPR and ELS are initialized with the solution of OMP. Three versions of SEA are studied: the cold-start version $SEA_0$, where SEA is initialized with the null vector, and the warm-start versions $SEA_{ELS}$ and $SEA_{OMP}$, where SEA is initialized with the solutions of ELS and OMP, respectively. We have also studied HTP and IHT initialized with OMP and ELS. They are called $HTP_{OMP}$, $HTP_{ELS}$, $IHT_{OMP}$ and $IHT_{ELS}$.

For all algorithms, each least-square projection for a fixed support, as in Line 7 of Algorithm 1, is solved using the conjugate gradient descent of SciPy [44]. For all algorithms, 256k iterations are performed. The results of HTP and to a lesser extent IHT and SEA depend on the choice of the step size. For the sake of fairness of the comparison with OMP, OMPR and ELS, we did not optimize the choice of the step size. The step size of SEA, HTP, and IHT is arbitrarily<fnref n="6" /> fixed to $\eta = \frac{1.8}{L}$, where $L$ is the spectral radius of $A$. The columns of $A$ are normalized before solving the problem. The sparse vector $x^* \in \mathbb{R}^n$ is random. Indexes of the support are randomly picked, uniformly without replacement. The non-zero entries of $x^*$ are drawn uniformly in $[-2, -1] \cup [1, 2]$ as in [24]. The noise $e$ is drawn uniformly using the same method as described in [7]. Their detailed descriptions are in the next two sections. For each experiment, the metrics used for performance evaluation are defined in the corresponding subsection. The code is implemented in Python 3 and is available in the git repository of the project <fnref n="7" />. As explained

<footnote n="6">⁶We do not report further experiments for $\eta \in \{\frac{2l}{L} | \text{ for } l = [-3, +3]\}$ that do not significantly alter the results in terms of running time, stability, performance, and do not impact our conclusions on phase transition diagrams. Variation of the step size for IHT and HTP in the deconvolution experiment is reported in Appendix E.10 and does not impact our conclusions.</footnote>

<footnote n="7">⁷https://gitlab.lis-lab.fr/valentin.emiya/sea-icml-2024</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2: Phase transition diagram: each curve is the threshold below which the related algorithm recovers at least 95% of the supports. $\zeta$ denotes the ratio between the number of rows and the number of columns in A while $\rho$ denotes the ratio between the sparsity and the number of rows in A. Matrix A have i.i.d. standard Gaussian entries and non-zero entries in $x^*$ are drawn uniformly in $[-2, -1] \cup [1, 2]$. $n = 500$ is fixed and results are obtained from 1000 runs.</figcaption></figure>

in Section 3.3 and in Appendix B.2, the computational cost of SEA mainly depends on the number of explored supports. The illustration related to the number of explored supports for a fixed number of iterations and the efficiency of the exploration can be found in Appendix E.3.

## 5.1 Phase Transition Diagram Experiment

Introduced by Donoho and Tanner [19] and used in compressed sensing [21, 26], phase transition diagrams show the recovery limits of an algorithm depending on the undersampling/indeterminacy $\zeta = \frac{m}{n}$ of A, and the sparsity/density $\rho = \frac{k}{m}$ of $x^*$.

We fix $n = 500$, $m$ takes 18 values in $\llbracket 1, n \rrbracket$ and $k$ all values in $\llbracket 1, 0.5m \rrbracket$. For each triplet $(m, n, k)$ and each algorithm, we run $r = 1000$ experiments (described below) to assess the success rate $\frac{s_{\zeta, \rho}}{r}$ of the algorithm, where $s_{\zeta, \rho}$ is the number of problems successfully solved. A problem is considered successfully solved if the support of the output of the algorithm is equal to $S^*$. For each run, the entries of $A \in \mathbb{R}^{m \times n}$ are drawn independently from the standard normal distribution. The restricted isometry constants are poor when $\zeta = \frac{m}{n}$ is small and improve when $m$ grows [3]. The noise $e$ is drawn uniformly from the sphere of radius $0.01 \|Ax^*\|_2$ in $\mathbb{R}^m$.

Figure 2 shows results from this experiment. Each curve indicates the threshold below which the algorithm has a success rate larger than 95%. The higher the curve, the better. We see that OMP, HTP and IHT achieve poor recovery successes. The smooth, decreasing part of the HTP and IHT curves on the left is an artifact due to the discrete values of $(m, n, k)$ and actually corresponds to a phase transition located at $k = 1$. $SEA_0$ outperforms OMP, HTP and IHT when $\frac{m}{n} < 0.6$. All the OMP-initialized algorithms (in blue) improve OMP performance except in the most coherent cases ($\frac{m}{n} < 0.2$) where $HTP_{OMP}$ and $IHT_{OMP}$ fail while $SEA_{OMP}$ exhibits the best improvement. Contrary to $HTP_{ELS}$ and $IHT_{ELS}$, $SEA_{ELS}$ (in red) improves further ELS performances and outperforms the other algorithms for all $\frac{m}{n}$. The main improvements are when $\frac{m}{n}$ is small ($\frac{m}{n} < 0.4$), i.e., for the most coherent matrices A. Thus, SEA refines a good support candidate into a better one by exploring new supports and achieves recovery for higher values of sparsity $k$ than competitors. The actual superiority of $SEA_{ELS}$ and $SEA_{OMP}$ for coherent matrices ($\frac{m}{n} < 0.3$) is a major conclusion from this experiment and illustrates its ability to successfully explore supports in difficult problems where competitors fail. We study the noiseless setup (i.e., $e = 0$) in Appendix D.

## 5.2 Deconvolution Experiment

Deconvolution purposes arise in many signal processing areas such as microscopy or remote sensing. Of particular interest is the deconvolution of sparse signals, also known as point source deconvolution [6] or
<figure><img src="image_5.png" /><figcaption>Figure 2: Phase transition diagram: each curve is the threshold below which the related algorithm recovers at least 95% of the supports. $\zeta$ denotes the ratio between the number of rows and the number of columns in A while $\rho$ denotes the ratio between the sparsity and the number of rows in A. Matrix A have i.i.d. standard Gaussian entries and non-zero entries in $x^*$ are drawn uniformly in $[-2, -1] \cup [1, 2]$. $n = 500$ is fixed and results are obtained from 1000 runs.</figcaption></figure>

in Section 3.3 and in Appendix B.2, the computational cost of SEA mainly depends on the number of explored supports. The illustration related to the number of explored supports for a fixed number of iterations and the efficiency of the exploration can be found in Appendix E.3.

## 5.1 Phase Transition Diagram Experiment

Introduced by Donoho and Tanner [19] and used in compressed sensing [21, 26], phase transition diagrams show the recovery limits of an algorithm depending on the undersampling/indeterminacy $\zeta = \frac{m}{n}$ of A, and the sparsity/density $\rho = \frac{k}{m}$ of $x^*$.

We fix $n = 500$, $m$ takes 18 values in $\llbracket 1, n \rrbracket$ and $k$ all values in $\llbracket 1, 0.5m \rrbracket$. For each triplet $(m, n, k)$ and each algorithm, we run $r = 1000$ experiments (described below) to assess the success rate $\frac{s_{\zeta, \rho}}{r}$ of the algorithm, where $s_{\zeta, \rho}$ is the number of problems successfully solved. A problem is considered successfully solved if the support of the output of the algorithm is equal to $S^*$. For each run, the entries of $A \in \mathbb{R}^{m \times n}$ are drawn independently from the standard normal distribution. The restricted isometry constants are poor when $\zeta = \frac{m}{n}$ is small and improve when $m$ grows [3]. The noise $e$ is drawn uniformly from the sphere of radius $0.01 \|Ax^*\|_2$ in $\mathbb{R}^m$.

Figure 2 shows results from this experiment. Each curve indicates the threshold below which the algorithm has a success rate larger than 95%. The higher the curve, the better. We see that OMP, HTP and IHT achieve poor recovery successes. The smooth, decreasing part of the HTP and IHT curves on the left is an artifact due to the discrete values of $(m, n, k)$ and actually corresponds to a phase transition located at $k = 1$. $SEA_0$ outperforms OMP, HTP and IHT when $\frac{m}{n} < 0.6$. All the OMP-initialized algorithms (in blue) improve OMP performance except in the most coherent cases ($\frac{m}{n} < 0.2$) where $HTP_{OMP}$ and $IHT_{OMP}$ fail while $SEA_{OMP}$ exhibits the best improvement. Contrary to $HTP_{ELS}$ and $IHT_{ELS}$, $SEA_{ELS}$ (in red) improves further ELS performances and outperforms the other algorithms for all $\frac{m}{n}$. The main improvements are when $\frac{m}{n}$ is small ($\frac{m}{n} < 0.4$), i.e., for the most coherent matrices A. Thus, SEA refines a good support candidate into a better one by exploring new supports and achieves recovery for higher values of sparsity $k$ than competitors. The actual superiority of $SEA_{ELS}$ and $SEA_{OMP}$ for coherent matrices ($\frac{m}{n} < 0.3$) is a major conclusion from this experiment and illustrates its ability to successfully explore supports in difficult problems where competitors fail. We study the noiseless setup (i.e., $e = 0$) in Appendix D.

## 5.2 Deconvolution Experiment

Deconvolution purposes arise in many signal processing areas such as microscopy or remote sensing. Of particular interest is the deconvolution of sparse signals, also known as point source deconvolution [6] or
<figure><img src="image_5.png" /><figcaption>Figure 3: Spike deconvolution: representation of an instance of $x^*$ and $y$ with the solutions provided by the algorithms when $k = 20$. This is a cropped version of a crowded area (spikes are close).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Spike deconvolution: average support distance between $S^*$ and the support of the solutions provided by several algorithms as a function of the sparsity level $k$.</figcaption></figure>

spike deconvolution [23, 22], assuming the linear operator is known (contrary to blind approaches [33]). The objective is to recover spike positions and amplitudes.

We set $n = 500$, a convolution matrix $A$ corresponding to a Gaussian filter with a standard deviation equal to 3. The coherence of matrix $A$ is $\max_{i \neq j} |A_i^T A_j| = 0.97$, resulting in very difficult problems for which the support recovery theorems do not apply. For each sparsity level $k \in [1, 50]$, every algorithm is tested on $r = 200$ distinct problems corresponding to different $k$-sparse vectors $x^*$. The maximal number of iterations is 1000, for all algorithms. The noise $e$ is drawn uniformly from the sphere of radius $0.1 \|Ax^*\|_2$ of $\mathbb{R}^m$, aiming for a signal-to-noise ratio of 20dB. Cases where sparsity is wrongly estimated and noise is stronger or applied differently are studied in Appendices E.6 to E.8.

Figure 3 illustrates the results for a 20-sparse vector $x^*$ restricted to a crowded area of the full signal (the later being depicted in Appendix E.1). Generally speaking, isolated spikes are recovered by almost all algorithms. However, algorithms often fail to accurately identify spikes when they are close to each other. For instance, ELS, OMP and OMPR falsely detect entries in the highest energetic part of the signal (around position 400) and are trapped in a local minimum. $\text{SEA}_0$, $\text{SEA}_{\text{OMP}}$, and $\text{SEA}_{\text{ELS}}$ recover the original signal with a better precision than its competitors. It is worth mentioning that only $\text{SEA}$ recovers perfectly this signal in the noiseless settings (see Appendix E.5.2). To illustrate the exploratory behavior of $\text{SEA}$, we show in Appendix E.2, the evolution of $\|Ax^t - y\|_2$ when $t$ and the number of explored supports varies, for the experiment of Figure 3.

On Figure 4, for each algorithm and for all $k \in [1, 50]$, we display the support distance metric [24] averaged over $r = 200$ runs and defined by

$$ \text{dist}_{\text{supp}}(x) = \frac{k - |S^* \cap \text{supp}(x)|}{k} \qquad (8) $$
<figure><img src="image_6.png" /><figcaption>Figure 3: Spike deconvolution: representation of an instance of $x^*$ and $y$ with the solutions provided by the algorithms when $k = 20$. This is a cropped version of a crowded area (spikes are close).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: Spike deconvolution: average support distance between $S^*$ and the support of the solutions provided by several algorithms as a function of the sparsity level $k$.</figcaption></figure>

spike deconvolution [23, 22], assuming the linear operator is known (contrary to blind approaches [33]). The objective is to recover spike positions and amplitudes.

We set $n = 500$, a convolution matrix $A$ corresponding to a Gaussian filter with a standard deviation equal to 3. The coherence of matrix $A$ is $\max_{i \neq j} |A_i^T A_j| = 0.97$, resulting in very difficult problems for which the support recovery theorems do not apply. For each sparsity level $k \in [1, 50]$, every algorithm is tested on $r = 200$ distinct problems corresponding to different $k$-sparse vectors $x^*$. The maximal number of iterations is 1000, for all algorithms. The noise $e$ is drawn uniformly from the sphere of radius $0.1 \|Ax^*\|_2$ of $\mathbb{R}^m$, aiming for a signal-to-noise ratio of 20dB. Cases where sparsity is wrongly estimated and noise is stronger or applied differently are studied in Appendices E.6 to E.8.

Figure 3 illustrates the results for a 20-sparse vector $x^*$ restricted to a crowded area of the full signal (the later being depicted in Appendix E.1). Generally speaking, isolated spikes are recovered by almost all algorithms. However, algorithms often fail to accurately identify spikes when they are close to each other. For instance, ELS, OMP and OMPR falsely detect entries in the highest energetic part of the signal (around position 400) and are trapped in a local minimum. $\text{SEA}_0$, $\text{SEA}_{\text{OMP}}$, and $\text{SEA}_{\text{ELS}}$ recover the original signal with a better precision than its competitors. It is worth mentioning that only $\text{SEA}$ recovers perfectly this signal in the noiseless settings (see Appendix E.5.2). To illustrate the exploratory behavior of $\text{SEA}$, we show in Appendix E.2, the evolution of $\|Ax^t - y\|_2$ when $t$ and the number of explored supports varies, for the experiment of Figure 3.

On Figure 4, for each algorithm and for all $k \in [1, 50]$, we display the support distance metric [24] averaged over $r = 200$ runs and defined by

$$ \text{dist}_{\text{supp}}(x) = \frac{k - |S^* \cap \text{supp}(x)|}{k} \qquad (8) $$
(the lower the distance, the better). For all considered sparsity values, $SEA_0$, $SEA_{OMP}$, and $SEA_{ELS}$ outperform the other algorithms. SEA improves OMP and ELS results while they are never enhanced by HTP nor IHT (curves are superimposed). Note that for small $k$, IHT shows poor performance because it assigns several neighboring elements of the support to the largest peak of $y$ and fails to correct this error afterward. As $k$ increases, due to the increasing difficulty of the problem, the algorithms are gradually becoming unable to recover $S^*$. Using a cold-start strategy, $SEA_0$ is here the best performing algorithm. The analyses conducted in Appendix E.3 indicate that the exploration carried out by SEA can be more efficient than the support element swaps performed by ELS. These experiments also suggest that a warm-start strategy, such as $SEA_{ELS}$ or $SEA_{OMP}$, may lead the algorithm to get trapped in a local minimum. The choice of the best strategy appears to depend on the quality of the initialization. We recommend selecting it based on empirical performance. The same conclusions are drawn when using additional metrics (Appendix E.4) and in the noiseless case (Appendix E.5).

# 6 Conclusions and Perspectives

In this article, we proposed SEA: a new principled algorithm for sparse support recovery, based on STE. Experiments show that SEA supplements state-of-the-art algorithms and outperforms them in particular when $A$ is coherent, thanks to its better exploration ability. Indeed, SEA initialized with the output of ELS is generally a good strategy to try to improve recovery results. Nonetheless, the cold-start strategy where SEA is initialized at 0 may also be profitable: it is the best setting in problems with very coherent matrices like in the deconvolution experiment. Understanding which strategy should be preferred remains an open question.

We established guarantees when the matrix $A$ satisfies the RIP, which we hope gives new insight on the STE. The theoretical guarantees involve conditions on $x^*$ that are not present for similar statements for other algorithms and that might restrict their applicability. Improving the theoretical analysis in the following directions are promising perspective. The algorithm perform well when $A$ is coherent: this is not explained by the current theoretical analysis which only applies to matrices satisfying the RIP. Checking explicitly RIP conditions being NP-hard [43], we will investigate theoretical guarantees based on mutual incoherence [18]. Also it would be interesting to adapt the strategy developed for obtaining the theoretical guarantees to other contexts, such as the optimization of quantized neural networks.

Finally, this paper opens up broader perspectives. The proposed STE is a deterministic approach for support exploration and may also be compared to or extended by the use of stochastic heuristics. Also, it would be of interest to study, either theoretically or numerically, the behavior of SEA in the compressed sensing setting. There are many perspectives of SEA and STE applications to sparse inverse problems such as sparse matrix factorization, tensor problems, as well as real-world applications such as in biology and astronomy.

# Acknowledgement

F. Malgouyres gratefully acknowledges the support of IRT Saint Exupéry and the DEEL project<fnref n="8" /> and thanks Franck Mamalet for all the discussions on the STE.

M. Mohamed was supported by a PhD grant from “Emploi Jeunes Doctorants (EJD)” plan which is funded by the French institution “Conseil régional Provence-Alpes-Côte d’Azur” and Euranova France. M. Mohamed gratefully acknowledges their financial support.

# References

[1] Alexander G. Anderson and Cory P. Berg. The high-dimensional geometry of binary neural networks. In *International Conference on Learning Representations*, 2018.

[2] Kyriakos Axiotis and Maxim Sviridenko. Sparse convex optimization via adaptively regularized hard thresholding. In *Proc. Int. Conf. Mach. Learn.*, volume 119 of *Proceedings of Machine Learning Research*, pages 452–462. PMLR, 13–18 Jul 2020.

[3] Bubacarr Bah and Jared Tanner. Bounds of restricted isometry constants in extreme asymptotics: formulae for Gaussian matrices. *Lin. Algebra Appl.*, 441:88–109, 2014.

<footnote n="8">⁸https://www.deel.ai/</footnote>
(the lower the distance, the better). For all considered sparsity values, $SEA_0$, $SEA_{OMP}$, and $SEA_{ELS}$ outperform the other algorithms. SEA improves OMP and ELS results while they are never enhanced by HTP nor IHT (curves are superimposed). Note that for small $k$, IHT shows poor performance because it assigns several neighboring elements of the support to the largest peak of $y$ and fails to correct this error afterward. As $k$ increases, due to the increasing difficulty of the problem, the algorithms are gradually becoming unable to recover $S^*$. Using a cold-start strategy, $SEA_0$ is here the best performing algorithm. The analyses conducted in Appendix E.3 indicate that the exploration carried out by SEA can be more efficient than the support element swaps performed by ELS. These experiments also suggest that a warm-start strategy, such as $SEA_{ELS}$ or $SEA_{OMP}$, may lead the algorithm to get trapped in a local minimum. The choice of the best strategy appears to depend on the quality of the initialization. We recommend selecting it based on empirical performance. The same conclusions are drawn when using additional metrics (Appendix E.4) and in the noiseless case (Appendix E.5).

# 6 Conclusions and Perspectives

In this article, we proposed SEA: a new principled algorithm for sparse support recovery, based on STE. Experiments show that SEA supplements state-of-the-art algorithms and outperforms them in particular when $A$ is coherent, thanks to its better exploration ability. Indeed, SEA initialized with the output of ELS is generally a good strategy to try to improve recovery results. Nonetheless, the cold-start strategy where SEA is initialized at 0 may also be profitable: it is the best setting in problems with very coherent matrices like in the deconvolution experiment. Understanding which strategy should be preferred remains an open question.

We established guarantees when the matrix $A$ satisfies the RIP, which we hope gives new insight on the STE. The theoretical guarantees involve conditions on $x^*$ that are not present for similar statements for other algorithms and that might restrict their applicability. Improving the theoretical analysis in the following directions are promising perspective. The algorithm perform well when $A$ is coherent: this is not explained by the current theoretical analysis which only applies to matrices satisfying the RIP. Checking explicitly RIP conditions being NP-hard [43], we will investigate theoretical guarantees based on mutual incoherence [18]. Also it would be interesting to adapt the strategy developed for obtaining the theoretical guarantees to other contexts, such as the optimization of quantized neural networks.

Finally, this paper opens up broader perspectives. The proposed STE is a deterministic approach for support exploration and may also be compared to or extended by the use of stochastic heuristics. Also, it would be of interest to study, either theoretically or numerically, the behavior of SEA in the compressed sensing setting. There are many perspectives of SEA and STE applications to sparse inverse problems such as sparse matrix factorization, tensor problems, as well as real-world applications such as in biology and astronomy.

# Acknowledgement

F. Malgouyres gratefully acknowledges the support of IRT Saint Exupéry and the DEEL project<sup>8</sup> and thanks Franck Mamalet for all the discussions on the STE.

M. Mohamed was supported by a PhD grant from “Emploi Jeunes Doctorants (EJD)” plan which is funded by the French institution “Conseil régional Provence-Alpes-Côte d’Azur” and Euranova France. M. Mohamed gratefully acknowledges their financial support.

# References

[1] Alexander G. Anderson and Cory P. Berg. The high-dimensional geometry of binary neural networks. In *International Conference on Learning Representations*, 2018.

[2] Kyriakos Axiotis and Maxim Sviridenko. Sparse convex optimization via adaptively regularized hard thresholding. In *Proc. Int. Conf. Mach. Learn.*, volume 119 of *Proceedings of Machine Learning Research*, pages 452–462. PMLR, 13–18 Jul 2020.

[3] Bubacarr Bah and Jared Tanner. Bounds of restricted isometry constants in extreme asymptotics: formulae for Gaussian matrices. *Lin. Algebra Appl.*, 441:88–109, 2014.

<footnote>8https://www.deel.ai/</footnote>
[4] Ramzi Ben Mhenni, Sébastien Bourguignon, and Jordan Ninin. Global optimization for sparse solution of least squares problems. *Optim. Methods Softw.*, 37(5):1740–1769, 2022.

[5] Yoshua Bengio, Nicholas Léonard, and Aaron Courville. Estimating or propagating gradients through stochastic neurons for conditional computation. *CoRR*, arXiv:1308.3432, 2013.

[6] Brett Bernstein and Carlos Fernandez-Granda. Deconvolution of point sources: A sampling theorem and robustness guarantees. *Comm. Pure Appl. Math.*, 72(6):1152–1230, 2019.

[7] Jeffrey D Blanchard and Jared Tanner. Performance comparisons of greedy algorithms in compressed sensing. *Numerical Linear Algebra with Applications*, 22(2):254–282, 2015.

[8] Thomas Blumensath and Mike E. Davies. Iterative hard thresholding for compressed sensing. *Appl. Comput. Harmon. Analysis*, 27(3):265–274, 2009.

[9] Thomas Blumensath and Mike E Davies. Normalized iterative hard thresholding: Guaranteed stability and performance. *IEEE J. Sel. Topics Signal Process.*, 4(2):298–309, 2010.

[10] T Tony Cai and Lie Wang. Orthogonal matching pursuit for sparse signal recovery with noise. *IEEE Trans. Inform. Theory*, 57(7):4680–4688, 2011.

[11] Emmanuel J Candès, Justin Romberg, and Terence Tao. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. *IEEE Trans. Inform. Theory*, 52(2):489–509, 2006.

[12] Emmanuel J Candès and Terence Tao. Decoding by linear programming. *IEEE Trans. Inform. Theory*, 51(12):4203–4215, 2005.

[13] Scott Shaobing Chen, David L Donoho, and Michael A Saunders. Atomic decomposition by basis pursuit. *SIAM Rev.*, 43(1):129–159, 2001.

[14] Pengyu Cheng, Chang Liu, Chunyuan Li, Dinghan Shen, Ricardo Henao, and Lawrence Carin. Straight-through estimator as projected wasserstein gradient flow. In *Third workshop on Bayesian Deep Learning (NeurIPS)*, 2018.

[15] Matthieu Courbariaux, Yoshua Bengio, and Jean-Pierre David. Binaryconnect: Training deep neural networks with binary weights during propagations. In *Adv. Neural Inf. Process. Syst.*, volume 28, pages 3123–3131, Montreal, Quebec, Canada, Dec. 7–12 2015.

[16] Wei Dai and Olgica Milenkovic. Subspace pursuit for compressive sensing signal reconstruction. *IEEE Trans. Inform. Theory*, 55(5):2230–2249, 2009.

[17] Geoff Davis, Stephane Mallat, and Marco Avellaneda. Adaptive greedy approximations. *Constr. Approx.*, 13(1):57–98, 1997.

[18] David Donoho and X. Huo. Uncertainty principles and ideal atomic decomposition. *IEEE Trans. Inform. Theory*, 47(7):2845–2862, 2001.

[19] David Donoho and Jared Tanner. Observed universality of phase transitions in high-dimensional geometry, with implications for modern data analysis and signal processing. *Phil. Trans. R. Soc. A*, 367(1906):4273–4293, 2009.

[20] David L Donoho. Compressed sensing. *IEEE Trans. Inform. Theory*, 52(4):1289–1306, 2006.

[21] David L Donoho, Arian Maleki, and Andrea Montanari. Message-passing algorithms for compressed sensing. *Proceedings of the National Academy of Sciences*, 106(45):18914–18919, 2009.

[22] Vincent Duval and Gabriel Peyré. Sparse spikes super-resolution on thin grids I: the LASSO. *Inverse Problems*, 33(5):055008, 2017.

[23] Vincent Duval and Gabriel Peyré. Exact support recovery for sparse spikes deconvolution. *Found. Comput. Math.*, 15(5):1315–1355, 2015.

[24] Michael Elad. *Sparse and Redundant Representations*. Springer New York, NY, 1 edition, 2010.
[4] Ramzi Ben Mhenni, Sébastien Bourguignon, and Jordan Ninin. Global optimization for sparse solution of least squares problems. *Optim. Methods Softw.*, 37(5):1740–1769, 2022.

[5] Yoshua Bengio, Nicholas Léonard, and Aaron Courville. Estimating or propagating gradients through stochastic neurons for conditional computation. *CoRR*, arXiv:1308.3432, 2013.

[6] Brett Bernstein and Carlos Fernandez-Granda. Deconvolution of point sources: A sampling theorem and robustness guarantees. *Comm. Pure Appl. Math.*, 72(6):1152–1230, 2019.

[7] Jeffrey D Blanchard and Jared Tanner. Performance comparisons of greedy algorithms in compressed sensing. *Numerical Linear Algebra with Applications*, 22(2):254–282, 2015.

[8] Thomas Blumensath and Mike E. Davies. Iterative hard thresholding for compressed sensing. *Appl. Comput. Harmon. Analysis*, 27(3):265–274, 2009.

[9] Thomas Blumensath and Mike E Davies. Normalized iterative hard thresholding: Guaranteed stability and performance. *IEEE J. Sel. Topics Signal Process.*, 4(2):298–309, 2010.

[10] T Tony Cai and Lie Wang. Orthogonal matching pursuit for sparse signal recovery with noise. *IEEE Trans. Inform. Theory*, 57(7):4680–4688, 2011.

[11] Emmanuel J Candès, Justin Romberg, and Terence Tao. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. *IEEE Trans. Inform. Theory*, 52(2):489–509, 2006.

[12] Emmanuel J Candès and Terence Tao. Decoding by linear programming. *IEEE Trans. Inform. Theory*, 51(12):4203–4215, 2005.

[13] Scott Shaobing Chen, David L Donoho, and Michael A Saunders. Atomic decomposition by basis pursuit. *SIAM Rev.*, 43(1):129–159, 2001.

[14] Pengyu Cheng, Chang Liu, Chunyuan Li, Dinghan Shen, Ricardo Henao, and Lawrence Carin. Straight-through estimator as projected wasserstein gradient flow. In *Third workshop on Bayesian Deep Learning (NeurIPS)*, 2018.

[15] Matthieu Courbariaux, Yoshua Bengio, and Jean-Pierre David. Binaryconnect: Training deep neural networks with binary weights during propagations. In *Adv. Neural Inf. Process. Syst.*, volume 28, pages 3123–3131, Montreal, Quebec, Canada, Dec. 7–12 2015.

[16] Wei Dai and Olgica Milenkovic. Subspace pursuit for compressive sensing signal reconstruction. *IEEE Trans. Inform. Theory*, 55(5):2230–2249, 2009.

[17] Geoff Davis, Stephane Mallat, and Marco Avellaneda. Adaptive greedy approximations. *Constr. Approx.*, 13(1):57–98, 1997.

[18] David Donoho and X. Huo. Uncertainty principles and ideal atomic decomposition. *IEEE Trans. Inform. Theory*, 47(7):2845–2862, 2001.

[19] David Donoho and Jared Tanner. Observed universality of phase transitions in high-dimensional geometry, with implications for modern data analysis and signal processing. *Phil. Trans. R. Soc. A*, 367(1906):4273–4293, 2009.

[20] David L Donoho. Compressed sensing. *IEEE Trans. Inform. Theory*, 52(4):1289–1306, 2006.

[21] David L Donoho, Arian Maleki, and Andrea Montanari. Message-passing algorithms for compressed sensing. *Proceedings of the National Academy of Sciences*, 106(45):18914–18919, 2009.

[22] Vincent Duval and Gabriel Peyré. Sparse spikes super-resolution on thin grids I: the LASSO. *Inverse Problems*, 33(5):055008, 2017.

[23] Vincent Duval and Gabriel Peyré. Exact support recovery for sparse spikes deconvolution. *Found. Comput. Math.*, 15(5):1315–1355, 2015.

[24] Michael Elad. *Sparse and Redundant Representations*. Springer New York, NY, 1 edition, 2010.
[25] Simon Foucart. Hard thresholding pursuit: an algorithm for compressive sensing. *SIAM J. Numer. Anal.*, 49(6):2543–2563, 2011.

[26] Simon Foucart, Holger Rauhut, Simon Foucart, and Holger Rauhut. *An invitation to compressive sensing*. Springer, 2013.

[27] Trevor Hastie, Robert Tibshirani, and Martin Wainwright. *Statistical Learning with Sparsity*. Chapman and Hall/CRC, May 2015.

[28] Koen Helwegen, James Widdicombe, Lukas Geiger, Zechun Liu, Kwang-Ting Cheng, and Roeland Nusselder. Latent weights do not exist: Rethinking binarized neural network optimization. *Advances in neural information processing systems*, 32, 2019.

[29] Geoffrey E. Hinton. Neural networks for machine learning. Coursera, video lectures, 2012. Lecture 15b.

[30] Itay Hubara, Matthieu Courbariaux, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio. Binarized neural networks. *Advances in neural information processing systems*, 29, 2016.

[31] Itay Hubara, Matthieu Courbariaux, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio. Binarized neural networks. *Adv. Neural Inf. Process. Syst.*, 29, Dec. 5–10 2016.

[32] Prateek Jain, Ambuj Tewari, and Inderjit Dhillon. Orthogonal matching pursuit with replacement. *Adv. Neural Inf. Process. Syst.*, 24, Dec. 12–17 2011.

[33] Han-Wen Kuo, Yenson Lau, Yuqian Zhang, and John Wright. Geometry and symmetry in short-and-sparse deconvolution. In *Proc. Int. Conf. Mach. Learn.*, volume 97 of *Proceedings of Machine Learning Research*, pages 3570–3580. PMLR, 09–15 Jun 2019.

[34] Hao Li, Soham De, Zheng Xu, Christoph Studer, Hanan Samet, and Tom Goldstein. Training quantized nets: A deeper understanding. In *Advances in Neural Information Processing Systems*, pages 5811–5821, 2017.

[35] Ziang Long, Penghang Yin, and Jack Xin. Learning quantized neural nets by coarse gradient method for nonlinear classification. *Research in the Mathematical Sciences*, 8:1–19, 2021.

[36] Stéphane G Mallat and Zhifeng Zhang. Matching pursuits with time-frequency dictionaries. *IEEE Trans. Signal Process.*, 41(12):3397–3415, 1993.

[37] Nicolai Meinshausen and Peter Bühlmann. High-dimensional graphs and variable selection with the Lasso. *Ann. Statist.*, 34(3):1436–1462, 2006.

[38] Deanna Needell and Joel A Tropp. CoSaMP: Iterative signal recovery from incomplete and inaccurate samples. *Appl. Comput. Harmon. Analysis*, 26(3):301–321, 2009.

[39] Yagyensh Chandra Pati, Ramin Rezaiifar, and Perinkulam Sambamurthy Krishnaprasad. Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition. In *Proc. Asilomar Conf. Signal Syst. Comput.*, volume 1, pages 40–44, Pacific Grove, CA, USA, 1993. IEEE.

[40] Ratshih Sayed, Haytham Azmi, Heba A. Shawkey, Alaa Hussein Khalil, and Mohamed Refky. A systematic literature review on binary neural networks. *IEEE Access*, 11:27546–27578, 2023.

[41] Shai Shalev-Shwartz, Nathan Srebro, and Tong Zhang. Trading accuracy for sparsity in optimization problems with sparsity constraints. *SIAM J. Optim.*, 20(6):2807–2832, 2010.

[42] Robert Tibshirani. Regression shrinkage and selection via the lasso. *J. R. Stat. Soc. Ser. B Stat. Methodol.*, 58(1):267–288, 1996.

[43] Andreas M. Tillmann and Marc E. Pfetsch. The computational complexity of the restricted isometry property, the nullspace property, and related concepts in compressed sensing. *IEEE Trans. Inform. Theory*, 60(2):1248–1259, 2014.
[25] Simon Foucart. Hard thresholding pursuit: an algorithm for compressive sensing. *SIAM J. Numer. Anal.*, 49(6):2543–2563, 2011.

[26] Simon Foucart, Holger Rauhut, Simon Foucart, and Holger Rauhut. *An invitation to compressive sensing*. Springer, 2013.

[27] Trevor Hastie, Robert Tibshirani, and Martin Wainwright. *Statistical Learning with Sparsity*. Chapman and Hall/CRC, May 2015.

[28] Koen Helwegen, James Widdicombe, Lukas Geiger, Zechun Liu, Kwang-Ting Cheng, and Roeland Nusselder. Latent weights do not exist: Rethinking binarized neural network optimization. *Advances in neural information processing systems*, 32, 2019.

[29] Geoffrey E. Hinton. Neural networks for machine learning. Coursera, video lectures, 2012. Lecture 15b.

[30] Itay Hubara, Matthieu Courbariaux, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio. Binarized neural networks. *Advances in neural information processing systems*, 29, 2016.

[31] Itay Hubara, Matthieu Courbariaux, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio. Binarized neural networks. *Adv. Neural Inf. Process. Syst.*, 29, Dec. 5–10 2016.

[32] Prateek Jain, Ambuj Tewari, and Inderjit Dhillon. Orthogonal matching pursuit with replacement. *Adv. Neural Inf. Process. Syst.*, 24, Dec. 12–17 2011.

[33] Han-Wen Kuo, Yenson Lau, Yuqian Zhang, and John Wright. Geometry and symmetry in short-and-sparse deconvolution. In *Proc. Int. Conf. Mach. Learn.*, volume 97 of *Proceedings of Machine Learning Research*, pages 3570–3580. PMLR, 09–15 Jun 2019.

[34] Hao Li, Soham De, Zheng Xu, Christoph Studer, Hanan Samet, and Tom Goldstein. Training quantized nets: A deeper understanding. In *Advances in Neural Information Processing Systems*, pages 5811–5821, 2017.

[35] Ziang Long, Penghang Yin, and Jack Xin. Learning quantized neural nets by coarse gradient method for nonlinear classification. *Research in the Mathematical Sciences*, 8:1–19, 2021.

[36] Stéphane G Mallat and Zhifeng Zhang. Matching pursuits with time-frequency dictionaries. *IEEE Trans. Signal Process.*, 41(12):3397–3415, 1993.

[37] Nicolai Meinshausen and Peter Bühlmann. High-dimensional graphs and variable selection with the Lasso. *Ann. Statist.*, 34(3):1436–1462, 2006.

[38] Deanna Needell and Joel A Tropp. CoSaMP: Iterative signal recovery from incomplete and inaccurate samples. *Appl. Comput. Harmon. Analysis*, 26(3):301–321, 2009.

[39] Yagyensh Chandra Pati, Ramin Rezaiifar, and Perinkulam Sambamurthy Krishnaprasad. Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition. In *Proc. Asilomar Conf. Signal Syst. Comput.*, volume 1, pages 40–44, Pacific Grove, CA, USA, 1993. IEEE.

[40] Ratshih Sayed, Haytham Azmi, Heba A. Shawkey, Alaa Hussein Khalil, and Mohamed Refky. A systematic literature review on binary neural networks. *IEEE Access*, 11:27546–27578, 2023.

[41] Shai Shalev-Shwartz, Nathan Srebro, and Tong Zhang. Trading accuracy for sparsity in optimization problems with sparsity constraints. *SIAM J. Optim.*, 20(6):2807–2832, 2010.

[42] Robert Tibshirani. Regression shrinkage and selection via the lasso. *J. R. Stat. Soc. Ser. B Stat. Methodol.*, 58(1):267–288, 1996.

[43] Andreas M. Tillmann and Marc E. Pfetsch. The computational complexity of the restricted isometry property, the nullspace property, and related concepts in compressed sensing. *IEEE Trans. Inform. Theory*, 60(2):1248–1259, 2014.
[44] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nature Methods*, 17:261–272, 2020.

[45] Martin J Wainwright. Sharp thresholds for high-dimensional and noisy sparsity recovery using $ℓ_1$-constrained quadratic programming (Lasso). *IEEE Trans. Inform. Theory*, 55(5):2183–2202, 2009.

[46] Penghang Yin, Jiancheng Lyu, Shuai Zhang, Stanley J. Osher, Yingyong Qi, and Jack Xin. Understanding straight-through estimator in training activation quantized neural nets. In *International Conference on Learning Representations*, 2019.

[47] Chunyu Yuan and Sos S. Agaian. A comprehensive review of binary neural network. *Artificial Intelligence Review*, pages 1–65, 2023.

[48] Xiaotong Yuan, Ping Li, and Tong Zhang. Exact recovery of hard thresholding pursuit. In *Adv. Neural Inf. Process. Syst.*, volume 29, pages 3558–3566, Barcelona, Spain, Dec. 5–10 2016.

[49] Peng Zhao and Bin Yu. On model selection consistency of Lasso. *J. Mach. Learn. Res.*, 7:2541–2563, 2006.
[44] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nature Methods*, 17:261–272, 2020.

[45] Martin J Wainwright. Sharp thresholds for high-dimensional and noisy sparsity recovery using $ℓ_1$-constrained quadratic programming (Lasso). *IEEE Trans. Inform. Theory*, 55(5):2183–2202, 2009.

[46] Penghang Yin, Jiancheng Lyu, Shuai Zhang, Stanley J. Osher, Yingyong Qi, and Jack Xin. Understanding straight-through estimator in training activation quantized neural nets. In *International Conference on Learning Representations*, 2019.

[47] Chunyu Yuan and Sos S. Agaian. A comprehensive review of binary neural network. *Artificial Intelligence Review*, pages 1–65, 2023.

[48] Xiaotong Yuan, Ping Li, and Tong Zhang. Exact recovery of hard thresholding pursuit. In *Adv. Neural Inf. Process. Syst.*, volume 29, pages 3558–3566, Barcelona, Spain, Dec. 5–10 2016.

[49] Peng Zhao and Bin Yu. On model selection consistency of Lasso. *J. Mach. Learn. Res.*, 7:2541–2563, 2006.
# A Problem Statement

The equivalence between problem (1) and problem (2) is established by the following proposition. Before stating the proposition, let us remind

$$F(x) = \frac{1}{2} \|Ax - y\|_2^2, \quad \forall x \in \mathbb{R}^n \qquad (9)$$

and

$$H(\mathcal{X}) \in \underset{\substack{x \in \mathbb{R}^n \\ \text{supp}(x) \subseteq \text{largest}_k(\mathcal{X})}}{\operatorname{argmin}} \frac{1}{2} \|Ax - y\|_2^2, \quad \forall \mathcal{X} \in \mathbb{R}^n. \qquad (10)$$

Let us also recall the optimization problem (1)

$$\underset{x \in \mathbb{R}^n}{\text{Minimize}} F(x) \quad \text{s.t.} \quad \|x\|_0 \le k \qquad (11)$$

and the optimization problem (2)

$$\underset{\mathcal{X} \in \mathbb{R}^n}{\text{Minimize}} F(H(\mathcal{X})). \qquad (12)$$

**Proposition A.1 (Optimization problem equivalence).** For all $m, n, k \in \mathbb{N}$, $A \in \mathbb{R}^{m \times n}$ and $y \in \mathbb{R}^m$. Problem (11) is equivalent to problem (12), in the sense that

1. for any solution $\mathcal{X}^* \in \operatorname{argmin}_{\mathcal{X} \in \mathbb{R}^n} F(H(\mathcal{X}))$ of (12), $H(\mathcal{X}^*)$ is solution of (11).

2. for any minimizer $x'$ of (11), we have $x' \in \operatorname{argmin}_{\mathcal{X} \in \mathbb{R}^n} F(H(\mathcal{X}))$, i.e., $x'$ is solution of (12).

*Proof.* To establish the first item, we consider a solution $\mathcal{X}^* \in \operatorname{argmin}_{\mathcal{X} \in \mathbb{R}^n} F(H(\mathcal{X}))$ of (12). By definition of $H$, in (10), $H(\mathcal{X}^*)$ is $k$-sparse. To prove that it minimizes (11), consider $x \in \mathbb{R}^n$ such that $\|x\|_0 \le k$, we have

$$F(H(\mathcal{X}^*)) \le F(H(x)) \le F(x), \qquad (13)$$

where the first inequality is due to the hypothesis on $\mathcal{X}^*$, and the last inequality to the definition of $H$. Finally, since $H(\mathcal{X}^*)$ is $k$-sparse and (13) holds for all $k$-sparse vector $x$, we conclude that $H(\mathcal{X}^*)$ is solution of (11).

To prove the second item, consider a minimizer $x'$ of (11) and $\mathcal{X} \in \mathbb{R}^n$. By definition of $H$, $H(\mathcal{X})$ is $k$-sparse. Using that $x'$ is solution of (11), we therefore have

$$F(x') \le F(H(\mathcal{X})). \qquad (14)$$

Moreover, since $x'$ is $k$-sparse, we have $\text{supp}(x') \subseteq \text{largest}_k(x')$, and by the definition of $H$,

$$F(H(x')) \le F(x').$$

Combining with (14), we obtain $F(H(x')) \le F(H(\mathcal{X}))$, for all $\mathcal{X} \in \mathbb{R}^n$, and conclude that $x'$ is solution of (12). □

# B Additional Algorithms

In this appendix, more details are given about SEA pseudo-code: the main differences with state-of-the-art algorithms HTP and IHT are discussed in Section B.1 and tricks for an efficient implementation of SEA are given in Section B.2.

## B.1 State-of-the-Art Algorithms

In terms of pseudo-code, SEA looks similar to Hard Thresholding Pursuit (Algorithm 3, [25]) and to a less extent to Iterative Hard Thresholding (Algorithm 4, [8]). In this section, we highlight the differences between these algorithms. In particular, in Algorithm 2, Algorithm 3 and Algorithm 4 distinctions are pointed out in red.

Both HTP and IHT are projected descent algorithms that alternate a gradient step at a sparse estimate $x^t$ and a projection of the resulting variable $\mathcal{X}^t$ onto the set of sparse vectors. The whole difference with SEA lies in the introduction of the support exploration variable $\mathcal{X}^t$ and its interaction with the
# A Problem Statement

The equivalence between problem (1) and problem (2) is established by the following proposition. Before stating the proposition, let us remind

$$F(x) = \frac{1}{2} \|Ax - y\|_2^2, \quad \forall x \in \mathbb{R}^n \qquad (9)$$

and

$$H(\mathcal{X}) \in \underset{\substack{x \in \mathbb{R}^n \\ \text{supp}(x) \subseteq \text{largest}_k(\mathcal{X})}}{\operatorname{argmin}} \frac{1}{2} \|Ax - y\|_2^2, \quad \forall \mathcal{X} \in \mathbb{R}^n. \qquad (10)$$

Let us also recall the optimization problem (1)

$$\underset{x \in \mathbb{R}^n}{\text{Minimize}} F(x) \quad \text{s.t.} \quad \|x\|_0 \le k \qquad (11)$$

and the optimization problem (2)

$$\underset{\mathcal{X} \in \mathbb{R}^n}{\text{Minimize}} F(H(\mathcal{X})). \qquad (12)$$

**Proposition A.1 (Optimization problem equivalence).** For all $m, n, k \in \mathbb{N}$, $A \in \mathbb{R}^{m \times n}$ and $y \in \mathbb{R}^m$. Problem (11) is equivalent to problem (12), in the sense that

1. for any solution $\mathcal{X}^* \in \operatorname{argmin}_{\mathcal{X} \in \mathbb{R}^n} F(H(\mathcal{X}))$ of (12), $H(\mathcal{X}^*)$ is solution of (11).

2. for any minimizer $x'$ of (11), we have $x' \in \operatorname{argmin}_{\mathcal{X} \in \mathbb{R}^n} F(H(\mathcal{X}))$, i.e., $x'$ is solution of (12).

*Proof.* To establish the first item, we consider a solution $\mathcal{X}^* \in \operatorname{argmin}_{\mathcal{X} \in \mathbb{R}^n} F(H(\mathcal{X}))$ of (12). By definition of $H$, in (10), $H(\mathcal{X}^*)$ is $k$-sparse. To prove that it minimizes (11), consider $x \in \mathbb{R}^n$ such that $\|x\|_0 \le k$, we have

$$F(H(\mathcal{X}^*)) \le F(H(x)) \le F(x), \qquad (13)$$

where the first inequality is due to the hypothesis on $\mathcal{X}^*$, and the last inequality to the definition of $H$. Finally, since $H(\mathcal{X}^*)$ is $k$-sparse and (13) holds for all $k$-sparse vector $x$, we conclude that $H(\mathcal{X}^*)$ is solution of (11).

To prove the second item, consider a minimizer $x'$ of (11) and $\mathcal{X} \in \mathbb{R}^n$. By definition of $H$, $H(\mathcal{X})$ is $k$-sparse. Using that $x'$ is solution of (11), we therefore have

$$F(x') \le F(H(\mathcal{X})). \qquad (14)$$

Moreover, since $x'$ is $k$-sparse, we have $\text{supp}(x') \subseteq \text{largest}_k(x')$, and by the definition of $H$,

$$F(H(x')) \le F(x').$$

Combining with (14), we obtain $F(H(x')) \le F(H(\mathcal{X}))$, for all $\mathcal{X} \in \mathbb{R}^n$, and conclude that $x'$ is solution of (12).

□

# B Additional Algorithms

In this appendix, more details are given about SEA pseudo-code: the main differences with state-of-the-art algorithms HTP and IHT are discussed in Section B.1 and tricks for an efficient implementation of SEA are given in Section B.2.

## B.1 State-of-the-Art Algorithms

In terms of pseudo-code, SEA looks similar to Hard Thresholding Pursuit (Algorithm 3, [25]) and to a less extent to Iterative Hard Thresholding (Algorithm 4, [8]). In this section, we highlight the differences between these algorithms. In particular, in Algorithm 2, Algorithm 3 and Algorithm 4 distinctions are pointed out in red.

Both HTP and IHT are projected descent algorithms that alternate a gradient step at a sparse estimate $x^t$ and a projection of the resulting variable $\mathcal{X}^t$ onto the set of sparse vectors. The whole difference with SEA lies in the introduction of the support exploration variable $\mathcal{X}^t$ and its interaction with the
sparse vector $x^t$. HTP and IHT perform a regular gradient step $\mathcal{X}^{t+1} \leftarrow x^t - \eta A^T (Ax^t - y)$ (where $\mathcal{X}$ denotes an intermediate variable here, not a support exploration variable) while SEA uses an STE update $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - \eta A^T (Ax^t - y)$ of the support exploration variable itself ($\mathcal{X}^t$) with a gradient computed at $x^t$. As a consequence, the vector $\mathcal{X}^t$ in HTP or IHT is always one gradient step away from sparse vector $x^t$. They do not explore much. This is not the case with SEA. The support exploration variable $\mathcal{X}^t$ is not expected to minimize the objective: it rather accumulates all the gradient iterates, where the gradient is computed at $x^t$. This is the whole point of the STE. In particular, $(\mathcal{X}^t)_{t \in \mathbb{N}}$ is not restricted to a small portion of $\mathbb{R}^n$ in the vicinity of sparse vectors. It can explore much more than in HTP and IHT. This explains why SEA has a different exploration/exploitation trade-off. It explores more. As can also be seen from the experiments in Appendix E.5.2, the loss oscillates a lot during SEA's iterative process, but SEA retains the best solution $x^{t_{BEST}}$ encountered during the exploration. SEA is not based on a descent principle as IHT, HTP and such.

Finally, one may also notice that HTP stops as soon as the gradient is small enough such that the support does not change during two successive iterations. On the contrary, SEA keeps accumulating gradients so that the support may remain unchanged for many iterations before a new support is explored. This is clearly visible in the illustrations of Appendix E.5.2.

<caption>**Algorithm 2 SEA**
(copy of Algorithm 1)</caption>

1: **Inputs:**
   noisy observation $y$,
   sampling matrix $A$,
   sparsity $k$,
   step size $\eta$
2: **Output:** sparse vector $x$
3: Initialize $\mathcal{X}^0$
4: $t \leftarrow 0$
5: **repeat**
6:   $S^t \leftarrow \text{largest}_k(\mathcal{X}^t)$
7:   $\begin{cases} x_i^t & \leftarrow 0 \text{ for } i \in \bar{S}^t \\ x_{S^t}^t & \leftarrow A_{S^t}^\dagger y \end{cases}$
8:   $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - \eta A^T (Ax^t - y)$
9:   $t \leftarrow t + 1$
10: **until** halting criterion is true
11:   $t_{BEST} \leftarrow \underset{t \in [0,t]}{\arg\min} \|Ax^t - y\|_2^2$
12: **return** $x^{t_{BEST}}$

<caption>**Algorithm 3 HTP**
[25]</caption>

1: **Inputs:**
   noisy observation $y$,
   sampling matrix $A$,
   sparsity $k$,
   step size $\eta$
2: **Output:** sparse vector $x$
3: Initialize $\mathcal{X}^0$
4: $t \leftarrow 0$
5: **repeat**
6:   $S^t \leftarrow \text{largest}_k(\mathcal{X}^t)$
7:   $\begin{cases} x_i^t & \leftarrow 0 \text{ for } i \in \bar{S}^t \\ x_{S^t}^t & \leftarrow A_{S^t}^\dagger y \end{cases}$
8:   $\mathcal{X}^{t+1} \leftarrow x^t - \eta A^T (Ax^t - y)$
9:   $t \leftarrow t + 1$
10: **until** halting criterion is true
11: **return** $x^t$

<caption>**Algorithm 4 IHT**
[8]</caption>

1: **Inputs:**
   noisy observation $y$,
   sampling matrix $A$,
   sparsity $k$,
   step size $\eta$
2: **Output:** sparse vector $x$
3: Initialize $\mathcal{X}^0$
4: $t \leftarrow 0$
5: **repeat**
6:   $S^t \leftarrow \text{largest}_k(\mathcal{X}^t)$
7:   $\begin{cases} x_i^t & \leftarrow 0 \text{ for } i \in \bar{S}^t \\ x_{S^t}^t & \leftarrow \mathcal{X}_{S^t}^\dagger \end{cases}$
8:   $\mathcal{X}^{t+1} \leftarrow x^t - \eta A^T (Ax^t - y)$
9:   $t \leftarrow t + 1$
10: **until** halting criterion is true
11: **return** $x^t$

## B.2 Efficient Implementation of SEA

Algorithm 1 is presented in a way that favors clarity and simplifies the theoretical analysis. In practice, one can notice that if the support $S^t$ does not change (line 6), then the sparse vector $x^t$ and the gradient $\eta A^T (Ax^t - y)$ do not change either. Algorithm 5 is an equivalent pseudo-code for a computationally-efficient implementation. The most expensive computations —the sparse projection at line 10 and the gradient at line 12— are only required when the support has never been explored before. Also, the best sparse vector can be memorized on the fly (line 16). Hence, the remaining operations, that are performed at each iteration, have a low computational cost: support extraction (line 8), search for a previous, identical support (line 9) and STE update (line 21). This computationally-efficient version of SEA has a larger spatial complexity due to the memorization of all the supports and gradients seen along the iterations. However, this overhead is limited since 1/ for each explored support, only two vectors are memorized, one of them being sparse; and 2/ the number of explored supports is generally much lower than the number of iterations. For instance, in the deconvolution experiment, on average, less than 1000 vectors of size 500 (including 500 $k$-sparse vectors) are stored during the running time of SEA.

In addition, solving the (unconstrained) linear system $A_{S^t}^T A_{S^t} x_{S^t} = A_{S^t}^T y$ can also be performed efficiently. While the pseudo-inverse is a convenient notation at line 10, the solution may be obtained more efficiently, e.g. in $O(k^2n)$ to compute $A_{S^t}^T A_{S^t}$ and $A_{S^t}^T y$, and apply the conjugate gradient algorithm. This complexity is a worst-case scenario so in practice, the solution is generally obtained more quickly.
sparse vector $x^t$. HTP and IHT perform a regular gradient step $\mathcal{X}^{t+1} \leftarrow x^t - \eta A^T (Ax^t - y)$ (where $\mathcal{X}$ denotes an intermediate variable here, not a support exploration variable) while SEA uses an STE update $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - \eta A^T (Ax^t - y)$ of the support exploration variable itself ($\mathcal{X}^t$) with a gradient computed at $x^t$. As a consequence, the vector $\mathcal{X}^t$ in HTP or IHT is always one gradient step away from sparse vector $x^t$. They do not explore much. This is not the case with SEA. The support exploration variable $\mathcal{X}^t$ is not expected to minimize the objective: it rather accumulates all the gradient iterates, where the gradient is computed at $x^t$. This is the whole point of the STE. In particular, $(\mathcal{X}^t)_{t \in \mathbb{N}}$ is not restricted to a small portion of $\mathbb{R}^n$ in the vicinity of sparse vectors. It can explore much more than in HTP and IHT. This explains why SEA has a different exploration/exploitation trade-off. It explores more. As can also be seen from the experiments in Appendix E.5.2, the loss oscillates a lot during SEA's iterative process, but SEA retains the best solution $x^{t_{BEST}}$ encountered during the exploration. SEA is not based on a descent principle as IHT, HTP and such.

Finally, one may also notice that HTP stops as soon as the gradient is small enough such that the support does not change during two successive iterations. On the contrary, SEA keeps accumulating gradients so that the support may remain unchanged for many iterations before a new support is explored. This is clearly visible in the illustrations of Appendix E.5.2.

<caption>**Algorithm 2 SEA**
(copy of Algorithm 1)</caption>

1: **Inputs:**
   noisy observation $y$,
   sampling matrix $A$,
   sparsity $k$,
   step size $\eta$
2: **Output:** sparse vector $x$
3: Initialize $\mathcal{X}^0$
4: $t \leftarrow 0$
5: **repeat**
6:   $S^t \leftarrow \text{largest}_k(\mathcal{X}^t)$
7:   $\left\{ \begin{array}{ll} x_i^t & \leftarrow 0 \text{ for } i \in \bar{S}^t \\ x_{S^t}^t & \leftarrow A_{S^t}^\dagger y \end{array} \right.$
8:   $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - \eta A^T (Ax^t - y)$
9:   $t \leftarrow t + 1$
10: **until** halting criterion is true
11:   $t_{BEST} \leftarrow \underset{t \in [0,t]}{\arg\min} \|Ax^t - y\|_2^2$
12: **return** $x^{t_{BEST}}$

<caption>**Algorithm 3 HTP**
[25]</caption>

1: **Inputs:**
   noisy observation $y$,
   sampling matrix $A$,
   sparsity $k$,
   step size $\eta$
2: **Output:** sparse vector $x$
3: Initialize $\mathcal{X}^0$
4: $t \leftarrow 0$
5: **repeat**
6:   $S^t \leftarrow \text{largest}_k(\mathcal{X}^t)$
7:   $\left\{ \begin{array}{ll} x_i^t & \leftarrow 0 \text{ for } i \in \bar{S}^t \\ x_{S^t}^t & \leftarrow A_{S^t}^\dagger y \end{array} \right.$
8:   $\mathcal{X}^{t+1} \leftarrow x^t - \eta A^T (Ax^t - y)$
9:   $t \leftarrow t + 1$
10: **until** halting criterion is true
11: **return** $x^t$

<caption>**Algorithm 4 IHT**
[8]</caption>

1: **Inputs:**
   noisy observation $y$,
   sampling matrix $A$,
   sparsity $k$,
   step size $\eta$
2: **Output:** sparse vector $x$
3: Initialize $\mathcal{X}^0$
4: $t \leftarrow 0$
5: **repeat**
6:   $S^t \leftarrow \text{largest}_k(\mathcal{X}^t)$
7:   $\left\{ \begin{array}{ll} x_i^t & \leftarrow 0 \text{ for } i \in \bar{S}^t \\ x_{S^t}^t & \leftarrow \mathcal{X}_{S^t}^t \end{array} \right.$
8:   $\mathcal{X}^{t+1} \leftarrow x^t - \eta A^T (Ax^t - y)$
9:   $t \leftarrow t + 1$
10: **until** halting criterion is true
11: **return** $x^t$

## B.2 Efficient Implementation of SEA

Algorithm 1 is presented in a way that favors clarity and simplifies the theoretical analysis. In practice, one can notice that if the support $S^t$ does not change (line 6), then the sparse vector $x^t$ and the gradient $\eta A^T (Ax^t - y)$ do not change either. Algorithm 5 is an equivalent pseudo-code for a computationally-efficient implementation. The most expensive computations —the sparse projection at line 10 and the gradient at line 12— are only required when the support has never been explored before. Also, the best sparse vector can be memorized on the fly (line 16). Hence, the remaining operations, that are performed at each iteration, have a low computational cost: support extraction (line 8), search for a previous, identical support (line 9) and STE update (line 21). This computationally-efficient version of SEA has a larger spatial complexity due to the memorization of all the supports and gradients seen along the iterations. However, this overhead is limited since 1/ for each explored support, only two vectors are memorized, one of them being sparse; and 2/ the number of explored supports is generally much lower than the number of iterations. For instance, in the deconvolution experiment, on average, less than 1000 vectors of size 500 (including 500 k-sparse vectors) are stored during the running time of SEA.

In addition, solving the (unconstrained) linear system $A_{S^t}^T A_{S^t} x_{S^t} = A_{S^t}^T y$ can also be performed efficiently. While the pseudo-inverse is a convenient notation at line 10, the solution may be obtained more efficiently, e.g. in $O(k^2n)$ to compute $A_{S^t}^T A_{S^t}$ and $A_{S^t}^T y$, and apply the conjugate gradient algorithm. This complexity is a worst-case scenario so in practice, the solution is generally obtained more quickly.
<caption>**Algorithm 5** Support Exploration Algorithm: efficient implementation</caption>

1: **Input:** noisy observation $y$, sampling matrix $A$, sparsity $k$, step size $\eta$
2: **Output:** sparse vector $x$
3: Initialize $\mathcal{X}^0$
4: $F^{\text{BEST}} \leftarrow +\infty$
5: $t \leftarrow 0$
6: $\mathbf{S} \leftarrow \{\}, \mathbf{g} \leftarrow \{\}$
7: **repeat**
8: &nbsp;&nbsp;&nbsp;&nbsp;$S \leftarrow \text{largest}_k(\mathcal{X}^t)$
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;{Compute sparse vector and gradient only for unseen supports}
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;9: **if** $S \notin \mathbf{S}$ **then**
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$x_i^S \leftarrow 0$ for $i \in \bar{S}$
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$x_S^S \leftarrow A_S^\dagger y$
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\text{loss}^S \leftarrow \frac{1}{2}\|Ax^S - y\|_2^2$
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$g^S \leftarrow \eta A^T (Ax^S - y)$
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;{Memorize support and gradient}
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{S} \leftarrow \mathbf{S} \cup \{S\}, \mathbf{g} \leftarrow \mathbf{g} \cup \{g^S\}$
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;{Memorize best iterate}
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;14: **if** $\text{loss}^S < F^{\text{BEST}}$ **then**
19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$F^{\text{BEST}} \leftarrow \text{loss}^S$
20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$x^{\text{BEST}} \leftarrow x^S$
21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\text{end if}$
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\text{else}$
23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Retrieve $g^S$ in $\mathbf{g}$
24: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**return** $x^{\text{BEST}}$
<caption>**Algorithm 5** Support Exploration Algorithm: efficient implementation</caption>

1: **Input:** noisy observation $y$, sampling matrix $A$, sparsity $k$, step size $\eta$
2: **Output:** sparse vector $x$
3: Initialize $\mathcal{X}^0$
4: $F^{\text{BEST}} \leftarrow +\infty$
5: $t \leftarrow 0$
6: $\mathbf{S} \leftarrow \{\}, \mathbf{g} \leftarrow \{\}$
7: **repeat**
8:   $S \leftarrow \text{largest}_k(\mathcal{X}^t)$
9:   \{Compute sparse vector and gradient only for unseen supports\}
10:   9:   **if** $S \notin \mathbf{S}$ **then**
11:     $x_i^S \leftarrow 0$ for $i \in \bar{S}$
12:     $x_S^S \leftarrow A_S^\dagger y$
13:     $\text{loss}^S \leftarrow \frac{1}{2} \|Ax^S - y\|_2^2$
14:     $g^S \leftarrow \eta A^T (Ax^S - y)$
15:     \{Memorize support and gradient\}
16:     $\mathbf{S} \leftarrow \mathbf{S} \cup \{S\}$, $\mathbf{g} \leftarrow \mathbf{g} \cup \{g^S\}$
17:     \{Memorize best iterate\}
18:     14:     **if** $\text{loss}^S < F^{\text{BEST}}$ **then**
19:       $F^{\text{BEST}} \leftarrow \text{loss}^S$
20:       $x^{\text{BEST}} \leftarrow x^S$
21:     17:     **end if**
22:     18:     **else**
23:     19:     Retrieve $g^S$ in $\mathbf{g}$
24:     20:     **end if**
21:     \{Update support exploration variable\}
22:     $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - g^S$
23:     $t \leftarrow t + 1$
24:     **until** halting criterion is true
25:     **return** $x^{\text{BEST}}$
# C Proofs and Complements of the Theoretical Analysis

The proof of Theorem 4.1 relies on the fact that when the hypotheses of the theorem are satisfied, the trajectory $(\mathcal{X}^t)_{t \in \mathbb{N}}$ is close to the trajectory of an algorithm that has access to an oracle update. The appendix first contains a description and an analysis of this algorithm in Appendix C.1. Then, in Appendix C.2, we analyze how much the STE-update can deviate from the Oracle Update Rule so that the true support $S^*$ is still recovered. Finally, we prove Theorem 4.1 in Appendix C.4. We prove Corollary 4.2 in Appendix C.5 and conclude with comments on the conditions $(RC_{RIP})$ and $(RC_{SRIP})$ in Appendix C.6.

## C.1 Support Exploration Algorithm Using the Oracle Update Rule

The theoretical analysis of SEA and the understanding of the underlying behavior of the algorithm rely on the introduction of an oracle case where the true solution $x^*$ and its support $S^*$ are known by the algorithm. In that case, at iteration $t$, we can use the oracle update rule $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - u^t$, using the direction $u^t$ defined for any index $i \in [1, n]$ by

$$ u_i^t = \begin{cases} -\eta x_i^* & i \in S^* \cap \overline{S^t} \\ 0 & i \in \overline{S^*} \cup S^t, \end{cases} \qquad (15) $$

where $S^t = \text{largest}_k(\mathcal{X}^t)$ contains the indices of the $k$ largest absolute entries in $\mathcal{X}^t$ and $\eta > 0$ is an arbitrary step size. The resulting pseudo-code is given by Algorithm 6 and we show the important supports in Figure 5.

<caption>**Algorithm 6** Support Exploration Algorithm using the Oracle Update Rule</caption>

1: **Input:** true solution $x^*$, true support $S^*$, sparsity $k$, step size $\eta$, noisy observation $y$, sampling matrix $A$
2: **Output:** sparse vector $x$
3: Initialize $\mathcal{X}^0$
4: $t \leftarrow 0$
5: **repeat**
6: &nbsp;&nbsp;&nbsp;&nbsp;$S^t \leftarrow \text{largest}_k(\mathcal{X}^t)$
7: &nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - u^t$
8: &nbsp;&nbsp;&nbsp;&nbsp;$t \leftarrow t + 1$
9: &nbsp;&nbsp;&nbsp;&nbsp;**until** $u^{t-1} = 0$
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \begin{cases} x_i \leftarrow 0 & \text{for } i \in \overline{S^{t-1}} \\ x_{S^{t-1}} \leftarrow A_{S^{t-1}}^\dagger y \end{cases} $$
11: **return** $x$

<figure><img src="image_7.png" /><figcaption>Figure 5: Visual representation of the main sets of indices encountered in the article.</figcaption></figure>

Notice $u_i^t$ is non-zero for indices $i$ from the true support $S^*$ but for which $|\mathcal{X}_i^t|$ is too small to be selected in $S^t$ at line 6. Whatever the initial content of $\mathcal{X}^0$, the oracle update rule always adds the same increment to $\mathcal{X}_i^t$, for $i \in S^* \cap \overline{S^t}$. This guarantees that, at some subsequent iteration $t' \ge t$, the true support $S^*$ is recovered among the $k$ largest absolute entries in $\mathcal{X}^{t'}$, i.e., $S^* \subseteq S^{t'}$, the intersection is empty, $u^{t'} = 0$ and Algorithm 6 stops.

In the following theorem, we formalize this statement and give an upper bound on the number of iterations required by the support exploration algorithm using the oracle update rule.
# C Proofs and Complements of the Theoretical Analysis

The proof of Theorem 4.1 relies on the fact that when the hypotheses of the theorem are satisfied, the trajectory $(\mathcal{X}^t)_{t \in \mathbb{N}}$ is close to the trajectory of an algorithm that has access to an oracle update. The appendix first contains a description and an analysis of this algorithm in Appendix C.1. Then, in Appendix C.2, we analyze how much the STE-update can deviate from the Oracle Update Rule so that the true support $S^*$ is still recovered. Finally, we prove Theorem 4.1 in Appendix C.4. We prove Corollary 4.2 in Appendix C.5 and conclude with comments on the conditions $(RC_{RIP})$ and $(RC_{SRIP})$ in Appendix C.6.

## C.1 Support Exploration Algorithm Using the Oracle Update Rule

The theoretical analysis of SEA and the understanding of the underlying behavior of the algorithm rely on the introduction of an oracle case where the true solution $x^*$ and its support $S^*$ are known by the algorithm. In that case, at iteration $t$, we can use the oracle update rule $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - u^t$, using the direction $u^t$ defined for any index $i \in [1, n]$ by

$$u_i^t = \begin{cases} -\eta x_i^* & i \in S^* \cap \overline{S^t} \\ 0 & i \in \overline{S^*} \cup S^t, \end{cases} \quad (15)$$

where $S^t = \text{largest}_k(\mathcal{X}^t)$ contains the indices of the $k$ largest absolute entries in $\mathcal{X}^t$ and $\eta > 0$ is an arbitrary step size. The resulting pseudo-code is given by Algorithm 6 and we show the important supports in Figure 5.

<caption>**Algorithm 6** Support Exploration Algorithm using the Oracle Update Rule</caption>

1: **Input:** true solution $x^*$, true support $S^*$, sparsity $k$, step size $\eta$, noisy observation $y$, sampling matrix $A$
2: **Output:** sparse vector $x$
3: Initialize $\mathcal{X}^0$
4: $t \leftarrow 0$
5: **repeat**
6: &nbsp;&nbsp;&nbsp;&nbsp;$S^t \leftarrow \text{largest}_k(\mathcal{X}^t)$
7: &nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - u^t$
8: &nbsp;&nbsp;&nbsp;&nbsp;$t \leftarrow t + 1$
9: &nbsp;&nbsp;&nbsp;&nbsp;**until** $u^{t-1} = 0$
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$x_i \leftarrow 0 \text{ for } i \in \overline{S^{t-1}}$$
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$x_{S^{t-1}} \leftarrow A_{S^{t-1}}^\dagger y$$
12: **return** $x$

<figure><img src="image_8.png" /><figcaption>Figure 5: Visual representation of the main sets of indices encountered in the article.</figcaption></figure>

Notice $u_i^t$ is non-zero for indices $i$ from the true support $S^*$ but for which $|\mathcal{X}_i^t|$ is too small to be selected in $S^t$ at line 6. Whatever the initial content of $\mathcal{X}^0$, the oracle update rule always adds the same increment to $\mathcal{X}_i^t$, for $i \in S^* \cap \overline{S^t}$. This guarantees that, at some subsequent iteration $t' \ge t$, the true support $S^*$ is recovered among the $k$ largest absolute entries in $\mathcal{X}^{t'}$, i.e., $S^* \subseteq S^{t'}$, the intersection is empty, $u^{t'} = 0$ and Algorithm 6 stops.

In the following theorem, we formalize this statement and give an upper bound on the number of iterations required by the support exploration algorithm using the oracle update rule.
**Theorem C.1** (Recovery - Oracle Update Rule). For all matrices A, error vectors e, initializations $\mathcal{X}^0$ and for all $\eta > 0$, there exists

$$t_s \le T_{max}^{oracle} = k \left( 1 + \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta \min_{i \in S^*} |x_i^*|} \right)$$

such that $S^* \subseteq S^{t_s}$, where $S^t$ is defined in Algorithm 6 line 6.

Moreover, $u^{t_s} = 0$ and Algorithm 6 returns $\underset{x \in \mathbb{R}^n}{\operatorname{argmin}} \|Ax - y\|_2^2$. $\operatorname{supp}(x) \subseteq S^{t_s}$

### C.1.1 Proof of Theorem C.1

We denote, for all $t \in \mathbb{N}^*$ and all $i \in [1, n]$, and $S^t$ defined in Algorithm 6, line 6

$$c_i^t = |\{t' \in [0, t-1] : i \in S^* \cap \overline{S^{t'}}\}|. \qquad (16)$$

We extend the definition to $t = 0$ and set, for all $i \in [1, n]$, $c_i^0 = 0$.

We can prove by induction on $t$ that, given the definition of $\mathcal{X}^t$ in Algorithm 6 and $u^t$ in (15), for all $t \in \mathbb{N}$,

$$\mathcal{X}^t = \mathcal{X}^0 + \eta c^t \odot x^*, \qquad (17)$$

where $\odot$ denotes the Hadamard product.

The following lemma states that if $c_i^t$ is large then $i$ is always selected by Algorithm 6.

**Lemma C.2.** For all $i \in S^*$ and all $t \in \mathbb{N}^*$

$$\text{if } c_i^t > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|} \quad \text{then} \quad \forall t' \ge t, \quad i \in S^{t'}. $$

*Proof.* Let $i \in S^*$ and $t \in \mathbb{N}^*$ be such that $c_i^t > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|}$. Consider $t' \ge t$. Since $t \mapsto c_i^t$ is non-decreasing, we have

$$c_i^{t'} \ge c_i^t > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|}.$$

Therefore, for all $j \in \overline{S^*}$,

$$|\mathcal{X}_i^{t'}| = |\mathcal{X}_i^0 + \eta c_i^{t'} x_i^*| \ge |\eta c_i^{t'} x_i^*| - |\mathcal{X}_i^0| > 2\|\mathcal{X}^0\|_{\infty} - |\mathcal{X}_i^0| \ge \|\mathcal{X}^0\|_{\infty} \ge |\mathcal{X}_j^0| = |\mathcal{X}_j^{t'}|.$$

Therefore $|\mathcal{X}_i^{t'}|$ is larger than at least $n-k$ elements of $\{|\mathcal{X}_j^{t'}| : j \in [1, n]\}$. Said differently, $i \in \text{largest}_k(\mathcal{X}^{t'}) = S^{t'}$.

This concludes the proof of Lemma C.2.

<figure><img src="image_8.png" /></figure>

This leads to the following upper bound.

**Lemma C.3.** For all $i \in S^*$ and all $t \in \mathbb{N}^*$

$$c_i^t \le \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|} + 1.$$

*Proof.* If Lemma C.3 is false, there exists $i \in S^*$ and $t \in \mathbb{N}^*$ such that $c_i^t > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|} + 1$. We denote

$$t' = \min\{t \in \mathbb{N}^* : c_i^t > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|} + 1\}.$$

We have $c_i^{t'} > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|} + 1 \ge 1$ and therefore $t' > 1$. As a consequence, $t' - 1 \in \mathbb{N}^*$ and $c_i^{t'-1}$ is defined by (16). Because of the definitions of $t'$ and $c_i^{t'}$, we must have $c_i^{t'} = c_i^{t'-1} + 1$. Therefore, $c_i^{t'-1} > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|}$. Since $i \in S^*$ and $t' - 1 \in \mathbb{N}^*$, using Lemma C.2, we conclude that $i \in S^{t'-1}$ and, using the definition of $c_i^{t'}$ in (16), that $c_i^{t'} = c_i^{t'-1}$.

This is impossible and we conclude that Lemma C.3 holds.

<figure><img src="image_9.png" /></figure>
**Theorem C.1** (Recovery - Oracle Update Rule). For all matrices A, error vectors e, initializations $\mathcal{X}^0$ and for all $\eta > 0$, there exists

$$t_s \le T_{max}^{oracle} = k \left( 1 + \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta \min_{i \in S^*} |x_i^*|} \right)$$

such that $S^* \subseteq S^{t_s}$, where $S^t$ is defined in Algorithm 6 line 6.

Moreover, $u^{t_s} = 0$ and Algorithm 6 returns $\underset{x \in \mathbb{R}^n}{\operatorname{argmin}} \|Ax - y\|_2^2$. $\operatorname{supp}(x) \subseteq S^{t_s}$

### C.1.1 Proof of Theorem C.1

We denote, for all $t \in \mathbb{N}^*$ and all $i \in [1, n]$, and $S^t$ defined in Algorithm 6, line 6

$$c_i^t = |\{t' \in [0, t-1] : i \in S^* \cap \overline{S^{t'}}\}|. \qquad (16)$$

We extend the definition to $t = 0$ and set, for all $i \in [1, n]$, $c_i^0 = 0$.

We can prove by induction on $t$ that, given the definition of $\mathcal{X}^t$ in Algorithm 6 and $u^t$ in (15), for all $t \in \mathbb{N}$,

$$\mathcal{X}^t = \mathcal{X}^0 + \eta c^t \odot x^*, \qquad (17)$$

where $\odot$ denotes the Hadamard product.

The following lemma states that if $c_i^t$ is large then $i$ is always selected by Algorithm 6.

**Lemma C.2.** For all $i \in S^*$ and all $t \in \mathbb{N}^*$

$$\text{if } c_i^t > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|} \quad \text{then} \quad \forall t' \ge t, \quad i \in S^{t'}. $$

*Proof.* Let $i \in S^*$ and $t \in \mathbb{N}^*$ be such that $c_i^t > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|}$. Consider $t' \ge t$. Since $t \mapsto c_i^t$ is non-decreasing, we have

$$c_i^{t'} \ge c_i^t > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|}.$$

Therefore, for all $j \in \overline{S^*}$,

$$|\mathcal{X}_i^{t'}| = |\mathcal{X}_i^0 + \eta c_i^{t'} x_i^*| \ge |\eta c_i^{t'} x_i^*| - |\mathcal{X}_i^0| > 2\|\mathcal{X}^0\|_{\infty} - |\mathcal{X}_i^0| \ge \|\mathcal{X}^0\|_{\infty} \ge |\mathcal{X}_j^0| = |\mathcal{X}_j^{t'}|.$$

Therefore $|\mathcal{X}_i^{t'}|$ is larger than at least $n-k$ elements of $\{|\mathcal{X}_j^{t'}| : j \in [1, n]\}$. Said differently, $i \in \text{largest}_k(\mathcal{X}^{t'}) = S^{t'}$.

This concludes the proof of Lemma C.2. □

This leads to the following upper bound.

**Lemma C.3.** For all $i \in S^*$ and all $t \in \mathbb{N}^*$

$$c_i^t \le \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|} + 1.$$

*Proof.* If Lemma C.3 is false, there exists $i \in S^*$ and $t \in \mathbb{N}^*$ such that $c_i^t > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|} + 1$. We denote

$$t' = \min\{t \in \mathbb{N}^* : c_i^t > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|} + 1\}.$$

We have $c_i^{t'} > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|} + 1 \ge 1$ and therefore $t' > 1$. As a consequence, $t' - 1 \in \mathbb{N}^*$ and $c_i^{t'-1}$ is defined by (16). Because of the definitions of $t'$ and $c_i^{t'}$, we must have $c_i^{t'} = c_i^{t'-1} + 1$. Therefore, $c_i^{t'-1} > \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|}$. Since $i \in S^*$ and $t' - 1 \in \mathbb{N}^*$, using Lemma C.2, we conclude that $i \in S^{t'-1}$ and, using the definition of $c_i^{t'}$ in (16), that $c_i^{t'} = c_i^{t'-1}$.

This is impossible and we conclude that Lemma C.3 holds. □
**Proof of Theorem C.1:** We denote

$$t_s = \min\{t \in \mathbb{N} : S^* \subseteq S^t\}. \qquad (18)$$

By convention, if for all $t \in \mathbb{N}$, $S^* \not\subseteq S^t$, we set $t_s = +\infty$. The first statement of Theorem C.1 is obvious if $t_s = 0$. We assume below that $t_s \ge 1$.

Consider $t \in [0, t_s - 1]$, using of the definition of $t_s$ in (18), there exists $i \in S^* \cap \overline{S^t}$. Using the definition of $c_i^{t+1}$ in (16), we obtain $c_i^{t+1} = c_i^t + 1$. Since for all $j \in [1, n]$, $t \mapsto c_j^t$ is non-decreasing, we conclude that

$$\text{for all } t \in [0, t_s - 1], \quad \sum_{i \in S^*} c_i^{t+1} \ge \sum_{i \in S^*} c_i^t + 1.$$

We therefore obtain

$$\begin{aligned} \sum_{i \in S^*} c_i^{t_s} &= \sum_{i \in S^*} \left( \sum_{i=0}^{t_s-1} (c_i^{t+1} - c_i^t) \right) \\ &= \sum_{i=0}^{t_s-1} \left( \left( \sum_{i \in S^*} c_i^{t+1} \right) - \left( \sum_{i \in S^*} c_i^t \right) \right) \\ &\ge \sum_{i=0}^{t_s-1} 1 = t_s. \end{aligned}$$

Using Lemma C.3, we obtain

$$t_s \le \sum_{i \in S^*} c_i^{t_s} \le \sum_{i \in S^*} \left( \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|} + 1 \right) \le k \left( 1 + \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta \min_{i \in S^*} |x_i^*|} \right).$$

To conclude the proof, we simply remark that, since $S^* \subseteq S^{t_s}$, by definition of $u^{t_s}$ in (15),

$$u^{t_s} = 0.$$

<figure><img src="image_10.png" /></figure>

Since $x^*$ and $S^*$ are not available in practice, we replace in Algorithm 6 the oracle update $u^t$ by the surrogate $\eta A^T(Ax^t - y)$ (line 8). The choice of this surrogate is an application of STE and is natural. For instance, one can show that $u^t = \eta A^T(Ax^t - y)$ in the simple case where A has orthonormal columns and the observation is noiseless (see Corollary C.7 and Lemma C.8 in Appendix C.3). In the general setting, SEA can be interpreted as a noisy version of the support exploration algorithm using the oracle update. Theorem 4.1 and its proof in Appendix C are based on the fact that $u^t - \eta A^T(Ax^t - y)$ is small, under suitable hypotheses on $x^*$ and the RIP constants of A.

## C.2 If the STE-Update is Sufficiently Close to the Oracle Update, SEA Visits $S^*$

To prove Theorem 4.1, we first provide a general recovery theorem here. The theorem states that if the discrepancy between the Oracle update and the STE update is sufficiently small, SEA visits $S^*$.

To do so, we define, for all $t \in \mathbb{N}$, the gradient noise: $b^t \in \mathbb{R}^n$ as

$$b^t = u^t - \eta A^T(Ax^t - y). \qquad (19)$$

We define the maximal gradient noise norm

$$\varepsilon = \sup_{t \in \mathbb{N}} \|b^t\|_{\infty} \in \mathbb{R}. \qquad (20)$$

We define the Recovery Condition (RC) as

$$\varepsilon < \frac{\eta}{2k} \min_{i \in S^*} |x_i^*|. \qquad (RC)$$

**Theorem C.4 (Recovery - General case).** If (RC) holds, then for all initializations $\mathcal{X}^0$ and all $\eta > 0$, there exists $t_s \le T_{max}$ such that $S^* \subseteq S^{t_s}$, where $S^t$ is defined in Algorithm 1 line 6 and

$$T_{max} = \frac{2k\|\mathcal{X}^0\|_{\infty} + (k+1)\eta \min_{i \in S^*} |x_i^*|}{\eta \min_{i \in S^*} |x_i^*| - 2k\varepsilon}. \qquad (21)$$
**Proof of Theorem C.1:** We denote

$$t_s = \min\{t \in \mathbb{N} : S^* \subseteq S^t\}. \qquad (18)$$

By convention, if for all $t \in \mathbb{N}$, $S^* \not\subseteq S^t$, we set $t_s = +\infty$. The first statement of Theorem C.1 is obvious if $t_s = 0$. We assume below that $t_s \ge 1$.

Consider $t \in [0, t_s - 1]$, using of the definition of $t_s$ in (18), there exists $i \in S^* \cap \overline{S^t}$. Using the definition of $c_i^{t+1}$ in (16), we obtain $c_i^{t+1} = c_i^t + 1$. Since for all $j \in [1, n]$, $t \mapsto c_j^t$ is non-decreasing, we conclude that

$$\text{for all } t \in [0, t_s - 1], \quad \sum_{i \in S^*} c_i^{t+1} \ge \sum_{i \in S^*} c_i^t + 1.$$

We therefore obtain

$$\begin{aligned} \sum_{i \in S^*} c_i^{t_s} &= \sum_{i \in S^*} \left( \sum_{i=0}^{t_s-1} (c_i^{t+1} - c_i^t) \right) \\ &= \sum_{i=0}^{t_s-1} \left( \left( \sum_{i \in S^*} c_i^{t+1} \right) - \left( \sum_{i \in S^*} c_i^t \right) \right) \\ &\ge \sum_{i=0}^{t_s-1} 1 = t_s. \end{aligned}$$

Using Lemma C.3, we obtain

$$t_s \le \sum_{i \in S^*} c_i^{t_s} \le \sum_{i \in S^*} \left( \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta|x_i^*|} + 1 \right) \le k \left( 1 + \frac{2\|\mathcal{X}^0\|_{\infty}}{\eta \min_{i \in S^*} |x_i^*|} \right).$$

To conclude the proof, we simply remark that, since $S^* \subseteq S^{t_s}$, by definition of $u^{t_s}$ in (15),

$$u^{t_s} = 0.$$

<figure><img src="image_9.png" /></figure>

Since $x^*$ and $S^*$ are not available in practice, we replace in Algorithm 6 the oracle update $u^t$ by the surrogate $\eta A^T(Ax^t - y)$ (line 8). The choice of this surrogate is an application of STE and is natural. For instance, one can show that $u^t = \eta A^T(Ax^t - y)$ in the simple case where A has orthonormal columns and the observation is noiseless (see Corollary C.7 and Lemma C.8 in Appendix C.3). In the general setting, SEA can be interpreted as a noisy version of the support exploration algorithm using the oracle update. Theorem 4.1 and its proof in Appendix C are based on the fact that $u^t - \eta A^T(Ax^t - y)$ is small, under suitable hypotheses on $x^*$ and the RIP constants of A.

## C.2 If the STE-Update is Sufficiently Close to the Oracle Update, SEA Visits $S^*$

To prove Theorem 4.1, we first provide a general recovery theorem here. The theorem states that if the discrepancy between the Oracle update and the STE update is sufficiently small, SEA visits $S^*$.

To do so, we define, for all $t \in \mathbb{N}$, the gradient noise: $b^t \in \mathbb{R}^n$ as

$$b^t = u^t - \eta A^T(Ax^t - y). \qquad (19)$$

We define the maximal gradient noise norm

$$\varepsilon = \sup_{t \in \mathbb{N}} \|b^t\|_{\infty} \in \mathbb{R}. \qquad (20)$$

We define the Recovery Condition (RC) as

$$\varepsilon < \frac{\eta}{2k} \min_{i \in S^*} |x_i^*|. \qquad (RC)$$

**Theorem C.4 (Recovery - General case).** If (RC) holds, then for all initializations $\mathcal{X}^0$ and all $\eta > 0$, there exists $t_s \le T_{max}$ such that $S^* \subseteq S^{t_s}$, where $S^t$ is defined in Algorithm 1 line 6 and

$$T_{max} = \frac{2k\|\mathcal{X}^0\|_{\infty} + (k+1)\eta \min_{i \in S^*} |x_i^*|}{\eta \min_{i \in S^*} |x_i^*| - 2k\varepsilon}. \qquad (21)$$
The proof is in Appendix C.2.1, right after the comments below.

The main interest of Theorem C.4 is to formalize quantitatively that, when $u^t - \eta A^T(Ax^t - y)$ is sufficiently small, SEA visits the correct support. However, the condition (RC) is difficult to use and interpret since it involves both $A$, $x^*$, and all the sparse iterates $x^t$. This is why we provide in Corollary C.7, Theorem 4.1 and Corollary 4.2 sufficient conditions on $A$, $e$ and $x^*$ guaranteeing that (RC) holds.

The conclusion of Theorem C.4 is that the iterative process of SEA visits the correct support at some iteration $t$. We have in general no guarantee that this time-step $t$ is equal to $t_{BEST}$. We are however guaranteed that SEA returns a sparse solution such that $\|Ax^{t_{BEST}} - y\|_2 \le \|Ax^{t_s} - y\|_2 \le \|Ax^* - y\|_2$. This does not give a guarantee on the support recovery but on the reconstruction error. In machine learning, this upper bound can be used to derive an upper bound of the risk.

Concerning the value of $T_{max}$, a quick analysis shows that $T_{max}$ increases with $\varepsilon$, when (RC) holds. In other words, the number of iterations required by the algorithm to provide the correct solution increases with the discrepancy between $u^t$ and $\eta A^T(Ax^t - y)$.

The initializations $\mathcal{X}^0 \ne 0$ have an apparent negative impact on the number of iterations required in the worst case. This is because in the worst-case $\mathcal{X}^0$ would be poorly chosen and SEA needs iterations to correct this poor choice.

Concerning $\eta$, notice that, since $u^t$ is proportional to $\eta > 0$, $\varepsilon$ is proportional to $\eta > 0$ and therefore (RC) is independent of $\eta$. When possible, any $\eta$ permits the recovery of $S^*$. The only influence of $\eta$ is on $T_{max}$. In this regard, since $\varepsilon$ is proportional to $\eta > 0$, the denominator of (21) is proportional to $\eta$. In the numerator, we see that the larger $\eta$ is, the faster SEA will override the initialization $\mathcal{X}^0$. The choice of $\eta$ is very much related to the question of the quality of the initialization discussed above.

### C.2.1 Proof of Theorem C.4

To prove Theorem C.4, we need to adapt a closed formula for the exploratory variable $\mathcal{X}^t$ already encountered in the proof of Theorem C.1. Then, we will study the properties of this closed formula through the counting vector $c^t$ in Appendix C.2.2. to find a sufficient condition of support recovery. Then we prove Theorem C.4 in Appendix C.2.3.

### C.2.2 Preliminaries

We remind Figure 5 on which the mains supports are represented and we remind, for each iteration $t \in \mathbb{N}$ and $i \in [1, n]$, the Oracle Update already defined in (15)

$$u_i^t = \begin{cases} -\eta x_i^* & i \in S^* \cap \overline{S}^t \\ 0 & i \in \overline{S}^* \cup S^t. \end{cases}$$

We also remind the gradient noise, already defined in (19), $b^t = u^t - \eta A^T(Ax^t - y)$.

We first remark that, for any $i \in S^t$,

$$b_i^t = 0. \qquad (22)$$

To prove this equality, we remark that, for all $i \in S^t$, $u_i^t = 0$ and prove that $(A^T(Ax^t - y))_i = 0$. Indeed, the latter holds because $i \in S^t$ and $x_{S^t}^t = A_{S^t}^\dagger y$ (see Algorithm 1, line 7). As is well-known for the Moore-Penrose inverse, $A_{S^t} A_{S^t}^\dagger$ is the orthogonal projector onto $\text{colspan}(A_{S^t})$. Therefore, $A_{S^t} A_{S^t}^\dagger y - y$ is orthogonal to $\text{colspan}(A_{S^t})$ and for all $x' \in \mathbb{R}^k$, $0 = \langle A_{S^t} A_{S^t}^\dagger y - y, A_{S^t} x' \rangle = \langle A_{S^t}^T (A_{S^t} A_{S^t}^\dagger y - y), x' \rangle$. Therefore, $A_{S^t}^T (A_{S^t} A_{S^t}^\dagger y - y) = 0$. Since, $x_i^t = 0$ for all $i \in \overline{S}^t$, we also have $A_{S^t} A_{S^t}^\dagger y = A_{S^t} x_{S^t}^t = Ax^t$ and we deduce that for all $i \in S^t$, $(A^T(Ax^t - y))_i = (A_{S^t}^T (A_{S^t} A_{S^t}^\dagger y - y))_j = 0$, where the line $j \in [1, k]$ of $A_{S^t}^T$ corresponds to the line $i \in S^t$ of $A^T$. This concludes the proof of (22).

As a consequence of the definition of $b^t$ and SEA, line 8, for any $t \in \mathbb{N}$,

$$\mathcal{X}^{t+1} = \mathcal{X}^t + b^t - u^t. \qquad (23)$$

The gradient noise $b^t$ is the error preventing the gradient from being in the direction of the oracle update $u^t$. At each iteration, this error is accumulating in $\mathcal{X}^t$. With $\beta^0 = 0$, for any $t \in \mathbb{N}^*$, we define this accumulated error by

$$\beta^t = \sum_{t'=0}^{t-1} b_{t'} \in \mathbb{R}^n. \qquad (24)$$

As already done in the proof of Theorem C.1 for the support sequence defined in Algorithm 6, we define counting vectors. However, this time they are defined for the sequence defined in Algorithm 1. We keep
The proof is in Appendix C.2.1, right after the comments below.

The main interest of Theorem C.4 is to formalize quantitatively that, when $u^t - \eta A^T(Ax^t - y)$ is sufficiently small, SEA visits the correct support. However, the condition (RC) is difficult to use and interpret since it involves both $A$, $x^*$, and all the sparse iterates $x^t$. This is why we provide in Corollary C.7, Theorem 4.1 and Corollary 4.2 sufficient conditions on $A$, $e$ and $x^*$ guaranteeing that (RC) holds.

The conclusion of Theorem C.4 is that the iterative process of SEA visits the correct support at some iteration $t$. We have in general no guarantee that this time-step $t$ is equal to $t_{BEST}$. We are however guaranteed that SEA returns a sparse solution such that $\|Ax^{t_{BEST}} - y\|_2 \le \|Ax^{t_s} - y\|_2 \le \|Ax^* - y\|_2$. This does not give a guarantee on the support recovery but on the reconstruction error. In machine learning, this upper bound can be used to derive an upper bound of the risk.

Concerning the value of $T_{max}$, a quick analysis shows that $T_{max}$ increases with $\varepsilon$, when (RC) holds. In other words, the number of iterations required by the algorithm to provide the correct solution increases with the discrepancy between $u^t$ and $\eta A^T(Ax^t - y)$.

The initializations $\mathcal{X}^0 \ne 0$ have an apparent negative impact on the number of iterations required in the worst case. This is because in the worst-case $\mathcal{X}^0$ would be poorly chosen and SEA needs iterations to correct this poor choice.

Concerning $\eta$, notice that, since $u^t$ is proportional to $\eta > 0$, $\varepsilon$ is proportional to $\eta > 0$ and therefore (RC) is independent of $\eta$. When possible, any $\eta$ permits the recovery of $S^*$. The only influence of $\eta$ is on $T_{max}$. In this regard, since $\varepsilon$ is proportional to $\eta > 0$, the denominator of (21) is proportional to $\eta$. In the numerator, we see that the larger $\eta$ is, the faster SEA will override the initialization $\mathcal{X}^0$. The choice of $\eta$ is very much related to the question of the quality of the initialization discussed above.

### C.2.1 Proof of Theorem C.4

To prove Theorem C.4, we need to adapt a closed formula for the exploratory variable $\mathcal{X}^t$ already encountered in the proof of Theorem C.1. Then, we will study the properties of this closed formula through the counting vector $c^t$ in Appendix C.2.2. to find a sufficient condition of support recovery. Then we prove Theorem C.4 in Appendix C.2.3.

### C.2.2 Preliminaries

We remind Figure 5 on which the mains supports are represented and we remind, for each iteration $t \in \mathbb{N}$ and $i \in [1, n]$, the Oracle Update already defined in (15)

$$u_i^t = \begin{cases} -\eta x_i^* & i \in S^* \cap \overline{S}^t \\ 0 & i \in \overline{S}^* \cup S^t. \end{cases}$$

We also remind the gradient noise, already defined in (19), $b^t = u^t - \eta A^T(Ax^t - y)$.

We first remark that, for any $i \in S^t$,

$$b_i^t = 0. \qquad (22)$$

To prove this equality, we remark that, for all $i \in S^t$, $u_i^t = 0$ and prove that $(A^T(Ax^t - y))_i = 0$. Indeed, the latter holds because $i \in S^t$ and $x_{S^t}^t = A_{S^t}^\dagger y$ (see Algorithm 1, line 7). As is well-known for the Moore-Penrose inverse, $A_{S^t} A_{S^t}^\dagger$ is the orthogonal projector onto $\text{colspan}(A_{S^t})$. Therefore, $A_{S^t} A_{S^t}^\dagger y - y$ is orthogonal to $\text{colspan}(A_{S^t})$ and for all $x' \in \mathbb{R}^k$, $0 = \langle A_{S^t} A_{S^t}^\dagger y - y, A_{S^t} x' \rangle = \langle A_{S^t}^T (A_{S^t} A_{S^t}^\dagger y - y), x' \rangle$. Therefore, $A_{S^t}^T (A_{S^t} A_{S^t}^\dagger y - y) = 0$. Since, $x_i^t = 0$ for all $i \in \overline{S}^t$, we also have $A_{S^t} A_{S^t}^\dagger y = A_{S^t} x_{S^t}^t = Ax^t$ and we deduce that for all $i \in S^t$, $(A^T(Ax^t - y))_i = (A_{S^t}^T (A_{S^t} A_{S^t}^\dagger y - y))_j = 0$, where the line $j \in [1, k]$ of $A_{S^t}^T$ corresponds to the line $i \in S^t$ of $A^T$. This concludes the proof of (22).

As a consequence of the definition of $b^t$ and SEA, line 8, for any $t \in \mathbb{N}$,

$$\mathcal{X}^{t+1} = \mathcal{X}^t + b^t - u^t. \qquad (23)$$

The gradient noise $b^t$ is the error preventing the gradient from being in the direction of the oracle update $u^t$. At each iteration, this error is accumulating in $\mathcal{X}^t$. With $\beta^0 = 0$, for any $t \in \mathbb{N}^*$, we define this accumulated error by

$$\beta^t = \sum_{t'=0}^{t-1} b_{t'} \in \mathbb{R}^n. \qquad (24)$$

As already done in the proof of Theorem C.1 for the support sequence defined in Algorithm 6, we define counting vectors. However, this time they are defined for the sequence defined in Algorithm 1. We keep
the same notations for simplicity. We set $c^0 = 0$, for any $t \in \mathbb{N}^*$ and $i \in [1, n]$, we also define the counting vector by

$$c_i^t = |\{t' \in [0, t-1] : i \in S^* \cap \overline{S^{t'}}\}|. \qquad (25)$$

We will use the recursive formula for $c^t$: For any $t \in \mathbb{N}$, $i \in [1, n]$

$$c_i^{t+1} = \begin{cases} c_i^t + 1 & \text{if } i \in S^* \cap \overline{S^t} \\ c_i^t & \text{if } i \in \overline{S^*} \cup S^t. \end{cases} \qquad (26)$$

For any $i \in [1, n]$, the sequence $(c_i^t)_{t \in \mathbb{N}}$ is non-decreasing.

Using (23), (24) and (25), we can establish by induction on $t$ that for any $t \in \mathbb{N}$,

$$\mathcal{X}^t = \mathcal{X}^0 + \eta c^t \odot x^* + \beta^t, \qquad (27)$$

where $\odot$ denotes the Hadamard product. This generalizes (17) to the noisy setting.

As can be seen from (27), the error accumulation $\beta^t$ is responsible for the exploration in the wrong directions. While $c^t \odot x^*$ encourages exploration in the direction of the missed components of $x^*$. Below, we provide important properties of $(c^t)_{t \in \mathbb{N}}$.

At each iteration of SEA, using (26) when $S^* \not\subseteq S^t$, there exists at least one $i \in S^*$ such that $c_i^{t+1} = c_i^t + 1$. Using also that, for all $i \in S^*$, $(c_i^t)_{t \in \mathbb{N}}$ is non-decreasing we obtain

$$\text{for all } t \in \mathbb{N} \text{ such that } S^* \not\subseteq S^t, \quad \sum_{i \in S^*} c_i^{t+1} \ge (\sum_{i \in S^*} c_i^t) + 1 \qquad (28)$$

We define the first recovery iterate<fnref n="9" /> $t_s$ as the smallest iteration $t$ such that $S^* \subseteq S^t$. More precisely,

$$t_s = \min \{t, S^* \subseteq S^t\} \in \mathbb{N}. \qquad (29)$$

By convention, if $S^*$ is never recovered, $t_s = +\infty$. By induction on $t$, using (28), we obtain a lower bound on $\sum_{i \in S^*} c_i^t$:

$$\text{For all } t \le t_s, \quad \sum_{i \in S^*} c_i^t \ge t. \qquad (30)$$

Let us now upper bound $\sum_{i \in S^*} c_i^t$. We first remind the definition of $\varepsilon$ in (20). We define the sharp Recovery Condition

$$\varepsilon < \frac{1}{2 \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}} \qquad (RC^*)$$

and

$$T'_{max} = \frac{\sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0|}{\eta |x_i^*|} + k + 1}{1 - 2\varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}}. \qquad (31)$$

If (RC) holds, we define for any $i \in S^*$, the $i^{th}$ counting threshold by

$$C_i = \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0| + 2T'_{max}\varepsilon}{\eta|x_i^*|}. \qquad (32)$$

**Proposition C.5 (Upper bound).** If (RC') holds, for any $i \in S^*$ and any $t \le T'_{max}$, we have $c_i^t \le C_i + 1$.

*Proof.* Assume (RC) holds. We have $T'_{max} > 0$. Let $i \in S^*$, we distinguish two cases:

1.**st case:** If for all $t \le T'_{max}$, $c_i^t \le C_i$: Then, obviously, for any $t \le T'_{max}$, $c_i^t \le C_i + 1$.

2.**nd case:** If there exists $t \le T'_{max}$, such that $c_i^t > C_i$:

We define $t_i = \min \{t \in \mathbb{N} : c_i^t > C_i\}$. We have $t_i \le T'_{max}$. The proof follows two steps:

1. We will prove that for all $t \in [t_i, T'_{max}]$, $c_i^t = c_i^{t_i}$. \qquad (33)

2. We will prove that for all $t \le T'_{max}$, $c_i^t \le C_i + 1$. \qquad (34)

<footnote n="9">⁹Again, a similar quantity is defined in the proof of Theorem C.1 for the supports $S^t$ defined in Algorithm 6. We use the same notation although this time the quantity is defined for the sets $S^t$ defined in Algorithm 1. It should not be ambiguous since the notations are used in different proofs and sections.</footnote>
the same notations for simplicity. We set $c^0 = 0$, for any $t \in \mathbb{N}^*$ and $i \in [1, n]$, we also define the counting vector by

$$c_i^t = |\{t' \in [0, t-1] : i \in S^* \cap \overline{S^{t'}}\}|. \qquad (25)$$

We will use the recursive formula for $c^t$: For any $t \in \mathbb{N}$, $i \in [1, n]$

$$c_i^{t+1} = \begin{cases} c_i^t + 1 & \text{if } i \in S^* \cap \overline{S^t} \\ c_i^t & \text{if } i \in \overline{S^*} \cup S^t. \end{cases} \qquad (26)$$

For any $i \in [1, n]$, the sequence $(c_i^t)_{t \in \mathbb{N}}$ is non-decreasing.

Using (23), (24) and (25), we can establish by induction on $t$ that for any $t \in \mathbb{N}$,

$$\mathcal{X}^t = \mathcal{X}^0 + \eta c^t \odot x^* + \beta^t, \qquad (27)$$

where $\odot$ denotes the Hadamard product. This generalizes (17) to the noisy setting.

As can be seen from (27), the error accumulation $\beta^t$ is responsible for the exploration in the wrong directions. While $c^t \odot x^*$ encourages exploration in the direction of the missed components of $x^*$. Below, we provide important properties of $(c^t)_{t \in \mathbb{N}}$.

At each iteration of SEA, using (26) when $S^* \not\subseteq S^t$, there exists at least one $i \in S^*$ such that $c_i^{t+1} = c_i^t + 1$. Using also that, for all $i \in S^*$, $(c_i^t)_{t \in \mathbb{N}}$ is non-decreasing we obtain

$$\text{for all } t \in \mathbb{N} \text{ such that } S^* \not\subseteq S^t, \quad \sum_{i \in S^*} c_i^{t+1} \ge (\sum_{i \in S^*} c_i^t) + 1 \qquad (28)$$

We define the first recovery iterate<fnref n="9" /> $t_s$ as the smallest iteration $t$ such that $S^* \subseteq S^t$. More precisely,

$$t_s = \min \{t, S^* \subseteq S^t\} \in \mathbb{N}. \qquad (29)$$

By convention, if $S^*$ is never recovered, $t_s = +\infty$. By induction on $t$, using (28), we obtain a lower bound on $\sum_{i \in S^*} c_i^t$:

$$\text{For all } t \le t_s, \quad \sum_{i \in S^*} c_i^t \ge t. \qquad (30)$$

Let us now upper bound $\sum_{i \in S^*} c_i^t$. We first remind the definition of $\varepsilon$ in (20). We define the sharp Recovery Condition

$$\varepsilon < \frac{1}{2 \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}} \qquad (RC^*)$$

and

$$T'_{max} = \frac{\sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0|}{\eta |x_i^*|} + k + 1}{1 - 2\varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}}. \qquad (31)$$

If (RC) holds, we define for any $i \in S^*$, the $i^{th}$ counting threshold by

$$C_i = \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0| + 2T'_{max}\varepsilon}{\eta|x_i^*|}. \qquad (32)$$

**Proposition C.5 (Upper bound).** If (RC') holds, for any $i \in S^*$ and any $t \le T'_{max}$, we have $c_i^t \le C_i + 1$.

*Proof.* Assume (RC) holds. We have $T'_{max} > 0$. Let $i \in S^*$, we distinguish two cases:

1.**st case:** If for all $t \le T'_{max}$, $c_i^t \le C_i$: Then, obviously, for any $t \le T'_{max}$, $c_i^t \le C_i + 1$.

2.**nd case:** If there exists $t \le T'_{max}$, such that $c_i^t > C_i$:

We define $t_i = \min \{t \in \mathbb{N} : c_i^t > C_i\}$. We have $t_i \le T'_{max}$. The proof follows two steps:

1. We will prove that for all $t \in [t_i, T'_{max}]$, $c_i^t = c_i^{t_i}$. \qquad (33)

2. We will prove that for all $t \le T'_{max}$, $c_i^t \le C_i + 1$. \qquad (34)

<footnote n="9">⁹Again, a similar quantity is defined in the proof of Theorem C.1 for the supports $S^t$ defined in Algorithm 6. We use the same notation although this time the quantity is defined for the sets $S^t$ defined in Algorithm 1. It should not be ambiguous since the notations are used in different proofs and sections.</footnote>
1. Let $t \in [t_i, T'_{max}]$, we have, using (27), the triangle inequality and the fact that $c_i^t \ge c_i^{t_i} > C_i$

$$ 
\begin{aligned}
|\mathcal{X}_i^t| &= |\mathcal{X}_i^0 + \eta c_i^t x_i^* + \beta_i^t| \\
&\geq \eta c_i^t |x_i^*| - |\mathcal{X}_i^0| - |\beta_i^t| \\
&> \eta C_i |x_i^*| - |\mathcal{X}_i^0| - |\beta_i^t|.
\end{aligned}
$$

Using the definition of $C_i$, in (32), we obtain

$$ 
\begin{aligned}
|\mathcal{X}_i^t| &> \eta \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0| + 2T'_{max} \varepsilon}{\eta |x_i^*|} |x_i^*| - |\mathcal{X}_i^0| - |\beta_i^t| \\
&= \max_{j \notin S^*} |\mathcal{X}_j^0| + 2T'_{max} \varepsilon - |\beta_i^t|.
\end{aligned}
$$

Since for any $j \in [1, n]$, $|\beta_j^t| \le \sum_{t'=0}^{t-1} |b_j^{t'}| \le t\varepsilon \le T'_{max} \varepsilon$, we have

$$ 
\begin{aligned}
|\mathcal{X}_i^t| &> \max_{j \notin S^*} |\mathcal{X}_j^0| + \max_{j \notin S^*} |\beta_j^t| + |\beta_i^t| - |\beta_i^t| \\
&\geq \max_{j \notin S^*} |\mathcal{X}_j^0 + \beta_j^t| \\
&= \max_{j \notin S^*} |\mathcal{X}_j^t|,
\end{aligned}
\quad (35) $$

where the last equality holds because of (27) and for all $j \notin S^*$, all $t \in \mathbb{N}$, $c_j^t = 0$.

Equation (35) implies that $|\mathcal{X}_i^t|$ is larger than $|\{j \notin S^*\}|$ elements of $\{|\mathcal{X}_j^t| \mid j \in [1, n]\}$ and, since $|S^*| \le k$, we have $|\{j \notin S^*\}| = n - |S^*| \ge n - k$. Finally, $|\mathcal{X}_i^t|$ is larger than $n - k$ elements of $\{|\mathcal{X}_j^t| \mid j \in [1, n]\}$ and $i \in \text{largest}_k(\mathcal{X}^t) = S^t$.

As a conclusion, for all $t \in [t_i, T'_{max}]$, $i \in S^t$. Using (26), this leads to $c_i^{t+1} = c_i^t$. Therefore, for all $t \in [t_i, T'_{max} + 1]$, $c_i^t = c_i^{t_i}$. This concludes the proof of the first step.

2. Since $t_i = \min\{t \in \mathbb{N} : c_i^t > C_i\}$ and since $c_i^0 = 0$, $t_i \ge 1$. Since by definition of $t_i$, $c_i^{t_i-1} \le C_i$ and $c_i^{t_i} \ne c_i^{t_i-1}$; we find that $c_i^{t_i} = c_i^{t_i-1} + 1 \le C_i + 1$.

Using (33), for all $t \in [t_i, T'_{max}]$, $c_i^t = c_i^{t_i} \le C_i + 1$. Finally, since $(c_i^t)_{t \in \mathbb{N}^*}$ is non-decreasing, it follows that for any $t \le t_i - 1$, $c_i^t \le c_i^{t_i-1} \le C_i$. This concludes the proof of (34).

□

### C.2.3 Proof of Theorem C.4

To prove Theorem C.4, we first prove a sharper, but difficult-to-interpret theorem.

**Theorem C.6 (Recovery - General case).** If (RC') holds, then for all initializations $\mathcal{X}^0$ and all $\eta > 0$, there exists $t_s \le T'_{max}$ such that $S^* \subseteq S^{t_s}$, where $S^t$ is defined in Algorithm 1 line 6 and

$$ T'_{max} = \frac{\sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0|}{\eta |x_i^*|} + k + 1}{1 - 2\varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}}. $$

*Proof.* We assume (RC') holds and prove Theorem C.6 using the results of Appendix C.2.2.

In order to do this, we first show that $T'_{max} = \sum_{i \in S^*} C_i + k + 1$, then we demonstrate that $t_s \le T'_{max}$. Since (RC') holds, using the definition of $T'_{max}$, we calculate

$$ 
\begin{aligned}
T'_{max} &= \frac{1}{1 - 2\varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}} \left( \sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0|}{\eta |x_i^*|} + k + 1 \right) \\
\left(1 - 2\varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}\right) T'_{max} &= \sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0|}{\eta |x_i^*|} + k + 1 \\
T'_{max} &= \sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0|}{\eta |x_i^*|} + k + 1 + 2T'_{max} \varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|} \\
T'_{max} &= \sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0| + 2T'_{max} \varepsilon}{\eta |x_i^*|} + k + 1.
\end{aligned}
$$
1. Let $t \in [t_i, T'_{max}]$, we have, using (27), the triangle inequality and the fact that $c_i^t \ge c_i^{t_i} > C_i$

$$ \begin{aligned} |\mathcal{X}_i^t| &= |\mathcal{X}_i^0 + \eta c_i^t x_i^* + \beta_i^t| \\ &\ge \eta c_i^t |x_i^*| - |\mathcal{X}_i^0| - |\beta_i^t| \\ &> \eta C_i |x_i^*| - |\mathcal{X}_i^0| - |\beta_i^t|. \end{aligned} $$

Using the definition of $C_i$, in (32), we obtain

$$ \begin{aligned} |\mathcal{X}_i^t| &> \eta \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0| + 2T'_{max} \varepsilon}{\eta |x_i^*|} |x_i^*| - |\mathcal{X}_i^0| - |\beta_i^t| \\ &= \max_{j \notin S^*} |\mathcal{X}_j^0| + 2T'_{max} \varepsilon - |\beta_i^t|. \end{aligned} $$

Since for any $j \in [1, n]$, $|\beta_j^t| \le \sum_{t'=0}^{t-1} |b_j^{t'}| \le t\varepsilon \le T'_{max} \varepsilon$, we have

$$ \begin{aligned} |\mathcal{X}_i^t| &> \max_{j \notin S^*} |\mathcal{X}_j^0| + \max_{j \notin S^*} |\beta_j^t| + |\beta_i^t| - |\beta_i^t| \\ &\ge \max_{j \notin S^*} |\mathcal{X}_j^0 + \beta_j^t| \\ &= \max_{j \notin S^*} |\mathcal{X}_j^t|, \end{aligned} \quad (35) $$

where the last equality holds because of (27) and for all $j \notin S^*$, all $t \in \mathbb{N}$, $c_j^t = 0$.

Equation (35) implies that $|\mathcal{X}_i^t|$ is larger than $|\{j \notin S^*\}|$ elements of $\{|\mathcal{X}_j^t| \mid j \in [1, n]\}$ and, since $|S^*| \le k$, we have $|\{j \notin S^*\}| = n - |S^*| \ge n - k$. Finally, $|\mathcal{X}_i^t|$ is larger than $n - k$ elements of $\{|\mathcal{X}_j^t| \mid j \in [1, n]\}$ and $i \in \text{largest}_k(\mathcal{X}^t) = S^t$.

As a conclusion, for all $t \in [t_i, T'_{max}]$, $i \in S^t$. Using (26), this leads to $c_i^{t+1} = c_i^t$. Therefore, for all $t \in [t_i, T'_{max} + 1]$, $c_i^t = c_i^{t_i}$. This concludes the proof of the first step.

2. Since $t_i = \min\{t \in \mathbb{N} : c_i^t > C_i\}$ and since $c_i^0 = 0$, $t_i \ge 1$. Since by definition of $t_i$, $c_i^{t_i-1} \le C_i$ and $c_i^{t_i} \ne c_i^{t_i-1}$; we find that $c_i^{t_i} = c_i^{t_i-1} + 1 \le C_i + 1$.

Using (33), for all $t \in [t_i, T'_{max}]$, $c_i^t = c_i^{t_i} \le C_i + 1$. Finally, since $(c_i^t)_{t \in \mathbb{N}^*}$ is non-decreasing, it follows that for any $t \le t_i - 1$, $c_i^t \le c_i^{t_i-1} \le C_i$. This concludes the proof of (34).

□

### C.2.3 Proof of Theorem C.4

To prove Theorem C.4, we first prove a sharper, but difficult-to-interpret theorem.

**Theorem C.6 (Recovery - General case).** If (RC') holds, then for all initializations $\mathcal{X}^0$ and all $\eta > 0$, there exists $t_s \le T'_{max}$ such that $S^* \subseteq S^{t_s}$, where $S^t$ is defined in Algorithm 1 line 6 and

$$ T'_{max} = \frac{\sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0|}{\eta |x_i^*|} + k + 1}{1 - 2\varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}}. $$

*Proof.* We assume (RC') holds and prove Theorem C.6 using the results of Appendix C.2.2.

In order to do this, we first show that $T'_{max} = \sum_{i \in S^*} C_i + k + 1$, then we demonstrate that $t_s \le T'_{max}$. Since (RC') holds, using the definition of $T'_{max}$, we calculate

$$ \begin{aligned} T'_{max} &= \frac{1}{1 - 2\varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}} \left( \sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0|}{\eta |x_i^*|} + k + 1 \right) \\ \left(1 - 2\varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}\right) T'_{max} &= \sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0|}{\eta |x_i^*|} + k + 1 \\ T'_{max} &= \sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0|}{\eta |x_i^*|} + k + 1 + 2T'_{max} \varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|} \\ T'_{max} &= \sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0| + 2T'_{max} \varepsilon}{\eta |x_i^*|} + k + 1. \end{aligned} $$
Using (32), we obtain $T'_{max} = \sum_{i \in S^*} C_i + k + 1$.

We finally prove Theorem C.4 by contradiction. Assume by contradiction that $t_s > T'_{max}$, where $t_s$ is defined in (29). Using (30) with $t = \lfloor T'_{max} \rfloor < t_s$, we have

$$ \sum_{i \in S^*} c_i^{\lfloor T'_{max} \rfloor} \geq \lfloor T'_{max} \rfloor = \lfloor \sum_{i \in S^*} C_i + k + 1 \rfloor > \sum_{i \in S^*} C_i + k. \quad (36) $$

However, using $|S^*| \leq k$ and Proposition C.5 for $t = \lfloor T'_{max} \rfloor$, we find

$$ \sum_{i \in S^*} C_i + k \geq \sum_{i \in S^*} (C_i + 1) \geq \sum_{i \in S^*} c_i^{\lfloor T'_{max} \rfloor} $$

This contradicts (36) and we can conclude that $t_s \leq T'_{max}$. This proves Theorem C.6. □

**Proof of Theorem C.4:**

If (RC) holds, that is $\varepsilon < \frac{\eta}{2k} \min_{i \in S^*} |x_i^*|$, since $\sum_{i \in S^*} \frac{1}{|x_i^*|} \leq \frac{k}{\min_{i \in S^*} |x_i^*|}$, we have

$$ \varepsilon < \frac{1}{2 \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}}. $$

Therefore, (RC) holds, and we can apply Theorem C.6. It ensures that for all initializations $\mathcal{X}^0$ and all $\eta > 0$, there exists $t_s \leq T'_{max}$ such that $S^* \subseteq S^{t_s}$.

To prove Theorem C.4, it suffices to prove that $T'_{max} \leq T_{max}$. Using $\sum_{i \in S^*} \frac{1}{|x_i^*|} \leq \frac{k}{\min_{i \in S^*} |x_i^*|}$, we obtain

$$ 1 - 2\varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|} \geq 1 - 2\varepsilon \frac{k}{\eta \min_{i \in S^*} |x_i^*|} $$

and using $\min_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0| \leq 2 \|\mathcal{X}^0\|_{\infty}$ we have

$$ T'_{max} = \frac{\sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0|}{\eta |x_i^*|} + k + 1}{1 - 2\varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}} \leq \frac{\frac{2k \|\mathcal{X}^0\|_{\infty}}{\eta \min_{i \in S^*} |x_i^*|} + k + 1}{1 - 2\varepsilon \frac{k}{\eta \min_{i \in S^*} |x_i^*|}} = T_{max}. $$

### C.3 Warm-Up: SEA Recovers the Correct Support when the Columns of A are Orthonormal

The following corollary particularizes Theorem C.4 to the noiseless and orthogonal case. In practice, a complicated algorithm like SEA is of course useless in such a case, and the state-of-the-art algorithms mentioned in the introduction have similar recovery properties. We give this corollary mostly to illustrate the diversity of links between the properties of the triplet $(A, x^*, e)$ and $\varepsilon$ and the behavior of SEA, where we remind the definitions of $b^t$ and $\varepsilon$ in (19) and (20). The following Corollary C.7 is not only a sanity check for the convergence of SEA under simplistic assumptions, but it also provides a helpful case to understand the proof of Theorem 4.1. Indeed, it gives a case where the oracle updates introduced in the proof coincides the surrogate, STE update, as mentioned in the last paragraph of Appendix C.1.

**Corollary C.7 (Recovery - Orthogonal case).** If A is a tall (or square) matrix with orthonormal columns ($A^T A = I_n$) and $\|e\|_2 = 0$, then $\varepsilon = 0$. As a consequence, for all $x^*$, for initialization $\mathcal{X}^0 = 0$ and all $\eta > 0$, if SEA performs more than $k + 1$ iterations, we have $S^* \subseteq S^{t_{BEST}}$ and $x^{t_{BEST}} = x^*$.

To prove Corollary C.7, we first show in Lemma C.8 that the gradient noise $b^t$ is null for all $t \in \mathbb{N}$. Then, we apply Theorem C.4 and prove that $S^* \subseteq S^{t_{BEST}}$ and $x^{t_{BEST}} = x^*$.

**Lemma C.8.** If A is a tall (or square) matrix with orthonormal columns ($A^T A = I_n$) and $\|e\|_2 = 0$, then for any $t \in \mathbb{N}$ and any $\eta > 0$,

$$ \eta A^T (Ax^t - y) = u^t, $$

i.e. $b^t = 0$.

*Proof.* Let $t \in \mathbb{N}$. Notice first that since $\|e\|_2 = 0$ and A is a tall (or square) matrix with orthonormal columns

$$ A^T (Ax^t - y) = A^T A (x^t - x^*) = x^t - x^*. \quad (37) $$
Using (32), we obtain $T'_{max} = \sum_{i \in S^*} C_i + k + 1$.

We finally prove Theorem C.4 by contradiction. Assume by contradiction that $t_s > T'_{max}$, where $t_s$ is defined in (29). Using (30) with $t = \lfloor T'_{max} \rfloor < t_s$, we have

$$ \sum_{i \in S^*} c_i^{\lfloor T'_{max} \rfloor} \geq \lfloor T'_{max} \rfloor = \lfloor \sum_{i \in S^*} C_i + k + 1 \rfloor > \sum_{i \in S^*} C_i + k. \quad (36) $$

However, using $|S^*| \leq k$ and Proposition C.5 for $t = \lfloor T'_{max} \rfloor$, we find

$$ \sum_{i \in S^*} C_i + k \geq \sum_{i \in S^*} (C_i + 1) \geq \sum_{i \in S^*} c_i^{\lfloor T'_{max} \rfloor} $$

This contradicts (36) and we can conclude that $t_s \leq T'_{max}$. This proves Theorem C.6. □

**Proof of Theorem C.4:**

If (RC) holds, that is $\varepsilon < \frac{\eta}{2k} \min_{i \in S^*} |x_i^*|$, since $\sum_{i \in S^*} \frac{1}{|x_i^*|} \leq \frac{k}{\min_{i \in S^*} |x_i^*|}$, we have

$$ \varepsilon < \frac{1}{2 \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}}. $$

Therefore, (RC) holds, and we can apply Theorem C.6. It ensures that for all initializations $\mathcal{X}^0$ and all $\eta > 0$, there exists $t_s \leq T'_{max}$ such that $S^* \subseteq S^{t_s}$.

To prove Theorem C.4, it suffices to prove that $T'_{max} \leq T_{max}$. Using $\sum_{i \in S^*} \frac{1}{|x_i^*|} \leq \frac{k}{\min_{i \in S^*} |x_i^*|}$, we obtain

$$ 1 - 2\varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|} \geq 1 - 2\varepsilon \frac{k}{\eta \min_{i \in S^*} |x_i^*|} $$

and using $\min_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0| \leq 2 \|\mathcal{X}^0\|_{\infty}$ we have

$$ T'_{max} = \frac{\sum_{i \in S^*} \frac{\max_{j \notin S^*} |\mathcal{X}_j^0| + |\mathcal{X}_i^0|}{\eta |x_i^*|} + k + 1}{1 - 2\varepsilon \sum_{i \in S^*} \frac{1}{\eta |x_i^*|}} \leq \frac{\frac{2k \|\mathcal{X}^0\|_{\infty}}{\eta \min_{i \in S^*} |x_i^*|} + k + 1}{1 - 2\varepsilon \frac{k}{\eta \min_{i \in S^*} |x_i^*|}} = T_{max}. $$

## C.3 Warm-Up: SEA Recovers the Correct Support when the Columns of A are Orthonormal

The following corollary particularizes Theorem C.4 to the noiseless and orthogonal case. In practice, a complicated algorithm like SEA is of course useless in such a case, and the state-of-the-art algorithms mentioned in the introduction have similar recovery properties. We give this corollary mostly to illustrate the diversity of links between the properties of the triplet $(A, x^*, e)$ and $\varepsilon$ and the behavior of SEA, where we remind the definitions of $b^t$ and $\varepsilon$ in (19) and (20). The following Corollary C.7 is not only a sanity check for the convergence of SEA under simplistic assumptions, but it also provides a helpful case to understand the proof of Theorem 4.1. Indeed, it gives a case where the oracle updates introduced in the proof coincides the surrogate, STE update, as mentioned in the last paragraph of Appendix C.1.

**Corollary C.7 (Recovery - Orthogonal case).** If A is a tall (or square) matrix with orthonormal columns ($A^T A = I_n$) and $\|e\|_2 = 0$, then $\varepsilon = 0$. As a consequence, for all $x^*$, for initialization $\mathcal{X}^0 = 0$ and all $\eta > 0$, if SEA performs more than $k + 1$ iterations, we have $S^* \subseteq S^{t_{BEST}}$ and $x^{t_{BEST}} = x^*$.

To prove Corollary C.7, we first show in Lemma C.8 that the gradient noise $b^t$ is null for all $t \in \mathbb{N}$. Then, we apply Theorem C.4 and prove that $S^* \subseteq S^{t_{BEST}}$ and $x^{t_{BEST}} = x^*$.

**Lemma C.8.** If A is a tall (or square) matrix with orthonormal columns ($A^T A = I_n$) and $\|e\|_2 = 0$, then for any $t \in \mathbb{N}$ and any $\eta > 0$,

$$ \eta A^T (Ax^t - y) = u^t, $$

i.e. $b^t = 0$.

*Proof.* Let $t \in \mathbb{N}$. Notice first that since $\|e\|_2 = 0$ and A is a tall (or square) matrix with orthonormal columns

$$ A^T (Ax^t - y) = A^T A (x^t - x^*) = x^t - x^*. \quad (37) $$
To prove the Lemma, we distinguish three cases: $i \in S^t$, $i \in \overline{S^*} \cap \overline{S^t}$ and $i \in S^* \cap \overline{S^t}$.

**1st case:** If $i \in S^t$, $\eta(A^T(Ax^t - y))_i = 0 = u_i^t$. The first equality is a consequence of the definition of $x^t$ in Algorithm 1, line 7. The second is due to the definition of $u^t$, in (15).

**2nd case:** If $i \in \overline{S^*} \cap \overline{S^t}$, taking the ith entry of (37) and using the support constraints of $x^t$ and $x^*$, we find

$$\eta(A^T(Ax^t - y))_i = 0 = u_i^t,$$

where the second equality is due to the definition of $u^t$, in (15).

**3rd case:** If $i \in S^* \cap \overline{S^t}$, the ith entry of (37) becomes

$$\eta(A^T(Ax^t - y))_i = -\eta x_i^* = u_i^t,$$

where again the second equality is due to the definition of $u^t$, in (15).

□

*Proof.* We now resume the proof of Corollary C.7 and assume that A is a tall (or square) matrix with orthonormal columns ($A^T A = I_n$), $\|e\|_2 = 0$ and $\mathcal{X}^0 = 0$. We remind the definition of $T_{max}$ in (21).

Using Lemma C.8, (20) and (19), we find that $\varepsilon = 0$. Therefore (RC) holds for all $x^*$ and Theorem C.4 implies that there exists $t_s \le T_{max}$ such that $S^* \subseteq S^{t_s}$. Since $\mathcal{X}^0 = 0$ and $\varepsilon = 0$, we find $T_{max} = k + 1$.

Since $\|e\|_2 = 0$, we know from Theorem C.4 and the definitions of $t_{BEST}$ and $x^t$ in Algorithm 1 that

$$\|Ax^{t_{BEST}} - y\|_2 \le \|Ax^{t_s} - y\|_2 \le \|Ax^* - y\|_2 = 0.$$

Using that A is a tall (or square) matrix with orthonormal columns ($A^T A = I_n$), and $\|e\|_2 = 0$, this leads to

$$\begin{aligned} 0 &= Ax^{t_{BEST}} - y \\ &= A^T A(x^{t_{BEST}} - x^*) \\ &= x^{t_{BEST}} - x^* \end{aligned}$$

Therefore, $S^* = \text{supp}(x^*) = \text{supp}(x^{t_{BEST}}) \subseteq S^{t_{BEST}}$.

This concludes the proof of Corollary C.7.

□

## C.4 Proof of Theorem 4.1

We remind in Appendix C.4.1 known properties of RIP matrices. We bound in Appendix C.4.2 the error made when approximating $x^*$ on a specific support $S$. This permits us to bound $b^t$ and apply Theorem C.4 to prove Theorem 4.1 in Appendix C.4.3. We finally apply Theorem 4.1 in Appendix C.5 to prove Corollary 4.2. Before that and to illustrate and quantify that the STE-update $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - \eta A^T(Ax^t - y)$ is a noisy version of the Oracle update $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - u^t$, we provide in the following theorem an upper bound on the discrepancy between the two updates. This bound is pivotal in the proof of Theorem 4.1. The statement of Theorem 4.1 is, up to the additional upper-bound (38), the same as the statement of the following theorem, which we prove in this section.

**Theorem C.9 (Recovery - RIP case).** Assume A satisfies the (2k + 1)-RIP and for all $i \in [1, n]$, $\|A_i\|_2 = 1$. Then, for all $t \in \mathbb{N}$,

$$\| \frac{u^t}{\eta} - A^T(Ax^t - y) \|_{\infty} \le \alpha_k^{\text{RIP}} \|x^*\|_2 + \gamma_k^{\text{RIP}} \|e\|_2. \quad (38)$$

If moreover $x^*$ satisfies (RC$_{RIP}$), then for all initializations $\mathcal{X}^0$ and all $\eta > 0$, there exists $t_s \le T_{RIP}$ such that $S^* \subseteq S^{t_s}$, where

$$T_{RIP} = \frac{2k \frac{\|\mathcal{X}^0\|_{\infty}}{\eta} + (k+1) \min_{i \in S^*} |x_i^*|}{\min_{i \in S^*} |x_i^*| - 2k (\alpha_k^{RIP} \|x^*\|_2 + \gamma_k^{RIP} \|e\|_2)} \quad (39)$$

If moreover, $x^*$ is such that

$$\min_{i \in S^*} |x_i^*| > \frac{2}{\sqrt{1 - \delta_{2k}}} \|e\|_2 \quad (40)$$

and SEA performs more than $T_{RIP}$ iterations, then $S^* \subseteq S^{t_{BEST}}$ and $\|x^{t_{BEST}} - x^*\|_2 \le \frac{2}{\sqrt{1 - \delta_k}} \|e\|_2$.
To prove the Lemma, we distinguish three cases: $i \in S^t$, $i \in \overline{S^*} \cap \overline{S^t}$ and $i \in S^* \cap \overline{S^t}$.

**1st case:** If $i \in S^t$, $\eta(A^T(Ax^t - y))_i = 0 = u_i^t$. The first equality is a consequence of the definition of $x^t$ in Algorithm 1, line 7. The second is due to the definition of $u^t$, in (15).

**2nd case:** If $i \in \overline{S^*} \cap \overline{S^t}$, taking the ith entry of (37) and using the support constraints of $x^t$ and $x^*$, we find

$$\eta(A^T(Ax^t - y))_i = 0 = u_i^t,$$

where the second equality is due to the definition of $u^t$, in (15).

**3rd case:** If $i \in S^* \cap \overline{S^t}$, the ith entry of (37) becomes

$$\eta(A^T(Ax^t - y))_i = -\eta x_i^* = u_i^t,$$

where again the second equality is due to the definition of $u^t$, in (15).

□

*Proof.* We now resume the proof of Corollary C.7 and assume that $A$ is a tall (or square) matrix with orthonormal columns $(A^T A = I_n)$, $\|e\|_2 = 0$ and $\mathcal{X}^0 = 0$. We remind the definition of $T_{max}$ in (21).

Using Lemma C.8, (20) and (19), we find that $\varepsilon = 0$. Therefore (RC) holds for all $x^*$ and Theorem C.4 implies that there exists $t_s \le T_{max}$ such that $S^* \subseteq S^{t_s}$. Since $\mathcal{X}^0 = 0$ and $\varepsilon = 0$, we find $T_{max} = k + 1$.

Since $\|e\|_2 = 0$, we know from Theorem C.4 and the definitions of $t_{BEST}$ and $x^t$ in Algorithm 1 that

$$\|Ax^{t_{BEST}} - y\|_2 \le \|Ax^{t_s} - y\|_2 \le \|Ax^* - y\|_2 = 0.$$

Using that $A$ is a tall (or square) matrix with orthonormal columns $(A^T A = I_n)$, and $\|e\|_2 = 0$, this leads to

$$\begin{aligned} 0 &= Ax^{t_{BEST}} - y \\ &= A^T A(x^{t_{BEST}} - x^*) \\ &= x^{t_{BEST}} - x^* \end{aligned}$$

Therefore, $S^* = \text{supp}(x^*) = \text{supp}(x^{t_{BEST}}) \subseteq S^{t_{BEST}}$.

This concludes the proof of Corollary C.7.

□

## C.4 Proof of Theorem 4.1

We remind in Appendix C.4.1 known properties of RIP matrices. We bound in Appendix C.4.2 the error made when approximating $x^*$ on a specific support $S$. This permits us to bound $b^t$ and apply Theorem C.4 to prove Theorem 4.1 in Appendix C.4.3. We finally apply Theorem 4.1 in Appendix C.5 to prove Corollary 4.2. Before that and to illustrate and quantify that the STE-update $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - \eta A^T(Ax^t - y)$ is a noisy version of the Oracle update $\mathcal{X}^{t+1} \leftarrow \mathcal{X}^t - u^t$, we provide in the following theorem an upper bound on the discrepancy between the two updates. This bound is pivotal in the proof of Theorem 4.1. The statement of Theorem 4.1 is, up to the additional upper-bound (38), the same as the statement of the following theorem, which we prove in this section.

**Theorem C.9 (Recovery - RIP case).** Assume $A$ satisfies the (2k + 1)-RIP and for all $i \in [1, n]$, $\|A_i\|_2 = 1$. Then, for all $t \in \mathbb{N}$,

$$\| \frac{u^t}{\eta} - A^T(Ax^t - y) \|_{\infty} \le \alpha_k^{\text{RIP}} \|x^*\|_2 + \gamma_k^{\text{RIP}} \|e\|_2. \quad (38)$$

If moreover $x^*$ satisfies $(RC_{RIP})$, then for all initializations $\mathcal{X}^0$ and all $\eta > 0$, there exists $t_s \le T_{RIP}$ such that $S^* \subseteq S^{t_s}$, where

$$T_{RIP} = \frac{2k \frac{\|\mathcal{X}^0\|_{\infty}}{\eta} + (k+1) \min_{i \in S^*} |x_i^*|}{\min_{i \in S^*} |x_i^*| - 2k (\alpha_k^{RIP} \|x^*\|_2 + \gamma_k^{RIP} \|e\|_2)} \quad (39)$$

If moreover, $x^*$ is such that

$$\min_{i \in S^*} |x_i^*| > \frac{2}{\sqrt{1 - \delta_{2k}}} \|e\|_2 \quad (40)$$

and SEA performs more than $T_{RIP}$ iterations, then $S^* \subseteq S^{t_{BEST}}$ and $\|x^{t_{BEST}} - x^*\|_2 \le \frac{2}{\sqrt{1 - \delta_k}} \|e\|_2$.