Fast Iterative Hard Thresholding Methods with Pruning Gradient Computations

We sincerely appreciate the reviewer’s positive feedback on our paper. As mentioned in the reviewer's comment, previous methods for accelerating IHT have focused on reducing the number of iterations. In contrast, our approach accelerates IHT by reducing the computation cost per iteration. Accelerating IHT is challenging due to its non-convexity, resulting from $\ell_0$ cardinality constraints in Problem (1), and there are relatively few acceleration methods for IHT compared to those for $\ell_1$ regularization, such as lasso. Therefore, we hope that the proposed method will offer a new direction for speeding up IHT.

We also thank the reviewer for the valuable suggestions on improving the presentation. We will revise the paper accordingly.

Below are our responses to the questions:

Would it be possible to combine your technique with techniques that aim to reduce the number of iterations? For instance, the techniques described in references 4 and 26, or 8 and 20?

Yes, it is relatively straightforward to combine these methods with our proposed method. This is because these methods utilize a linear combination of parameters to reduce the number of iterations. For example, the acceleration methods in [8] and [20] can be integrated with our method as follows:

1. Compute the linear combination of parameters at the $t$-th iteration $\phi^{t}$, by following Equation (21) in [8] or Equation (19) in [20].

2. Set $\theta^{t}=\phi ^{t}$. Then, we can compute the upper and lower bounds using Equations (4) and (5) in our paper, respectively.

3. Perform Algorithm 4 in our paper based on these upper and lower bounds.

We are currently working on combining these acceleration methods with our approach and appreciate the valuable feedback.

Why did you not include the methods from references 4, 26, 8, and 20 in your experiments? Are they superseded by references 2 and 19?

We sincerely thank the reviewer for the insightful comment. References [4, 26] utilize acceleration techniques known as double over-relaxation. However, the theoretical guarantees for these methods are limited. Specifically, the method in [26] lacks theoretical guarantee, and the method in [4] is guaranteed to converge only when the acceleration technique is applied under conditions where the objective function decreases. In contrast, [8, 20] provide convergence proofs when using Nesterov’s acceleration [24] and squared loss. Additionally, [19] examines the properties of the momentum step in Nesterov’s acceleration. Since the algorithms in [8, 20, 19] are almost the same, we consider [8, 20] to be superseded by [19] in our paper.

Unlike the above methods, [2] adds squared $\ell_2$ regularization to Problem (1), and is completely different from the above papers. Therefore, we have added it to the baseline in our experiment.

We would like to incorporate the above explanation into the related work section.

We sincerely appreciate the reviewer's thoughtful review of our paper. Below are our answers to the questions and weaknesses raised by the reviewer.

I find the terminology "pruning" a bit confusing. My understanding is that it refers to setting a particular entry $z_{j}^{t}$ to zero depending on whether or not $\bar{z}^{t}_{j}$ is less than a certain threshold. Is this correct? If yes, then maybe "thresholding" is a better terminology?

A1. Yes, that is correct. We chose the term "pruning" instead of "thresholding" to distinguish our procedure from the thresholding procedure of the hard thresholding operator $P_{k}$​. The thresholding of $P_{k}$ sets unnecessary $z_{j}^{t}$ to zero after computing the gradient and updating all the parameters. In contrast, our method's pruning sets unnecessary $z_{j}^{t}$ to zero before computing the gradient and updating the parameters.

By "pruning is safe", do you mean that the pruning step does not compromise the accuracy of IHT at all? Can you prove this more explicitly?

A2. Yes, our pruning method is theoretically guaranteed to maintain the same accuracy as the original IHT. This is demonstrated in Theorem 2, with detailed proof provided in the appendix. Below, we outline the proof:

1. The original IHT selects the top-$k$ parameters in absolute values. In contrast, our method uses a candidate set $\mathbb{D}^{t}$ with cardinality $k$, as described in Definition 2, and updates $\mathbb{D}^{t}$ to include the top-$k$ parameters during optimization.

