On Iterative Hard Thresholding Methods for High-dimensional M-Estimation

Introduction

Modern statistical estimation is routinely faced with real world problems where the number of parameters $p$ handily outnumbers the number of observations $n$. In general, consistent estimation of parameters is not possible in such a situation. Consequently, a rich line of work has focused on models that satisfy special structural assumptions such as sparsity or low-rank structure. Under these assumptions, several works (for example, see ) have established that consistent estimation is information theoretically possible in the “$n\ll p$” regime as well.

The question of efficient estimation, however, is faced with feasibility issues since consistent estimation routines often end-up solving NP-hard problems. Examples include sparse regression which requires loss minimization with sparsity constraints and low-rank regression which requires dealing with rank constraints which are not efficiently solvable in general .

Interestingly, recent works have demonstrated that these hardness results can be avoided by assuming certain natural conditions over the loss function being minimized such as restricted strong convexity (RSC) and restricted strong smoothness (RSS). The estimation routines proposed in these works typically make use of convex relaxations or greedy methods which do not suffer from infeasibility issues.

Despite this, certain limitations have precluded widespread use of these techniques. Convex relaxation-based methods typically suffer from slow rates as they solve non-smooth optimization problems apart from being hard to analyze in terms of global guarantees. Greedy methods, on the other hand, are slow in situations with non-negligible sparsity or relatively high rank, owing to their incremental approach of adding/removing individual support elements.

Instead, the methods of choice for practical applications are actually projected gradient (PGD) methods, also referred to as iterative hard thresholding (IHT) methods. These methods directly project the gradient descent update onto the underlying (non-convex) feasible set. This projection can be performed efficiently for several interesting structures such as sparsity and low rank. However, traditional PGD analyses for convex problems viz. do not apply to these techniques due to the non-convex structure of the problem.

An exception to this is the recent work that demonstrates that PGD with non-convex regularization can offer consistent estimates for certain high-dimensional problems. However, the work in is only able to analyze penalties such as SCAD, MCP and capped $L_{1}$. Moreover, their framework cannot handle commonly used penalties such as $L_{0}$ or low-rank constraints.

All existing results analyzing these methods require the condition number of the loss function, restricted to sparse vectors, to be smaller than a universal constant. The best known such constant is due to the work of that requires a bound on the RIP constant $\delta_{2k}\leq 0.5$ (or equivalently a bound $\tfrac{1+\delta_{2k}}{1-\delta_{2k}}\leq 3$ on the condition number). In contrast, real-life high dimensional statistical settings, wherein pairs of variables can be arbitrarily correlated, routinely require estimation methods to perform under arbitrarily large condition numbers. In particular if two variates have a covariance matrix like ${\left[\begin{matrix}1&1-\epsilon\\ 1-\epsilon&1\end{matrix}\right]}$, then the restricted condition number (on a support set of size just $2$) of the sample matrix cannot be brought down below $1/\epsilon$ even with infinitely many samples. In particular when $\epsilon<1/6$, none of the existing results for hard thresholding methods offer any guarantees. Moreover, most of these analyses consider only the least squares objective. Although recent attempts have been made to extend this to general differentiable objectives , the results continue to require that the restricted condition number be less than a universal constant and remain unsatisfactory in a statistical setting.

2 Overview of Results

Our main contribution in this work is an analysis of PGD/IHT-style methods in statistical settings. Our bounds are tight, achieve known minmax lower bounds , and hold for arbitrary differentiable, possibly even non-convex functions. Our results hold even when the underlying condition number is arbitrarily large and only require the function to satisfy RSC/RSS conditions. In particular, this reveals that these iterative methods are indeed applicable to statistical settings, a result that escaped all previous works.

Our first result shows that the PGD/IHT methods achieve global convergence if used with a relaxed projection step. More formally, if the optimal parameter is $s^{*}$-sparse and the problem satisfies RSC and RSS constraints $\alpha$ and $L$ respectively (see Section 2), then PGD methods offer global convergence so long as they employ projection to an $s$-sparse set where $s\geq 4(L/\alpha)^{2}s^{*}$. This gives convergence rates that are identical to those of convex relaxation and greedy methods for the Gaussian sparse linear model. We then move to a family of efficient “fully corrective” methods and show as before, that for arbitrary functions satisfying the RSC/RSS properties, these methods offer global convergence.

Next, we show that these results allow PGD-style methods to offer global convergence in a variety of statistical estimation problems such as sparse linear regression and low rank matrix regression. Our results effortlessly extend to the noisy setting as a corollary and give bounds similar to those of that relies on solving an $L_{1}$ regularized problem.

