Thank you for the reviews! We'd like to address your concerns as follows.

• On the scope of the paper. We focus mainly on the case where $\mathrm{IE}=2$ because $\mathrm{IE}=2$ is special in the following ways: (1) the IE of even functions is at least $2$, (2) most link functions/activations have IE at most $2$, and (2') the IE of many functions, including all finite degree polynomials, can be reduced to $2$ by a proper label transformation.

  Moreover, we believe our hierarchical feature learning analysis can potentially shed some light on the much harder non-orthogonal version of this problem. In the non-orthogonal case, the minimizers of the correlation loss, truncated at any finite order, are not aligned with the ground truth directions. However, as the order increases, they will become localized around the ground-truth directions, as the influence decay faster (for example, the "infinite-order" terms are nonzero only if the learner neuron is fully aligned with a ground truth direction). Intuitively, one might be able to use the lower order terms first to obtain warm starts, and then gradually move to higher order terms to move toward the ground truth directions. Formalizing this or finding conditions under which this is possible is challenging, and the orthogonal $h_2 + h_L + \cdots$ setting considered in this paper provides a simple example where this type of argument is doable.

  Finally, our results can be (partially) extended to other settings as well (polynomially bounded link functions, including linear terms, non-isotropic, and/or misspecified). You could check the responses to Reviewer 4V3d, x9Fa and uJL9 for details.

• $P d P^{L-2}$ vs $d P^{L-1}$. We write $P d P^{L-2}$ instead of $d P^{L-1}$ to emphasize the different sources of $P$ and $P^{L-1}$. The first $P$ comes from the fact we have $P$ different directions, while the second part $d P^{L-1}$ comes from learning a single-index model in a $P$-dimensional subspace with noises scaling with $d$.

Thank you for the reviews! We'd like to address your concerns as follows.

• Recovering the ground truth directions. We show in Lemma B.8 that at the end of Stage 1.2, every ground truth direction will be recovered by a neuron, in the sense that the correlation is at least $1 - \epsilon$. The only thing Stage 2 does is picking out these neurons, or, more accurately, find a combination of the neurons that is at least as good as the one consists of these neurons. If one really wants to recover the directions explicitly, one may add a $l^1$ regularizer $\lambda \sum_{k=1}^m |a_k|$ and do lasso regression, which returns a sparse solution. Alternatively, this step can also be done with $l^1$-regularized SGD via a local convergence argument similar to the one in [Chi20].

• On the scope of the paper. Our results can be extended to the misspecified case, where the link functions along different directions are different and unknown to the learner, as long as the activation of the learner has some universal approximation property, and the target link functions still have isotropic 2nd order terms, and the condition number of the $L$-th order terms is constant. This is similar to the situation in [BBSS22]. Recovering the subspace/directions does not require the activation to be exactly the same as the target link function, and after learning the directions, fitting the target link function can be done with a one-dimensional random feature argument.

  Our results can be (partially) extended to other settings as well (polynomially bounded link functions, including linear terms, and/or non-isotropic). You could check the responses to Reviewer 4V3d and x9Fa for details.

• Issue of IE/GE. We are not claiming that the IE/GE theory is wrong. Instead, our point is that the vanilla IE/GE theory, which considers only the lowest order, is insufficient to characterize the richer behavior of multi-index models. For example, we cannot distinguish between $h_2$ and $h_2 + h_4$ when looking only at the lowest order. On the other hand, SGD's behaviors for these target functions are different, as directional convergence is possible only in the 2nd case. We propose to look beyond the lowest order and show that, in the setting considered in this work, this leads to a sharper characterization.

  In addition, we believe this point of view might shed some light on the much harder non-orthogonal setting. In this case, the minimizers of the correlation loss, truncated at any finite order, are not aligned with the ground truth directions. However, as the order increases, they will become localized around the ground-truth directions, as the influence decay faster (for example, the "infinite-order" terms are nonzero only if the learner neuron is fully aligned with a ground truth direction). Intuitively, one might be able to use the lower order terms first to obtain warm starts, and then gradually move to higher order terms to move toward the ground truth directions. Formalizing this or finding conditions under which this is possible is challenging, and the orthogonal $h_2 + h_L + \cdots$ setting considered in this paper provides a simple example where this type of argument is doable.

[Chi20] Lenaic Chizat. Sparse Optimization on Measures with Over-parameterized Gradient Descent. 2020.

[BBSS22] Alberto Bietti, Joan Bruna, Clayton Sanford, Min Jae Song. Learning Single-Index Models with Shallow Neural Networks. 2022.

Thank you for your feedback! We'd like to address your concerns as follows.

• The condition on $P$. The condition $P \ge \log^C d$ is mostly an artifact of our analysis. First, we want to ensure certain concentration properties hold in the relevant subspace (such as the norm in the relevant subspace). We do not want to explicitly track the failure probability of these events, so we assume $P \ge \log^C d$, so that the failure is inverse super-polynomially small. In addition, this also creates a cleaner separation between the two stages. For example, at the end of the first stage, the 2nd order terms have order $\| v\_{\le P} \|\_2^2 \approx 1$, while the higher order terms have order $\| v\_{\le P} \|\_L^L$. The higher order terms are much smaller than the 2nd order terms when $P \gg 1$, and therefore can be safely ignored. We believe that this condition can be removed at the cost of a messier proof.

• Comparison with [DLS22]. First, we wish to stress that our sample complexity of the second stage is also tighter than the one in [DLS22]. They use a random feature argument in Stage 2, which leads to the sample complexity $d r^q$, where $q$ is the degree of the polynomial. In our case, we show that SGD can do feature learning, and the sample complexity is $d r^L$, where $L$ is the information exponent, which can be much smaller than the degree. In addition, our results can be extended to any link functions with a polynomially bounded derivative (see the response to Reviewer 4V3d), while their results only deal with polynomials. Moreover, their width complexity is $P^q$ and ours is $P \log P$ (cf. the statement of Theorem 2.1 in Appendix D). In some sense, [DLS22] shows that SGD can learn the subspace using $d^2$ samples, while our results show that, at least in a more restricted but still non-trivial setting, SGD can automatically do hierarchical feature learning, which yields sample/width complexities that are (almost) tight, in the sense that they match the complexities of separately learning the subspace, and then independently learning $P$ single-index models in the relevant subspace.

• Relaxing the orthonormality condition. First, we believe that our results can be extended to the near orthogonal case. However, we expect the non-orthogonal case to be hard, as it is known that non-orthogonal tensor decomposition is hard in general, whence we cannot expect SGD to learn the directions with polynomial samples/iterations, even when the relevant subspace is given.

  Things are trickier if we consider only the 2nd order terms (which is always "orthogonal" because of the spectral theorem), but do not assume isotropicity. In this case, correlation loss will fail, since all neurons will simply converge to the largest eigendirection (similar to what will happen in the matrix power iteration). This can be fixed by using the MSE loss and an initialization whose scale matches the target's. Our analysis does not directly apply to this case due to the interaction between the learner neurons, but the behaviors should be similar. That is, the 2nd order terms drive the distribution of learner neurons to (close to) the Gaussian in the relevant subspace, with the correct covariance matrix, and then the higher order terms become non-negligible and allow the neurons to recover the exact directions.

I thank the authors for their response and I do not have further questions. I would like to maintain my score.

Thank you for the reviews and suggestions! We'd like to answer your questions as follows.

• Width complexity. We realized that the width can be as small as $O(P \log P)$ after submitting the main text. This is tight up to the logarithmic term and matches the coupon collector bound. We have included this in the appendix. See Lemma A.7 and the statement of the main theorem in Appendix D.

• Relaxing Assumption 1(c). Thank you this question! It makes us realize that polynomial boundedness is sufficient. Assumption 1(c) is only used to get a tail bound on the per-sample gradient, which amounts to controlling the tail of $f_* = \sum_p \phi(v_p^* \cdot x)$ (times the derivative of a single neuron). We can use the polynomial boundedness to bound the $L^p$ norm of each $\phi(v_p^* \cdot x)$ and then show that they are sub-Weibull. Since $v_p^\*$ are orthogonal, $\phi(v_p^\* \cdot x)$ are independent, and therefore, the sum of them is still sub-Weibull and this leads to the desired tail bound. We will add this into the revision!

• The diverging width assumption. The condition $P \ge \log^C d$ is mostly an artifact of our analysis. First, we want to ensure certain concentration properties hold in the relevant subspace (such as the norm in the relevant subspace). We do not want to explicitly track the failure probability of these events, so we assume $P \ge \log^C d$, so that the failure is inverse super-polynomially small. In addition, this also creates a cleaner separation between the two stages. For example, at the end of the first stage, the 2nd order terms have order $\| v\_{\le P} \|\_2^2 \approx 1$, while the higher order terms have order $\| v\_{\le P} \|\_L^L$. The higher order terms are much smaller than the 2nd order terms when $P \gg 1$, and therefore can be safely ignored. We believe that this condition can be removed at the cost of a messier proof.

• Allowing linear terms. First, we wish to stress that we consider the case IE=2 since it is the lowest possible when the link function is even, or equivalently, both $v_k^\*$ and $-v_k^\*$ appear in the teacher network. Now, back to the question of the influence of linear components. You are right that there will be a spike (which is the all 1 direction if we assume w.l.o.g. that the targets are $\{ e_p \}_p$) and the all neurons will be attracted by that direction.

  To fix this issue, the easiest approach is to first train one neuron to fit that direction to remove it. Due to the linearity, learning the linear term is roughly equivalent to learning a single-index model with IE=1, and it costs only $O(d)$ samples.

  The more natural but technically more challenging approach is to use MSE loss. Similar to the 2nd vs 4th order terms argument, one should be able to show that the first order terms will get learned before the higher terms become large. The tricky thing here is that after learning the first order term, the learner distribution will be shifted to be centered at the all-one point (assuming $v_k^* = e_p$ for simplicity), instead of one single neuron learning that direction with all other neurons staying close to the initialization, as the gradient of the linear term depends only on the average of the neurons and not on any individual neuron. The analysis become tricky after that, since we can no longer ignore the interaction between learner neurons, though in principle, a finite-width mean-field type analysis similar to the one in [LMZ20] can potentially be used to analyze the dynamics of neurons in those later stages.