2. At the beginning of each iteration, our method initializes $\mathbb{D}^{t}$ to include the indices of the current top-$k$ parameters (line 103). For the parameters corresponding to $\mathbb{D}^{t}$, we compute their gradients to obtain the exact $z_{i}^{t}$ for $i$ in $\mathbb{D}^{t}$ (line 13 in Algorithm 4). We then set $z_{i_{\rm min}}^{t}$ as the minimum $z_{i}^{t}$​ in $\mathbb{D}^{t}$, as shown in Lemma 3.

3. Next, we consider the upper bound of $z_{j}^{t}$​ for parameters not included in $\mathbb{D}^{t}$. This upper bound $\bar{z}^{t}\_{j}$ is given by Definition 1 and Lemma 1. If $\bar{z}^{t}\_{j} \leq z\_{i\_{\rm min}}^{t}$ holds,​ $z_{j}^{t}\leq z\_{i\_{\rm min}}^{t}$ also holds.​ Therefore, the $j$-th parameter must not be included in $\mathbb{D}^{t}$. In this case, the $j$-th parameter​ can be safely pruned without computing the gradient (Lemma 3). For parameters that do not satisfy the condition, our method exactly computes $z_{j}^{t}$​ as shown in line 8 of Algorithm 2. Consequently, $\mathbb{D}^{t}$ is exactly the same as the top-$k$ parameters of the original IHT (Lemma 4).

As described above, since our method prunes parameters that must be zero, pruning the corresponding gradient computations does not affect the final top-$k$ parameters. Additionally, our method not only achieves the same final accuracy as the original IHT but also ensures that the parameters during optimization match perfectly with those of the original IHT. We would like to highlight that this property is not found in previous methods such as AccIHT[19].

I think X^TX takes O(n^2m) computational cost instead of O(mn)?

A3. In the part of $X^T(X\theta^{t}-y)$ in Equation (2), we first compute $h=X\theta^{t}-y$ at $\mathcal{O}(mn)$ time, then compute $X^T h$ at $\mathcal{O}(mn)$ time. Consequently, Equation (2) requires $\mathcal{O}(mn)$ time. We would like to incorporate this explanation into our paper.

Did you mean to use f(\theta) instead of 1/2||y-X\theta||_2^2 in (1)?

A4. Yes, for the sake of simplicity, we use the expression as stated in line 58.

I'd suggest the author add one or two sentences in the final Section to discuss the limitations of this work.

We thank the reviewer for this comment. The limitation of our method is that the speedup factor decreases for large $k$ as described in line 270-272. However, as demonstrated in Theorem 1, the worst time complexity of our method is the same as that of the original IHT. We will incorporate this explanation into the paper.

To weaknesses:

It is unclear to me why the pruning strategy is defined the way it is.

Kindly refer to our response in A2, which details the pruning strategy of our method and explains why this pruning strategy is safe. The point is that our method prunes unnecessary parameters before computing gradients by using upper bounds, which have small time complexity.

Consider adding more literature review: it is not clear how big of a gap exists in the literature that this work is trying to bridge. Section 4 Related Work feels out of place. Could consider moving this to the beginning of the paper.

We would like to move Section 4 to the beginning of the paper. As explained in lines 245-246, while previous acceleration methods for IHT have reduced the number of iterations, there has not yet been a method to accelerate IHT by reducing the computation cost per iteration. This paper aims to fill this gap based on the pruning strategy.

A general suggestion: consider adding more explanation on the rationale behind each technical definition/ result and why it is defined/ stated in that particular way. For instance, it may be helpful to mention that Lemma 1 is derived from the triangle inequality + Holders inequality.

We plan to explain the sketch of the proof and the reason for defining it before and after theorems and definitions, respectively.

I think some definitions are stated in the form of lemmas which they should not be, e.g. Lemmas 3 and 8. In general, I think some lemmas are unnecessary or can be combined.

We will revise the manuscript to write Lemmas 3 and 8 as normal paragraphs or definitions.

Presentation is a bit long-winded at places, e.g. third paragraph

We intend to shorten redundant expressions, such as those found in the third paragraph.

To minor comments:

We plan to revise the paper by incorporating these helpful suggestions, including the notations.

We appreciate the reviewer’s constructive comments on our paper. To begin, we address the questions raised by the reviewer below.