Our proofs are able to exploit that even though hard-thresholding is not the prox-operator for any convex prox function, it still provides strong contraction when projection is performed onto sets of sparsity $s\gg s^{*}$. This crucial observation allows us to provide the first unified analysis for hard thresholding based gradient descent algorithms. Our empirical results confirm our predictions with respect to the recovery properties of IHT-style algorithms on badly-conditioned sparse recovery problems, as well as demonstrate that these methods can be orders of magnitudes faster than their $L_{1}$ and greedy counterparts.

3 Organization of the Paper

Section 2 sets the notation and the problem statement. Section 3 introduces the PGD/IHT algorithm that we study and proves that the method guarantees recovery assuming the RSC/RSS property. We also generalize our guarantees to the problem of low-rank matrix regression. Section 4 then provides crisp sample complexity bounds and statistical guarantees for the PGD/IHT estimators. Section 5 extends our analysis to a broad family of “fully-corrective” hard thresholding methods compressive sensing algorithms that include the so-called two-stage hard thresholding and partial hard thresholding algorithms and provide similar results for them as well. We present some empirical results in Section 6 and conclude in Section 7.

Problem Setup and Notations

For this problem the RSC and RSS properties for $f$ are defined similarly as in Definition 1, 2 except that the $L_{0}$ norm is replaced by the rank function.

Iterative Hard-thresholding Method

In this section we study the popular projected gradient descent (a.k.a iterative hard thresholding) method for the case of the feasible set being the set of sparse vectors (see Algorithm 1 for pseudocode). The projection operator $P_{s}(\bm{z})$, can be implemented efficiently in this case by projecting $\bm{z}$ onto the set of $s$-sparse vectors by selecting the $s$ largest elements (in magnitude) of $\bm{z}$. The standard projection property implies that $\|P_{s}(\bm{z})-\bm{z}\|_{2}^{2}\leq\|\bm{\theta}^{\prime}-\bm{z}\|_{2}^{2}$ for all $\|\bm{\theta}^{\prime}\|_{0}\leq s$. However, it turns out that we can prove a significantly stronger property of hard thresholding for the case when $\|\bm{\theta}^{\prime}\|_{0}\leq s^{*}$ and $s^{*}\ll s$. This property is key to analysing IHT and is formalized below.

Our analysis combines the above observation with the RSC/RSS properties of $f$ to provide geometric convergence rates for the IHT procedure below.

Let $f$ have RSC and RSS parameters given by $L_{2s+s^{*}}(f)=L$ and $\alpha_{2s+s^{*}}(f)=\alpha$ respectively. Let Algorithm 1 be invoked with $f$, $s\geq 32\left(\frac{L}{\alpha}\right)^{2}s^{*}$ and $\eta=\frac{2}{3L}$. Also let $\bm{\theta}^{*}=\arg\min_{\bm{\theta},\|\bm{\theta}\|_{0}\leq s^{*}}f(\bm{\theta})$. Then, the $\tau$-th iterate of Algorithm 1, for $\tau=O(\frac{L}{\alpha}\cdot\log(\frac{f(\bm{\theta}^{0})}{\epsilon}))$ satisfies:

(Sketch) Let $S^{t}=supp(\bm{\theta}^{t})$, $S^{*}=supp(\bm{\theta}^{*})$, $S^{t+1}=supp(\bm{\theta}^{t+1})$ and $I^{t}=S^{*}\cup S^{t}\cup S^{t+1}$. Using the RSS property and the fact that $supp(\bm{\theta}^{t})\subseteq I^{t}$ and $supp(\bm{\theta}^{t+1})\subseteq I^{t}$, we have:

where $\zeta_{1}$ follows from an application of Lemma 1 with $I=I^{t}$ and the Pythagoras theorem. The above equation has three critical terms. The first term can be bounded using the RSS condition. Using $f(\bm{\theta}^{t})-f(\bm{\theta}^{*})\leq\langle\bm{g}^{t}_{S^{t}\cup S^{*}},\bm{\theta}^{t}-\bm{\theta}^{*}\rangle-\frac{\alpha}{2}\|\bm{\theta}^{t}-\bm{\theta}^{*}\|_{2}^{2}\leq\frac{1}{2\alpha}\|\bm{g}^{t}_{S^{t}\cup S^{*}}\|_{2}^{2}$ bounds the third term in (3). The second term is more interesting as in general elements of $\bm{g}^{t}_{\overline{S^{*}}}$ can be arbitrarily small. However, elements of $\bm{g}^{t}_{I^{t}\backslash(S^{t}\cup S^{*})}$ should be at least as large as $\bm{g}^{t}_{S^{*}\backslash S^{t+1}}$ as they are selected by hard-thresholding. Combining this insight with bounds for $\bm{g}^{t}_{S^{*}\backslash S^{t+1}}$ and with (3), we obtain the theorem. See Appendix A for a detailed proof. ∎

We now generalize our previous analysis to a projected gradient descent (PGD) method for low-rank matrix regression. Formally, we study the following problem:

The hard-thresholding projection step for low-rank matrices can be solved using SVD i.e.

where $W=U\Sigma V^{T}$ is the singular value decomposition of $W$. $U_{s},V_{s}$ are the top-$s$ singular vectors (left and right, respectively) of $W$ and $\Sigma_{s}$ is the diagonal matrix of the top-$s$ singular values of $W$. To proceed, we first note a property of the above projection similar to Lemma 1.

Let $W=U\Sigma V^{T}$ be the singular value decomposition of $W$. Now, $\|PM_{s}(W)-W\|_{F}^{2}=\sum_{i=s+1}^{|I^{t}|}\sigma_{i}^{2}=\|P_{s}(diag(\Sigma))-diag(\Sigma)\|_{2}^{2}$, where $\sigma_{1}\geq\dots\geq\sigma_{|I^{t}|}\geq 0$ are the singular values of $W$. Using Lemma 1, we get:

where the last step uses the von Neumann’s trace inequality ($Tr(A\cdot B)\leq\sum_{i}\sigma_{i}(A)\sigma_{i}(B)$). ∎

The following result for low-rank matrix regression immediately follows from Lemma 4.

Let $f$ have RSC and RSS parameters given by $L_{2s+s^{*}}(f)=L$ and $\alpha_{2s+s^{*}}(f)=\alpha$. Replace the projection operator $P_{s}$ in Algorithm 1 with its matrix counterpart $PM_{s}$ as defined in (5). Suppose we invoke it with $f,s\geq 32\left(\frac{L}{\alpha}\right)^{2}s^{*},\eta=\frac{2}{3L}$. Also let $W^{*}=\arg\min_{W,rank(W)\leq s^{*}}f(W)$. Then the $\tau$-th iterate of Algorithm 1, for $\tau=O(\frac{L}{\alpha}\cdot\log(\frac{f(W^{0})}{\epsilon})$ satisfies:

A proof progression similar to that of Theorem 1 suffices. The only changes that need to be made are: firstly Lemma 2 has to be invoked in place of Lemma 1. Secondly, in place of considering vectors restricted to a subset of coordinates viz. $\bm{\theta}_{S},\bm{g}^{t}_{I}$, we would need to consider matrices restricted to subspaces i.e. $W_{S}=U_{S}U_{S}^{T}W$ where $U_{S}$ is a set of singular vectors spanning the range-space of $S$. ∎

High Dimensional Statistical Estimation

This section elaborates on how the results of the previous section can be used to give guarantees for IHT-style techniques in a variety of statistical estimation problems. We will first present a generic convergence result and then specialize it to various settings. Suppose we have a sample of data points $Z_{1:n}$ and a loss function $\mathcal{L}(\bm{\theta};Z_{1:n})$ that depends on a parameter $\bm{\theta}$ and the sample. Then we can show the following result. (See Appendix B for a proof.)

Let $\bar{\bm{\theta}}$ be any $s^{*}$-sparse vector. Suppose $\mathcal{L}(\bm{\theta};Z_{1:n})$ is differentiable and satisfies RSC and RSS at sparsity level $s+s^{*}$ with parameters $\alpha_{s+s^{*}}$ and $L_{s+s^{*}}$ respectively, for $s\geq 32\left(\tfrac{L_{2s+s^{*}}}{\alpha_{2s+s^{*}}}\right)^{2}s^{*}$. Let $\bm{\theta}^{\tau}$ be the $\tau$-th iterate of Algorithm 1 for $\tau$ chosen as in Theorem 1 and $\varepsilon$ be the function value error incurred by Algorithm 1. Then we have

Note that the result does not require the loss function to be convex. This fact will be crucially used later. We now apply the above result to several statistical estimation scenarios.

where $\epsilon$ is the function value error incurred by Algorithm 1.

Fully-corrective Methods

In this section, we study a variety of “fully-corrective” methods. These methods keep the optimization objective fully minimized over the support of the current iterate. To this end, we first prove a fundamental theorem for fully-corrective methods that formalizes the intuition that for such methods, a large function value should imply a large gradient at any sparse $\bm{\theta}$ as well. This result is similar to Lemma 1 of but holds under RSC/RSS conditions (rather than the RIP condition as in ), as well as for the general loss functions.

Consider a function $f$ with RSC parameter given by $L_{2s+s^{*}}(f)=L$ and RSS parameter given by $\alpha_{2s+s^{*}}(f)=\alpha$. Let $\bm{\theta}^{*}=\arg\min_{\bm{\theta},\|\bm{\theta}\|_{0}\leq s^{*}}f(\bm{\theta})$ with $S^{*}=supp(\bm{\theta}^{*})$. Let $S^{t}\subseteq[p]$ be any subset of co-ordinates s.t. $|S^{t}|\leq s$. Let $\bm{\theta}^{t}=\arg\min_{\bm{\theta},supp(\bm{\theta})\subseteq S^{t}}f(\bm{\theta})$. Then, we have:

Here we will concentrate on a family of two-stage fully corrective methods that contains popular compressive sensing algorithms like CoSaMP and Subspace Pursuit (see Algorithm 2 for pseudocode). These algorithms have thus far been analyzed only under RIP conditions for the least squares objective. Using our analysis framework developed in the previous sections, we present a generic RSC/RSS-based analysis for general two-stage methods for arbitrary loss functions. Our analysis shall use the following key observation that the the hard thresholding step in two stage methods does not increase the objective function a lot.

Let $Z_{t}\subseteq[n]$ and $|Z_{t}|\leq q$. Let $\bm{\beta}^{t}=\arg\min_{\bm{\beta},supp(\bm{\beta})\subseteq Z_{t}}f(\bm{\beta})$ and $\widehat{\bm{\theta}}^{t}=P_{q}(\bm{\beta}^{t})$. Then, the following holds:

Let $\bm{v}^{t}=\nabla_{\bm{\theta}}f(\bm{\beta}^{t})$. Then, using the RSS property we get:

where $w$ is any vector such that $w_{\overline{Z_{t}}}=0$ and $\|w\|_{0}\leq s^{*}$. $\zeta_{1}$ follows by observing $\bm{v}^{t}_{Z_{t}}=0$ and by noting that $supp(\widehat{\bm{\theta}}^{t})\subseteq Z_{t}$. $\zeta_{2}$ follows by Lemma 1 and the fact that $\|w\|_{0}\leq s^{*}$. Now, using RSS property and the fact that $\nabla_{\bm{\theta}}f(\bm{\beta}^{t})=0$, we have:

The result now follows by combining (7) and (8). ∎

2 Partial Hard Thresholding Methods

We now study Partial Hard Thresholding methods (PHT), a family of fully-corrective iterative methods introduced by . This family is known to provide the best known RIP guarantees in the compressive sensing setting, but the proof is restricted to the RIP setting, and for the least-squares objective. An interesting member of this family is Orthogonal Matching Pursuit with Replacement (OMPR), which is also a Forward-Backward Greedy Selection method but performs one forward-backward step per iteration.

Let $f$, $s$ be supplied to Algorithm 3 and let the RSC and RSS parameters of $f$ be given by $L_{2s+s^{*}}(f)=L$ and $\alpha_{2s+s^{*}}(f)=\alpha$ respectively. Let $s\geq 4\left(\frac{L}{\alpha}\right)^{2}s^{*}$. Then, either $f(S^{t})=f(S^{*})$ or $|S^{t+1}\backslash S^{t}|\geq 1$. That is, at least one new element is added at each iteration of Algorithm 3.

where we have used the fact that $\bm{g}^{t}_{S^{t}}=0$ and $\bm{\theta}^{t}_{\overline{S^{t}}}=0$.

Using Lemma 3 and the above equation, we have:

The lemma now follows by observing that $\frac{|S^{t}\backslash S^{*}|}{|S^{*}\backslash S^{t}|}\geq\frac{s-s^{*}}{s^{*}}\geq\frac{\gamma^{2}}{\eta^{2}}\geq\frac{1}{\eta^{2}\alpha^{2}}$, by the choice of $s$. ∎

Let $f,s$ be supplied to Algorithm 3. Also, let the RSC and RSS parameters of $f$ be given by $\alpha_{2s+s^{*}}(f)=\alpha$ and $L_{2s+s^{*}}(f)=L$ respectively. Let $s\geq 4\left(\frac{L}{\alpha}\right)^{2}s^{*}$ and let $\eta=\frac{1}{2L}$. Then, the $\tau$-th iterate of Algorithm 3 satisfies:

where $\bm{\theta}^{*}=\arg\min_{\bm{\theta},\|\bm{\theta}\|_{0}\leq s^{*}}f(\bm{\theta})$.

Experiments

We conducted simulations on high dimensional sparse linear regression problems to verify our predictions. Our experiments demonstrate that hard thresholding and projected gradient techniques can not only offer recovery in stochastic setting, but offer much more scalable routines for the same.

Data: Our problem setting is identical to the one described in the previous section. We fixed a parameter vector $\bar{\bm{\theta}}$ by choosing $s^{*}$ random coordinates and setting them randomly to $\pm 1$ values. Data samples were generated as $Z_{i}=(X_{i},Y_{i})$ where $X_{i}\sim\mathcal{N}(0,I_{p})$ and $Y_{i}=\langle\bar{\bm{\theta}},X_{i}\rangle+\xi_{i}$ where $\xi_{i}\sim\mathcal{N}(0,\sigma^{2})$. We studied the effect of varying dimensionality $p$, sparsity $s^{*}$, sample size $n$ and label noise level $\sigma$ on the recovery properties of the various algorithms as well as their run times. We chose baseline values of $p=20000,s^{*}=100,\sigma=0.1,n=f_{o}\cdot s^{*}\log p$ where $f_{o}$ is the oversampling factor with default value $f_{o}=2$. Keeping all other quantities fixed, we varied one of the quantities and generated independent data samples for the experiments.

Algorithms: We studied a variety of hard-thresholding style algorithms including HTP , GraDeS (or IHT ), CoSaMP , OMPR and SP . We compared them with a standard implementation of the L1 projected scaled sub-gradient technique for the lasso problem and a greedy method FoBa for the same.

Evaluation Metrics: For the baseline noise level $\sigma=0.1$, we found that all the algorithms were able to recover the support set within an error of $2\%$. Consequently, our focus shifted to running times for these experiments. In the experiments where noise levels were varied, we recorded, for each method, the number of undiscovered support set elements.

Results: Figure1 describes the results of our experiments in graphical form. For sake of clarity we included only HTP, GraDeS, L1 and FoBa results in these graphs. Graphs for the other algorithms CoSaMP, SP and OMPR can be seen in the supplementary material. The graphs indicate that whereas hard thresholding techniques are equally effective as L1 and greedy techniques for recovery in noisy settings, as indicated by Figure1(a), the former can be much more efficient and scalable than the latter. For instance, as Figure1(b), for the base level of $p=20000$, HTP was $150\times$ faster than the L1 method. For higher values of $p$, the runtime gap widened to more than $350\times$. We also note that in both these cases, HTP actually offered exact support recovery whereas L1 was unable to recover $2$ and $4$ support elements respectively.

Although FoBa was faster than L1 on Figure1(b) experiments, it was still slower than HTP by $50\times$ and $90\times$ for $p=20000$ and $25000$ respectively. Moreover, due to its greedy and incremental nature, FoBa was found to suffer badly in settings with larger true sparsity levels. As Figure 1(c) indicates, for even moderate sparsity levels of $s^{*}=300$ and $500$, FoBa is $60-75\times$ slower than HTP. As mentioned before, the reason for this slowdown is the greedy approach followed by FoBa: whereas HTP took less than $5$ iterations to converge for these two problems, FoBa spend $300$ and $500$ iterations respectively. GraDeS was found to offer much lesser run times in comparison being slower than HTP by $30-40\times$ for larger values of $p$ and $2-5\times$ slower for larger values of $s^{*}$.

Experiments on badly conditioned problems. We also ran experiments to verify the performance of IHT algorithms in high condition number setting. Values of $p,s^{*}$ and $\sigma$ were kept at baseline levels. After selecting the optimal parameter vector $\bar{\bm{\theta}}$, we selected $s^{*}/2$ random coordinates from its support and $s^{*}/2$ random coordinates outside its support and constructed a covariance matrix with heavy correlations between these chosen coordinates. The condition number of the resulting matrix was close to $50$. Samples were drawn from this distribution and the recovery properties of the different IHT-style algorithms was observed as the projected sparsity levels $s$ were increased. Our results (see Figure 1(d)) corroborate our theoretical observation that these algorithms show a remarkable improvement in recovery properties for ill-conditioned problems with an enlarged projection size.

Discussion and Conclusions

In our work we studied iterative hard thresholding algorithms and showed that these techniques can offer global convergence guarantees for arbitrary, possibly non-convex, differentiable objective functions, which nevertheless satisfy Restricted Strong Convexity/Smoothness (RSC/RSM) conditions. Our results apply to a large family of algorithms that includes existing algorithms such as IHT, GraDeS, CoSaMP, SP and OMPR. Previously the analyses of these algorithms required stringent RIP conditions that did not allow the (restricted) condition number to be larger than universal constants specific to these algorithms.

Our basic insight was to relax this stringent requirement by running these iterative algorithms with an enlarged support size. We showed that guarantees for high-dimensional M-estimation follow seamlessly from our results by invoking results on RSC/RSM conditions that have already been established in the literature for a variety of statistical settings. Our theoretical results put hard thresholding methods on par with those based on convex relaxation or greedy algorithms. Our experimental results demonstrate that hard thresholding methods outperform convex relaxation and greedy methods in terms of running time, sometime by orders of magnitude, all the while offering competitive or better recovery properties.

Our results apply to sparsity and low rank structure, arguably two of the most commonly used structures in high dimensional statistical learning problems. In future work, it would be interesting to generalize our algorithms and their analyses to more general structures. A unified analysis for general structures will probably create interesting connections with existing unified frameworks such as those based on decomposability and atomic norms .

References

Appendix A Proofs for Section 3

Without loss of generality, assume that we have reordered coordinates such that $|\bm{z}_{1}|\geq|\bm{z}_{2}|\geq\ldots\geq|\bm{z}_{I}|$. Since the projection operator $P_{s}(\cdot)$ operates by selecting the largest elements by magnitude, we have $\bm{\theta}_{1}=\bm{z}_{1},\ldots,\bm{\theta}_{s}=\bm{z}_{s}$ and $\bm{\theta}_{s+1}=\bm{\theta}_{s+2}=\ldots=\bm{\theta}_{|I|}=0$.

Also define $\bm{\theta}^{\bm{z}}=P_{s^{*}}(\bm{z})$. By the above argument, we have $\bm{\theta}^{\bm{z}}_{1}=\bm{z}_{1},\ldots,\bm{\theta}^{\bm{z}}_{s^{*}}=\bm{z}_{s^{*}}$ and $\bm{\theta}^{\bm{z}}_{s^{*}+1}=\bm{\theta}^{\bm{z}}_{s^{*}+2}=\ldots=\bm{\theta}^{\bm{z}}_{|I|}=0$. Now we have

since the coordinates of $\bm{z}$ are arranged in decreasing order of magnitude. Combining the above with the observation that, due to the projection property $\|\bm{\theta}^{*}-\bm{z}\|\geq\|\bm{\theta}^{\bm{z}}-\bm{z}\|$, proves the result. ∎

Recall that $\bm{\theta}^{t+1}=P_{s}(\bm{\theta}^{t}-\frac{\eta^{\prime}}{L}\bm{g}^{t})$ where $\eta^{\prime}=\tfrac{2}{3}<1$. Let $S^{t}=supp(\bm{\theta}^{t})$, $S^{*}=supp(\bm{\theta}^{*})$, and $S^{t+1}=supp(\bm{\theta}^{t+1})$. Also, let $I^{t}=S^{*}\cup S^{t}\cup S^{t+1}$.

Now, using the RSS property and the fact that $supp(\bm{\theta}^{t})\subseteq I^{t}$ and $supp(\bm{\theta}^{t+1})\subseteq I^{t}$, we have:

As $supp(\bm{\theta}^{t})=S^{t}$, $supp(\bm{\theta}^{t+1})=S^{t+1}$ and $S^{t}\backslash S^{t+1},S^{t+1}$ are disjoint, we have:

where the equality $\zeta_{1}$ follows from the gradient step, i.e., $\bm{\theta}^{t+1}_{S^{t+1}}=\bm{\theta}^{t}_{S^{t+1}}-\frac{\eta^{\prime}}{L}\bm{g}^{t}_{S^{t+1}}$. The inequality $\zeta_{2}$ follows using the fact that $\bm{\theta}^{t+1}$ is obtained using hard thresholding and the fact that $|S^{t}\backslash S^{t+1}|=|S^{t+1}\backslash S^{t}|$, as follows:

The equality $\zeta_{3}$ follows from $\|\bm{g}^{t}_{S^{t+1}}\|_{2}^{2}=\|\bm{g}^{t}_{S^{t+1}\backslash S^{t}}\|_{2}^{2}+\|\bm{g}^{t}_{S^{t}\cap S^{t+1}}\|_{2}^{2}$.

Next, let us try to upper bound the first two terms on the right hand side above. Since $I^{t}\backslash(S^{t}\cup S^{*})=S^{t+1}\backslash(S^{t}\cup S^{*})\subseteq S^{t+1}$, we have $\bm{\theta}^{t+1}_{I^{t}\backslash(S^{t}\cup S^{*})}=\bm{\theta}^{t}_{I^{t}\backslash(S^{t}\cup S^{*})}-\frac{\eta^{\prime}}{L}\bm{g}^{t}_{I^{t}\backslash(S^{t}\cup S^{*})}$. However, as $\bm{\theta}^{t}_{I^{t}\backslash S^{t}}=0$, we actually have $\bm{\theta}^{t+1}_{I^{t}\backslash(S^{t}\cup S^{*})}=-\frac{\eta^{\prime}}{L}\bm{g}^{t}_{I^{t}\backslash(S^{t}\cup S^{*})}$. Now let us choose a set $R\subseteq S^{t}\backslash S^{t+1}$ such that $|R|=|S^{t+1}\backslash(S^{t}\cup S^{*})|$. Such a choice is possible since $|S^{t+1}\backslash(S^{t}\cup S^{*})|=|S^{t}\backslash S^{t+1}|-|(S^{t+1}\cap S^{*})\backslash S^{t}|$ (which itself is a consequence of the fact that $|S^{t+1}|=|S^{t}|$). Moreover, since $\bm{\theta}^{t+1}$ is obtained by hard-thresholding $\left(\bm{\theta}^{t}-\frac{\eta^{\prime}}{L}\bm{g}^{t}\right)$, for any choice of $R$ made above, we have:

Using above equation, and the fact that $\bm{\theta}^{t+1}_{R}=0$ (since $R\subseteq\overline{S^{t+1}}$), we have:

We can bound the size of $I^{t}\backslash R$ as $|I^{t}\backslash R|\leq|S^{t+1}|+|(S^{t}\backslash S^{t+1})\backslash R|+|S^{*}|\leq s+|(S^{t+1}\cap S^{*})\backslash S^{t}|+s^{*}\leq s+2s^{*}.$ Also, since $S^{t+1}\subseteq(I^{t}\backslash R)$, we have $\bm{\theta}^{t+1}_{I^{t}\backslash R}=P_{s}(\bm{\theta}^{t}_{I^{t}\backslash R}-\frac{\eta^{\prime}}{L}\bm{g}^{t}_{I^{t}\backslash R})$.

Using the above observation with (16) and Lemma 1, we get:

where the inequality $\zeta_{1}$ follows by $|I^{t}\backslash R|\leq s+2s^{*}$ as shown earlier and the observation that $\frac{x-a}{x-b}$ is a positive and increasing function on the interval $x\geq a$ if $a\geq b\geq 0$. Note that since we have $S^{t+1}\subseteq(I^{t}\backslash R)$, we get $|I^{t}\backslash R|\geq s$. The inequality $\zeta_{2}$ follows by using RSC.

Using (14), (17), and using $S^{t+1}\backslash(S^{t}\cup S^{*})\subseteq(S^{t+1}\cup S^{t})$, we get:

We now set $\eta^{\prime}=2/3$ as per our earlier choice and set $s=32\left(\frac{L}{\alpha}\right)^{2}s^{*}$, so that we have $\frac{2s^{*}}{s+s^{*}}\leq\frac{\alpha^{2}}{16L(L-\eta^{\prime}\alpha)}$. Since $L\geq\alpha$, we also have $\frac{\alpha^{2}}{16L(L-\eta^{\prime}\alpha)}\leq\frac{3}{16}$. Using these inequalities, we now rearrange the terms in (18) above.

Splitting $\|\bm{g}^{t}_{I^{t}}\|_{2}^{2}=\|\bm{g}^{t}_{S^{t}\cup S^{*}}\|_{2}^{2}+\|\bm{g}^{t}_{S^{t+1}\backslash(S^{t}\cup S^{*})}\|_{2}^{2}$ gives us

where the last inequality above follows using Lemma 6. The result now follows by observing that $\frac{2s^{*}}{s+s^{*}}\geq 0$. ∎

where the first inequality above follows from (21). ∎

Appendix B Proofs for Section 4

Let $\bm{\theta}^{*}$ be the empirical loss minimizer over the set of $s$-sparse vectors. Then invoking Theorem 1 with $f=\mathcal{L}(\cdot;Z_{1:n})$, we get

where the 2nd inequality is by definition of $\bm{\theta}^{*}$ and 3rd is by RSC (since $\bm{\theta}^{*},\bm{\theta}^{\tau}$ are $s^{*},s$ sparse). Duality gives us the upper bound

Combining the last two inequalities and rearranging gives a quadratic inequality in $\|\bar{\bm{\theta}}-\bm{\theta}^{\tau}\|_{2}$:

Appendix C Proofs for Section 5

We will start by proving a more general result of which the claimed result will be a corollary. More specifically, we shall prove that for any $\gamma\geq\frac{1}{\alpha}$, we have

Setting $\gamma=\frac{1}{\alpha}$ will yield the claimed result. It is easy to see that the following inequality holds trivially since $\gamma\geq\frac{1}{\alpha}$

For the second inequality, we first use the RSC condition to obtain:

Now let $MD_{t}=S^{*}\backslash S^{t}$ be the set of true support elements missing from $\bm{\theta}^{t}$ and $FA_{t}=S^{t}\backslash S^{*}$ be the set of incorrect elements included in the support of $\bm{\theta}^{t}$. Since $\bm{\theta}^{t}$ is obtained by a “fully corrective” process (recall $\bm{\theta}^{t}=\arg\min_{\bm{\theta},supp(\bm{\theta})\subseteq S^{t}}f(\bm{\theta})$), we have $\bm{g}^{t}_{S^{t}}=\bm{0}$. Thus $\langle\bm{\theta}^{*}-\bm{\theta}^{t},\bm{g}^{t}\rangle=\langle\bm{\theta}^{*}_{MD_{t}},\bm{g}^{t}_{MD_{t}}\rangle$.

Putting this into the above expansion gives

We now present some simple inequalities that will help us get our desired bounds. Firstly, we have

since the first expression is a norm. Next, since $MD_{t}\cap FA_{t}=\emptyset$, we have

We finish off the proof by noticing that since $\bm{g}^{t}_{S^{t}}=\bm{0}$, we have $\|\bm{g}^{t}_{MD_{t}}\|_{2}^{2}=\|\bm{g}^{t}_{S^{t}\cup S^{*}}\|_{2}^{2}$ ∎

Let $\bm{z}^{t}_{S^{t}}=\bm{\theta}^{t}_{S^{t}}$, $\bm{z}^{t}_{Z^{t}\backslash S^{t}}=-\frac{1}{L}\bm{g}^{t}_{Z^{t}\backslash S^{t}}$, and $\bm{z}^{t}_{\overline{Z^{t}}}=0$.

Now, using Lemma 4 and (27) along with $f(\bm{\theta}^{t+1})\leq f(\widetilde{\bm{\theta}}^{t})$ and $f(\bm{\beta}^{t})\leq f(\bm{z}^{t})$, we have:

where $\zeta_{1}$ follows by observing that $\bm{g}^{t}_{S^{t}}=0$ and $\bm{v}^{t}_{S^{t+1}\backslash S^{t}}=-\eta\bm{g}^{t}_{S^{t+1}\backslash S^{t}}$. $\zeta_{2}$ follows by the property of PHT operator which ensures that each element of $\bm{v}^{t}_{S^{t+1}\backslash S^{t}}$ is bigger than $\bm{\theta}^{t}_{S^{t}\backslash S^{t+1}}$ and by using $|S^{t+1}\backslash S^{t}|=|S^{t}\backslash S^{t+1}|$. $\zeta_{3}$ follows by using $\bm{v}^{t}_{S^{t+1}\backslash S^{t}}=-\eta\bm{g}^{t}_{S^{t+1}\backslash S^{t}}$.

Similar to the analysis given in , we divide the analysis in three mutually exclusive cases. The lemma then follows by combining (29) with the case-by-case analyses below and observing that $f(\bm{\theta}^{t+1})\leq f(\bm{v}^{t})$ because of the fully corrective step.

Using $\bm{z}^{t}_{S^{*}\backslash S^{t}}=-\eta\bm{g}^{t}_{S^{*}\backslash S^{t}}$, $\bm{z}^{t}_{S^{t}}=\bm{\theta}^{t}_{S^{t}}$ and $(S^{t+1}\cap S^{t})\backslash S^{*}\subseteq S^{t}\backslash S^{*}$, we have:

Using the above equation and Lemma 3, we have:

Using $s\geq 4\left(\frac{L}{\alpha}\right)^{2}s^{*}$, we have:

Case 3: $|S^{t+1}\backslash S^{t}|\geq|S^{*}\backslash S^{t}|$. Now, as $\bm{v}^{t}_{S^{t+1}\backslash S^{t}}$ is obtained by selecting $|S^{t+1}\backslash S^{t}|$ largest elements of $\bm{z}^{t}_{\overline{S^{t}}\cup bot^{t}}$. Hence, using Lemma 3, we have:

The lemma now follows by combining (29), (31), (35), and (36) ∎

Appendix D Supplementary Experimental Results

Below we present plots that were not included in the main text.