1. Line 42-45 is duplicated with the abstract.

We thank the reviewer for this comment. We will revise those sentences accordingly.

2. Definition 2 is less informative. It generally says nothing about the construction of the candidate set $\mathbb{D}^{t}$.

We are grateful to the reviewer for this comment. We will revise the manuscript to write it as a paragraph rather than a definition.

3. The performance seems highly sensitive to the selection of $t^*$.

We greatly appreciate the reviewer for highlighting this critical point. The problem of sensitivity of $t^*$ is solved by the automatic determination of $t^*$ as described in lines 194-215. Specifically, the proposed method automatically adjusts $t^*$ during optimization to enhance performance, thereby eliminating the need to select $t^*$. According to Lemma 11, the error bounds of the upper and lower bounds may increase depending on the value of $t^*$, potentially reducing the pruning rate described in Definition 4 and consequently affecting performance. However, our method monitors the pruning rate and automatically adjusts $t^*$ to prevent a decrease in the pruning rate. Specifically, if the pruning rate decreases during optimization, the interval until the next $t^*$ is reduced following Equation (8). From Lemma 11, we know that the error bound of the upper and lower bounds tends to decrease when the interval until the next $t^*$ is small. Thus, this automatic adjustment of $t^*$ helps maintain the pruning rate.

4. The proposed method seems only suitable when the learning rate is small and the parameters update gradually. When a larger learning rate is used or the momentum is introduced as in the related work (Section 4), the pruning may no longer take effect.

Although the pruning rate might decrease due to the introduction of a large learning rate or momentum, the proposed method can adjust the pruning rate by automatically determining $t^*$ during optimization, as mentioned in our previous response.

To demonstrate this, we conducted an additional experiment with a learning rate on the gisette dataset. The results are shown in PDF file in the global response. In the experiment, we increased the learning rate to 10 times that used in the experiment of Figure 1 (a) in our paper. The results exhibit a similar trend to Figure 1(a). Our method could accelerate IHT even with a large learning rate. We would like to highlight that all the objective values of AccIHT, which utilizes the momentum, diverged with the setting of $k=1280$, and we could not evaluate the processing times. This suggests that the method with momentum tends to fail when a large learning rate is used in our setting. In contrast, the methods without momentum, including our proposed method, converged in all settings. We would like to incorporate this finding into our paper.

Finally, we would like to mention that even if the pruning rate of the proposed method is low, its worst time complexity remains the same as the original IHT (lines 218-224).

To Weaknesses:

The proposed upper/lower bounds seem to be very restricted to the structure of the sparse linear regression task. And the proposed method does not work when general convex loss functions are considered.

As the reviewer pointed out, the proposed method is specialized for sparse linear regression with $\ell_0$ cardinality constraints. We agree that extending our method to general convex loss functions is important, and we plan to address this in future work. We would like to mention that the sparse linear regression with $\ell_0$ cardinality constraints is widely used in various fields, such as feature selection, sparse coding, dictionary learning, and compressed sensing, as described in lines 13--15. Although we performed a feature selection task in our experiment, we conducted an additional experiment on an image compressive sensing recovery task in response to Reviewer Hc1f. Our method completed the experiment in approximately 10 hours, while all other methods could not complete it within the author response period (4 days). We think that our method can contribute to reduce the times required for many tasks using IHT.

For the sparse linear regression task, there are already efficient algorithms. For example, the cordinate descent algorithm uses in glmnet, which only needs to compute one dimension of marginal gradient in per iteration. The comparation beyond the IHT-type algorithm is ignored in the paper.

While it is true that there are many algorithms for accelerating sparse linear regression, the non-convexity due to $\ell_0$ cardinality constraints in Problem (1) of our paper limits the number of available acceleration methods, which is a crucial motivation for proposing our method. For instance, the coordinate descent method of glmnet, mentioned by the reviewer, is specialized for convex optimization with $\ell_1$ regularization, such as lasso and elastic net, and therefore cannot be used for the non-convex optimization with $\ell_0$ cardinality constraints in Problem (1). Similarly, screening [a] and working set algorithm [b] also accelerate sparse linear regression tasks but rely on the computation of a duality gap based on convex optimization and are not applicable to Problem (1). Although lines 245-252 of our paper contain a similar explanation, we believe that incorporating the above discussion will enhance the paper's readability, and we would like to reflect this in the revised manuscript.

[a] O. Fercoq, A. Gramfort, J. Salmon, "Mind the duality gap: safer rules for the Lasso", ICML, 2015.

[b] M. Massias, J. Salmon, and A. Gramfort, "Celer: a Fast Solver for the Lasso with Dual Extrapolation", ICML, 2018.

To Limitations:

We intend to incorporate the above discussion into our paper. We sincerely appreciate the reviewer’s helpful comments.

The authors' responses sufficiently address my questions. In response, I would like to raise my score. Overall, it is an interesting work from the computation aspect of the sparse regression problem.

We extend our sincere gratitude to the reviewer for working on our paper. We address the weaknesses and questions raised by the reviewer below.

The importance of the work seems to be not clearly conveyed.

A1. As mentioned in lines 13—15, IHT is widely used in various fields, such as feature selection, sparse coding, dictionary learning, and compressed sensing. The original paper of IHT [6] has thereby been cited more than 2800 times. However, if design matrix $X$ is large, IHT suffers from high computation costs of gradient computations as described in 25—32. For instance, in our experiment of feature selection task, IHT requires over 24 hours for the total processing time, including hyperparameter search on the ledger and epsilon datasets, as shown in Figure 1 (c) and (e). Furthermore, in our answer of A3 below, we conducted an additional experiment on the image compressive sensing recovery task, and the baselines could not complete the task within 4 days. Our method completes all experiments, including additional experiments, within 24 hours, and we believe it can contribute to reduce the processing times required for many tasks using IHT. 

From a technical point of view, although accelerating IHT is the critical issue as described above, the number of approaches to accelerate IHT has been limited due to the non-convexity caused by $\ell_0$ cardinality constraints in Problem (1). As mentioned in lines 245-246, the previous approaches have reduced the number of iterations to accelerate IHT. In contrast, our approach reduces the computation cost per iteration by pruning unnecessary computations. Therefore, our method will provide a new direction for the acceleration of IHT. As illustrated in Figure 1, our approach significantly enhances the speedup factor compared to baselines using the previous approaches. We would like to incorporate this explanation into the beginning of our paper.

The evaluation is also limited. While providing significant speedup, for the scale of the problem considered, the wall time of the baselines still seems to be within a tolerable range.

A2. We thank the reviewer for highlighting that issue. Regarding the comparison of processing times in Figure 1, the processing time for each $k$ and hyperparameter falls within a tolerable range. However, in our experiments, the total processing time for each baseline on the ledger and epsilon datasets exceeded 24 hours. While the acceptable range of processing time is subjective, many practitioners would likely find more than a day impractical. In contrast, our method took less than 24 hours, even for the ledger and epsilon datasets. We would like highlight that we conducted an additional experiment on the image compressive sensing recovery task, and the baselines could not complete the task within 4 days as we will explain in our answer of A3 below.

Is it possible to include more experiments from the other use cases mentioned in the paper such as sparse coding, dictionary learning, or compressed sensing, preferably of a larger scale?

A3. We are grateful for the constructive feedback. We conducted an additional experiment on the image compressive sensing recovery task by following [a]. Although our experimental setting was almost the same as in Section 5 of [a], we used a large image with 2728$\times$2728 resolution. Consequently, the number of image blocks was 348,843. We investigated $k = [1, 8, 16, 24, 32]$ and evaluated the total processing times, including hyperparameter search. Our method completed the experiment in approximately 10 hours. Unfortunately, all other methods could not complete the experiment within the author's response period (4 days). The strong baselines of AccIHT and RegIHT also could not finish the experiment because they require additional hyperparameters for the weight step size and the momentum step size, respectively. We would like to highlight that our method does not require any additional hyperparameters, in contrast to AccIHT and RegIHT.

[a] J. Zhang, C. Zhao, D. Zhao, and W. Gao, “Image compressive sensing recovery using adaptively learned sparsifying basis via L0 minimization”. Signal Process. 103: 114-126, 2014.