Sparse high-dimensional functions have arisen as a rich framework to study the behavior of gradient-descent methods using shallow neural networks, showcasing their ability to perform feature learning beyond linear models.
Amongst those functions, the simplest are single-index models $f(x) = \phi( x \cdot \theta^*)$, where the labels are generated by an arbitrary non-linear scalar link function $\phi$ applied to an unknown one-dimensional projection $\theta^*$ of the input data. By focusing on Gaussian data, several recent works have built a remarkable picture, where the so-called information exponent (related to the regularity of the link function) controls the required sample complexity. In essence, these tools exploit the stability and spherical symmetry of Gaussian distributions. In this work, building from the framework of [9], we explore extensions of this picture beyond the Gaussian setting, where both stability or symmetry might be violated. Focusing on the planted setting where $\phi$ is known, our main results establish that Stochastic Gradient Descent can efficiently recover the unknown direction $\theta^*$ in the high-dimensional regime, under assumptions that extend previous works  [41, 40].

§ INTRODUCTION

Over the past years, there has been sustained effort to enlarge our mathematical understanding of high-dimensional learning, particularly when using neural networks trained with gradient-descent methods — highlighting the interplay between algorithmic, statistical and approximation questions. An essential, distinctive aspect of such models is their ability to perform feature learning, or to extract useful low-dimensional features out of high-dimensional observations.

An appealing framework to rigorously analyze this behavior are sparse functions of the form $f(x) = \phi( \Theta_*^\top x)$, where the labels are generated by a generic non-linear, low-dimensional function $\phi: \mathbb{R}^k \to \mathbb{R}$ of linear features $\Theta^\top_* x$, with $\Theta_* \in \mathbb{R}^{d \times k}$ with $k \ll d$. While the statistical and approximation aspects of such function classes are by now well-understood [8, 6], the outstanding challenge remains computational, in particular in understanding the ability of gradient-descent methods to succeed. Even in the simplest setting of single-index models ($k=1$), and assuming that $\phi$ is known,  the success of gradient-based learning depends on an intricate interaction between the data distribution $x \sim \nu$ and the `link' function $\phi$; and in fact computational lower bounds are known for certain such choices [41, 35, 24, 16, 33].

Positive results thus require to make specific assumptions, either about the data, or about the link function, or both. On one end, there is a long literature, starting at least with [26, 32, 25], that exploits certain properties of $\phi$, such as invertibility or monotonicity, under generic data distributions satisfying mild anti-concentration properties [34, 21, 41, 40]. On the other end, by focusing on canonical high-dimensional measures such as the Gaussian distribution, the seminal works [9, 20] built a harmonic analysis framework of SGD, resulting in a fairly complete picture of the sample complexity required to learn generic link functions $\phi$, and revealing a rich asymptotic landscape beyond the proportional regime $n \asymp d$, characterized by the number of vanishing moments, or information exponent $s$ of $\phi$,
whereby $n \asymp d^{s-1}$ samples are needed for recovery. Since then, several authors have built and enriched this setting to multi-index models [1, 2, 14, 4], addressing the semi-parametric learning of the link function [12], as well as exploring SGD-variants [10, 7, 11, 13]. This harmonic analysis framework relies on two key properties of the Gaussian measure and their interplay with SGD: its spherical symmetry and its stability by linear projection. Together, they provide an optimization landscape that is well-behaved in the limit of infinite data, and enable SGD to escape the `mediocrity' of initialisation, where the initial direction $\theta_0$, in the high-dimensional setting, has vanishingly small correlation $| \theta_0 \cdot \theta^*| \simeq 1/\sqrt{d}$ with the planted direction $\theta^*$.

In this work, we study to what extent the `Gaussian picture' is robust to perturbations, focusing on the planted setting where $\phi$ is known. Our motivation comes from the fact that real data is rarely Gaussian, yet amenable to being approximately Gaussian via CLT-type arguments.
We establish novel positive results along two main directions: (i) when spherical symmetry is preserved but stability is lost, and (ii) when spherical symmetry is lost altogether. In the former, we show that spherical harmonics can be leveraged to provide a benign optimization landscape for SGD under mild regularity assumptions, for initialisations that can be reached with constant probability with the same sample complexity as in the Gaussian case.
In the latter, we quantify the lack of symmetry with robust projected Wasserstein distances, and show that for `quasi-symmetric' measures with small distance to the Gaussian reference, SGD efficiently succeeds for link functions with information exponent $s \leq 2$. Finally, using Stein's method, we address substantially `non-symmetric' distributions, demonstrating the strength and versatility of the harmonic analysis framework.

§ PRELIMINARIES AND PROBLEM SETUP

The focus of this work is to understand regression problems with input/output data $(x, y) \in \R^d \times \R $ generated by single-index models. This is a class of problems where the data labels are produced by a non-linear map of a one-dimensional projection of the input, that is

\begin{align}
\label{eq:y_model}
y = \phi(x \cdot \theta^*),
\end{align}

where $\phi:\mathbb{R} \to \mathbb{R}$ is also known as the link function, and $\theta^* \in \S$, the sphere of $\R^d$, is the hidden direction that the models wants to learn. Quite naturally, the learning is made through the family of generalized linear predictors $\mathcal{H} = \{\phi_{\theta} : x \to \phi( x \cdot \theta), \ \text{for }\, \theta\! \in \S \}$, built upon the link function (which is assumed known) and parametrized by the sphere.

Loss function.
We assume that the input data is distributed according to a probability $\nu \in \mathcal{P}(\R^d)$. Equation (<ref>) implies that the target function that produces the labels, $\phi_{\theta^*},$ lies in this parametric class. The overall loss classically corresponds to the average over all the data of the square penalisation $l(\theta, x) := (\phi_\theta(x) - \phi_{\theta^*}(x))^2$ so that the population loss writes

\begin{align}
\label{eq:mainsetup}
L(\theta) := \mathbb{E}_{\nu} \left[\big(\phi(x \cdot \theta) - \phi(x \cdot \theta^*)\big)^2\right] = \left\|\phi_{\theta} - \phi_{\theta^*}\right\|^2_{L^2_\nu},
\end{align}

where we used the notation $\|f\|^p_{L^p_\nu} = \E_\nu[|f|^p]$, valid for all $p \in \N^*$. Let us put emphasis on the fact that the loss $L$ is a non-convex function of the parameter $\theta$, hence it is not a priori guaranteed that gradient-based method are able to retrieve the ground-truth $\theta^*$. This often requires a precise analysis of the loss landscape, and where the high-dimensionality can play a role of paramount importance: we place ourselves in this high-dimensional setting for which the dimension is fixed but considered very large $d \gg 1$. Finally, we assume throughout the article that $\phi_\theta$ belongs to the weighted Sobolev space $W^{1,4}_\nu:=\{ \phi, \text{ such that }
\sup_{\theta \in \S} \left[\|\phi_\theta\|_{L^4_\nu} + \|\phi_\theta'\|_{L^4_\nu}\right] <  \infty \}$.

Stochastic gradient descent.
To recover the signal given by $\theta^* \in \S$, we run online stochastic gradient descent (SGD) on the sphere $\S$. This corresponds to having at each iteration $t \in \mathbb{N}^*$ a fresh sample $x_t$ drawn from $\nu$ and independent of the filtration $\mathcal{F}_t = \sigma(x_1, \dots, x_{t-1})$ and performing a spherical gradient step, with step-size $\delta>0$, with respect to $\theta \to l(\theta, x_t)$:

\begin{align}
\label{eq:SGD}
\theta_{t+1} &= \frac{\theta_t - \delta \nabla_\theta^\mathcal{S}  l(\theta_t, x_t)}{\left|\theta_t - \delta \nabla_\theta^\mathcal{S}  l(\theta_t, x_t)\right|}, \qquad \text{ with initialization } \theta_{0} \sim \mathrm{Unif}(\S),
\end{align}

An important scalar function that enables to track the progress of the SGD iterates is the correlation with the signal $m_\theta:=\theta \cdot \theta^* \in [-1,1]$. We will drop the subscript in case there is no ambiguity. Note that, due to the high-dimensionality of the setup, we have the following lemma:

For all $a > 0$, we have $ \P_{\theta_0}(m_{\theta_0} \geq a/\sqrt{d}) \leq a^{-1}e^{-a^2/4}  $. Additionally, for any $\delta>0$ such that $ \max\{a, \delta\} \leq \sqrt{d}/4$, we have the lower bound:
$\P_{\theta_0}(m_{\theta_0} \geq a/\sqrt{d}) \geq \frac{\delta}{4} e^ {- (a + \delta)^2} $.

This fact implies that, when running the algorithm in practice, it is initialized with high probability near the equator of $\S$, or at least in a band of typical size $1/\sqrt{d}$ (see Figure <ref> for a schematic illustration of this fact). Finally, we use the notation $\nabla_\theta^\mathcal{S}$ to denote the spherical gradient, that is $ \nabla_\theta^\mathcal{S} l(\theta, x) = \nabla_\theta  l(\theta, x) - (\nabla_\theta  l(\theta, x) \cdot \theta) \theta$. As $\nabla_\theta^\mathcal{S}  l(\cdot, x_t)$ is an unbiased estimate of $\nabla_\theta^\mathcal{S}  L$, it is expected that the latter gradient field rules how the SGD iterates travel across the loss landscape.

Loss landscape in the Gaussian case. As stressed in the introduction, this set-up has been studied by [20, 9] in the case where $\nu$ is the standard Gaussian, noted as $\gamma$ here to avoid any confusion for later. Let us comment a bit this case to understand what can be the typical landscape of this single-index problem.
Thanks to the spherical symmetry, the loss admits a scalar summary statistic, given precisely by the correlation $m_\theta$. Moreover, the loss admits an explicit
representation in terms of the Hermite decomposition of the link function $\phi$: if $\{h_j\}_j$ denotes the orthonormal basis of Hermite polynomials of $L^2_\gamma$, then $L(\theta) = 2\sum_j |\langle \phi, h_j\rangle|^2(1 - m^j) := \bar{\ell}(m)$. As a result,

the gradient field projected along the signal is a (locally simple) positive function of the correlation that behaves similarly to

\begin{align}
\label{eq:loss_property_gaussian}
- \nabla_\theta^\mathcal{S}  L(\theta) \cdot \theta^* \simeq C m^{s-1} (1 - m),
\end{align}

where $s \in \N^*$ is the index of the first non-zero of the Hermite coefficients $\{\langle \phi, h_j\rangle\}_j$.
This has at least three important consequences for the gradient flow: (i) if initialized positively, the correlation is an increasing function along the dynamics and there is no bad local minima in the loss landscape,  (ii) the parameter $s \in \N^*$ controls the flatness of the loss landscape near the origin and therefore controls the optimization speed of SGD in this region (iii) as soon as $m$ is large enough, the contractive term $1-m$ makes the dynamics converge exponentially fast.

Loss landscape in general cases. Obviously for general distributions $\nu$, the calculation presented in Eq.(<ref>) is no-longer valid. However, the crux of the present paper is that properties (i)-(ii)-(iii) are robust to the change of distribution and can be shown to be preserved under small adaptations. More precisely, we have the following definition.

We say that $L$ has the local polynomial growth of order $k \in \N^*$ and scale $b \geq 0$, if there exists $C > 0$ such that for all $m_\theta \geq b$,

\begin{equation}
- \nabla_\theta^\mathcal{S}  L(\theta) \cdot \theta^* \geq C (1 - m_\theta) \left(m_\theta - b\right)^{k-1}~.
\end{equation}

In such a case we say that $L$ satisfies $\LPG(k,b)$.

In this definition, and as showed later in specific examples given in Section <ref>, the scale parameter $b$ should be thought as a small parameter  proportional to $1/\sqrt{d}$.

If $\nu$ is Gaussian, we can rewrite Eq.(<ref>) and show that $L$ verifies $\LPG(s,0)$ for $s \in \N^*$, referred to as the information exponent of the problem in [9].
An important consequence of satisfying $\LPG(k,b)$ is that the the population landscape is free of bad local minima outside the equatorial band $\Sigma_b:= \{ \theta \in \S \ , \   m_\theta \leq b\}$. Therefore, when $b$ is of scale $1/\sqrt{d}$, Lemma <ref> indicates that one can efficiently produce initializations that avoid it.
Hence, this property is the fundamental ingredient that enables the description of the path taken by SGD. We derive it thoroughly in next Section. Section <ref> is devoted to showcasing generic examples when this property is satisfied.

§ STOCHASTIC GRADIENT DESCENT UNDER LPG

In this section, we derive the main results on the trajectory of the stochastic gradient descent. They state that property $\LPG(s,b/\sqrt{d})$ is in fact sufficient to recover the same quantitative guarantees as the one depicted in [9], despite the lack of Gaussianity of the distribution $\nu$. Recall that the recursion satisfied by the SGD iterates is given by Eq.(<ref>). To describe their movement, let us introduce the following notations: for all $t \in \mathbb{N}^*$, we denote the normalization by $r_t := \left|\theta_t - \delta \nabla_\theta^\mathcal{S}  \ell(\theta_t, x_t)\right|$ and the martingale induced by the stochastic gradient descent as $M_t := l(\theta_t, x_t) - \E_\nu [ l(\theta_t, x)]$.

Tracking the correlation. Recall that the relevant signature of the dynamics is the one-dimensional correlation $m_{t} = \theta_t \cdot \theta^*$. For infinitesimal step-sizes $\delta \to 0$, it is expected that $\theta_t$ follow the spherical gradient flow $\dot{\theta}_t = - \nabla_\theta^\mathcal{S} L(\theta_t)$, that translates naturally on the summary statistics $m_t$ as the following time evolution

\begin{align}
\label{eq:idealized_ODE}
\dot{m_t} = - \nabla_\theta^\mathcal{S} L(\theta_t) \cdot \theta^*.
\end{align}

The main idea behind the result of this section is to show that, even if the energy landscape near $m = m_0$ is rough at scale $1/\sqrt{d}$, the noise induced by SGD does not prevent $m$ to grow as the idealized dynamics described by the ODE (<ref>). Let us write the iterative recursion followed by $(m_t)_{t \geq 0}$: with the notation recalled above, for $t \in \mathbb{N}^*$, we have

\begin{align}
\label{eq:dynamics_main_m}
m_{t+1} =  \frac{1}{r_t}\left( m_t - \delta\nabla_\theta^S L(\theta_t) \cdot \theta^* - \delta\nabla_\theta^S M_t \cdot \theta^*  \right).
\end{align}

With this dynamics at hand, the proof consists in controlling both the discretization error through $r_t$ and the directional martingale induced by the term $\delta\nabla_\theta^S M_t \cdot \theta^*$.

*Weak recovery. As it is the case for the gradient flow, most of the time spent by the SGD dynamics is near the equator, or more precisely in a band of the type $\Sigma_{b,c} = \{ \theta \in \S, \ b/\sqrt{d} \leq m_\theta \leq c/\sqrt{d} \}$, where $b<c$ are constants independent of the dimension. Hence, the real first step of the dynamics is to go out any of these bands. This is the reason why it is natural to define $S_a := \{\theta \in \S, \, m_\theta \geq a\}$, the spherical cap of level $a \in (0,1)$ as well as the hitting time

\begin{align}
\tau^+_a := \inf\{ t \geq 0, \, m_{\theta_t} \geq a  \},
\end{align}

which corresponds to the first time $(\theta_t)_{t \geq 0}$ enters in $S_a$. We arbitrarily choose a numerical constant independent of dimension, say $a = 1/2$, and refer to the related hitting time $\tau^+_{1/2}$ as the weak recovery time of the algorithm.

Let $(\theta_t)_{t \geq 0}$ follow the SGD dynamics of Eq.(<ref>) and let $L$ satisfy $\LPG(s,b/\sqrt{d})$, with $b >0 $ and $s \in \N^*$, then, conditionally on the fact that $m_0 \geq 4b/\sqrt{d} $, for any $\varepsilon > 0$, we have

\begin{align}
\quad \tau^+_{1/2} \leq \left\{
\begin{array}{llll}
d \cdot C/\varepsilon & \text{ when $s = 1,\ \ \ $ and with the choice} & \delta = \varepsilon/d \\
d \log^2(d) \cdot C/\varepsilon & \text{ when $s = 2,\ \ \ $ and with the choice} & \delta = \varepsilon/(d \log d) \\
d^{s-1} \cdot C/\varepsilon & \text{ when $s \geq 3,\ \ \ $ and with the choice} & \delta = \varepsilon d^{-s/2}
\end{array}
\right. ~,
\end{align}

with probability at least $1 - C\varepsilon$, for some $C >0$ that depends solely on the link function $\phi$.

Let us comment on this result. It says that that the integer $s$ coming from the growth condition controls the hardness of exiting the equator of the sphere. Indeed, as can be seen in $\LPG(s,b/\sqrt{d})$, the larger the $s$, the smaller the gradient projection is and hence the less information the SGD dynamics has to move from the initialization. This result can be seen as an extension of <cit.> valid only in the Gaussian case ($b = 0$). Furthermore, the Gaussian case shows that Theorem <ref> is tight up to $\log(d)$ factors. Finally, note that the result is conditional to the fact that the initialization is larger that some constant factor of $1/\sqrt{d}$, which has at least constant probability to happen in virtue of Lemma <ref>. This probability can be lowered by any constant factor by sampling offline a constant factor (independent of $d$) of i.i.d. initializations and keeping the one that maximizes its correlation with $\theta^*$.

Finally, note that in all the cases covered by the analysis, it is possible to keep track of the constant $C$, and show that overall it depends only (i) on the property of $\nu$ w.r.t. the Gaussian on the one hand and (ii) on the Sobolev norm of the link function $\|\phi\|_{W_\nu^{1,4}}$ on the other hand.

*Strong recovery. We place ourselves after the weak recovery time described in the previous question and want to understand if the $(\theta_t)_{t \geq 0}$ dynamics goes to $\theta^*$ and if yes, how fast it does so. This is what we refer to as the strong recovery question, captured by the fact that the one-dimensional summary statistics  $m$ go towards $1$. Thanks to the Markovian property of the SGD dynamics, we have the equality between all time $s > 0$ marginal laws of

\begin{align*}
\left(\theta_{\tau^+_{\nicefrac{1}{2}} + s}\  \bigg| \  \tau^+_{1/2},\, \theta_{\tau^+_{1/2}}\right) \overset{\text{Law}}{=} \left(\theta_{s} \  \bigg| \  \theta_{s} = \theta_{\tau^+_{1/2}} \right),
\end{align*}

and hence the strong recovery question is equivalent to study the dynamics with initialization that has already weakly recovered the signal, i.e. such that $m_\theta  = 1/2$. We show that this part of the loss landscape is very different that the equator band in which the dynamics spends most of its times: in all the cases, we can choose stepsizes independent of the dimension and show that the time to reach the vicinity of $\theta^*$ will be independent of $d$.

Let $(\theta_t)_{t \geq 0}$ follow the SGD dynamics of Eq.(<ref>) and let $L$ satisfy $\LPG(s,b/\sqrt{d})$, with $b >0 $ and $s \in \N^*$, then, for any $\varepsilon > 0$, taking $\delta = \varepsilon \log(1/\varepsilon)^{-1}$,   we have that there exists a time $T > 0$, such that

\begin{align}
|1 - m_T| \leq \varepsilon, \text{ and }\quad |T - \tau^+_{1/2}| \leq C \log(1/\varepsilon)\varepsilon^{-1}
\end{align}

with probability at least $1 - C\varepsilon$, for some $C >0$ that depends solely of the link function $\phi$.

As introduced above, the important messages of this theorem are that (i) there is no difference between the different parameters setups captured by the information exponent $s$, and (ii) the time it takes to reach an $\varepsilon$-vicinity of $\theta^*$ is independent of the dimension. We remark that we decided to present Theorem <ref> resetting the step-size $\delta$ to put emphasis on the intrinsic difference between the two phases. Yet, we could have kept the same stepsize as in the weak recovery case: this obviously would slow down unnecessarily the second phase.

(-4,0) arc (180:360:4 and 1);
[dashed] (-4,0) arc (180:0:4 and 1);

(-3.95,.5) arc (180:360:3.95 and 1);
[dashed] (-3.95,.5) arc (180:0:3.95 and 1);

[thick] (-4,0) arc (180:0:4 and 4);

[teal, thick] plot [smooth, tension=0.5] coordinates
(2,-1.25,-1.75) (2.0653304520112862, -1.0432818908769334, -1.7570889677693344) (1.653304520112862, -1.35432818908769334, -1.7570889677693344) (1.53304520112862, -1.15432818908769334, -1.7570889677693344) (1.653304520112862, -1.15432818908769334, -1.7570889677693344) (1.653304520112862, -1.15432818908769334, -1.7570889677693344) (1.53304520112862, -1.15432818908769334, -1.7570889677693344) (1.353304520112862, -1.5432818908769334, -1.7570889677693344) (1.253304520112862, -1.15432818908769334, -1.7570889677693344)  (-0.8643467206473194, -1.2741589636098296, -1.713346571687929) (-0.5854241939363898, -1.5603375176592302, -1.7026887285749814) (-3.10167004625995113, -1.0207262123095036, -1.7403753036061693) (-4.0015938188224307037, -1.076955164501522, -1.797513647919381) (0.03958282309186134, -1.4126458731857201, -1.8082986090450741) (0.4304065824253258, -1.5714797424606084, -1.7826566209500503) (-0.18290568047143796, -1.3714851280644934, -1.665058121198839) (-0.2141140616704278, -1.262736795657368, -1.7735105695338182) (-0.7359662339324988, -1.1933133412385735, -1.7500513117700844) (-1.080017844455942, -1.4440976104037713, -1.7202310698906578) (-0.10790460267119584, -1.36709352843529, -1.7694202444428277) (-0.5717959946699571, -1.6709912933384345, -1.7287576853165842) (-1.0810227730873538, -1.2678003390731818, -1.7070830602634035) (-1.0534246646049377, -1.1545710992352378, -1.7360956651453512) (-0.7765052753289565, -1.289319675143908, -1.7077561390856055) (-3.3548928604060262, -1.1867863636986526, -1.781372913021331) (-4.29110656221205145, -0.6034746092966507, -1.7284470803877188) (-2.8352886387760792, -1.4516559311541479, -1.7389164213779167)
(-1.86715707903012986, -1.132604490560212, -1.7959886616293295)  (-2.36715707903012986, -0.132604490560212, -1.7959886616293295) (-2.36715707903012986, 1.132604490560212, -1.7959886616293295) (-1.96715707903012986, 2.132604490560212, -1.7959886616293295) (-0.96715707903012986, 3.132604490560212, -1.4959886616293295) (0, 4.0, 0)

(-1.86715707903012986, -1.42604490560212, -1.7959886616293295) node[anchor=south] ()  •;

(-1.46715707903012986, -0.802604490560212, -1.7959886616293295) node()  $\theta_{\mathsf{w}}$;

[thick, thick, ->] (0,0) – ++(0,4);
(0.25,4.35) node () $\theta^*$;

[thick, thick, ->] (0,0,0) – ++(2,-1.25,-1.75);
(2.2,-1.2,-2) node () $\theta_0$;

[black, thick, <->] (0,0) ++(172.5:4.25cm) arc (172.5:180:4.25cm) node [midway, left] () $b/\sqrt{d}\ $;

Sketch of the SGD dynamics. After a long time spent in a band of typical size $1/\sqrt{d}$, the dynamics escapes at weak recovery point $\theta_{\mathsf{w}}$ and then goes rapidly to the ground-truth $\theta^*$.

§ TYPICAL CASES OF LOSS LANDSCAPE WITH PROPERTY PROPERTY

In this section, we showcase two prototypical cases where the $\LPG$ holds true: the section <ref> deals with the spherically symmetric setting, whereas the section <ref> describes a perturbative regime where the distribution is approximately Gaussian is a quantitative sense.

§.§ The symmetric case

We start our analysis with the spherically symmetric setting. We show that a spherical harmonic decomposition provides a valid extension of the Hermite decomposition in the Gaussian case, leading to essentially the same quantitative performance up to constant (in dimension) factors.

*Spherical Harmonic Representation of the Population Loss.
We express the data distribution $\nu$ as a mixture of uniform measures $\nu = \int_0^\infty \tau_{r,d}\, \rho(dr) $, where $\tau_{r,d}=\mathrm{Unif}(r \S)$.
Let $\tau_d = \tau_{1,d}$ and $u_d\in \mathcal{P}([-1,1])$ be the projection of $\tau_d$ onto one direction, with density given in close form by $u_d(dt) = Z^{-1}(1-t^2)^{(d-3)/2} \mathbf{1}(|t|\leq 1) dt $, where $Z$ is a normarlizing factor.

Let $\{P_{j,d}\}_{j \in \mathbb{N}}$ be the orthogonal basis of Gegenbauer polynomials  of $L^2_{u_d}([-1,1])$, normalized such that $P_{j,d}(1) = 1$ for all $j, d$. For each $r>0$, consider

$$l_r(\theta) := \langle\phi_\theta, \phi_{\theta^*}  \rangle_{\tau_{r,d}}  =   \langle \phi^{(r)}_\theta, \phi^{(r)}_{\theta^*}\rangle_{\tau_{d}}~,$$

where we define $\phi^{(r)}: [-1,1] \to \mathbb{R}$ such that $\phi^{(r)}(t):= \phi( r t)$. We write its decomposition in $L^2_{u_d}([-1,1])$ as
$\phi^{(r)} = \sum_{j} \alpha_{j,r,d} P_{j,d}~,~\text{with }~\alpha_{j,r,d} = \frac{\langle \phi^{(r)}, P_{j,d} \rangle}{\| P_{j,d} \|^2}~.$

Let $\Omega_{d}$ be the Lebesgue measure of $\S$, and $N(j,d) = \frac{2d+j-2}{d}\binom{d+j-3}{d-1}$ the so-called dimension on the spherical harmonics of degree $j$ in dimension $d$.
From the Hecke-Funk representation formula <cit.> and the chosen normalization $P_{j, d} (1) = 1$, we have $\| P_{j,d} \|_2 = \left(\frac{\Omega_{d-1}}{\Omega_{d-2} N(j,d)} \right)^{1/2}$ <cit.> and obtain finally

\begin{align}
l_r(\theta) &=
%\langle\phi^{(r)}_\theta, \phi^{(r)}_{\theta^*}\rangle_{\tau_{d}} = \frac{\Omega_{d-2}}{\Omega_{d-1}}\sum_j \frac{\Omega_{d-1}}{\Omega_{d-2}} N(j,d)^{-1} \alpha_{j,r,d}^2 P_{j,d}( \theta \cdot \theta^*) \\
\sum\nolimits_{j} \bar{\alpha}_{j,r,d}^2 P_{j,d}( \theta \cdot \theta^*)~,
\end{align}

where we defined for convenience $\bar{\alpha}_{j,r,d} = \alpha_{j,r,d} / \sqrt{N(j,d)}$. As a result, it follows that the overall loss writes as solely the correlaton $m_\theta = \theta \cdot \theta^*$ as

\begin{align}
\label{eq:lossmain}
%&= \mathbb{E}_\nu [ |\phi_\theta - \phi_{\theta^*}|^2] \nonumber
%&= C_\phi - 2 \langle \phi_\theta, \phi_{\theta^*} \rangle_{\nu} \nonumber = C_\phi - 2\int_0^\infty \langle \phi_\theta, \phi_{\theta^*} \rangle_{\tau_{r,d}} \rho(dr) \nonumber  \\
%&= 2 \mathbb{E}_\eta [\phi^2] - 2 \sum_{j} \left(\int_0^\infty \bar{\alpha}_{j,r,d}^2 \rho(dr) \right) P_{j,d}(\theta \cdot \theta^*)
= 2\|\phi\|^2_{L^2_\nu} -2 \sum\nolimits_{j} {\beta}_{j,d} P_{j,d}(m_\theta)~
:= \ell(m) ~,
\end{align}

where ${\beta}_{j,d} = \int_0^\infty \bar{\alpha}_{j,r,d}^2 \rho(dr) \geq 0$.

Unsurprisingly, we observe that, thanks to the spherical symmetry and analogous to the Gaussian case, the loss still admits $m=\theta \cdot \theta^*$ as a summary statistics. Yet, it is represented in terms of Gegenbauer polynomials, rather than monomials as in the Gaussian case. The monomial representation is a consequence of the stability of the Gaussian measure, as seen by the fact that
$\langle \phi_\theta, \phi_{\tilde{\theta}} \rangle_{\gamma_d} = \langle \phi, \mathsf{A}_{\theta \cdot \tilde{\theta}} \phi \rangle_{\gamma}~,$
where $(\mathsf{A}_m f)(t) = \mathbb{E}_{z \sim \gamma} [f ( mt + \sqrt{1-m^2}z)]$ has a semi-group structure (it is even known in fact as the Ornstein-Ulhenbeck semi-group).

Let $\eta \in \mathcal{P}(\mathbb{R})$ be the marginal of $\nu$ along any direction. The following proposition gives a closed form formula of the coefficients $\beta_{j,d}$, represented as integrals over the radial distribution $\rho$ and projections of the link function $\phi$:

The $\beta_{j,d}$ defined in (<ref>) have the integral representation

\begin{equation}
\beta_{j,d} = \langle \phi, \mathcal{K}_j \phi \rangle_{L^2_\eta} ~,
\end{equation}

where $\mathcal{K}_j$ is a positive semi-definite integral operator of $L^2_\eta$ that depends solely on $\rho$ and $\phi$.

Note that a closed form expression of $\mathcal{K}_j$ can be found in Appendix <ref>. The above proposition is in fact the stepping stone to calculate properly the information exponent that plays a crucial role in the property $\LPG$. This is given through the link between the spectrum $\beta_{j,d}$ and the decomposition of $\phi$ in the $L^2_\eta$ orthogonal basis of polynomials, that we denote by $\{ q_j\}_j$.

Let $s = \inf\{ j; \beta_{j,d} > 0\}$ and $\tilde{s} = \inf\{ j; \langle \phi, q_j\rangle_{\eta} \neq 0\}$. Then $s \leq \tilde{s}$.

Thus, the number of vanishing moments of $\phi$ with respect to the data marginal $\eta$ provides an upper bound on the `effective' information exponent of the problem $s$, as we will see next.

Local Polynomial Growth.
From (<ref>), and as
$\nabla_\theta^{\mathcal{S}} L(\theta) = \ell'(m) (\theta^* - m \theta)~,$ we directly obtain

\begin{equation}
-\nabla_\theta^{\mathcal{S}} L(\theta) \cdot \theta^* = -(1 - m^2) \ell'(m) =  2 (1 - m^2) \sum\nolimits_{j} \beta_{j,d} P'_{j,d}(m)~,
\end{equation}

which is the quantity we want to understand to exhibit the property $\LPG$ in this case. Hence, we now turn into the question of obtaining sufficient guarantees on the coefficients $(\beta_{j,d})_j$ that ensure local polynomial growth.
Since the typical scale of initialization for $m$ is $\Theta(1/\sqrt{d})$,
our goal is to characterize sufficient conditions of local polynomial growth with $b = O(1/\sqrt{d})$.

For that purpose, let us define two key quantities of Gegenbauer polynomials:

\begin{align}
\upsilon_{j,d} & := -\min_{t \in (0,1)} P_{j,d}(t)~, ~~&\text{(smallest value)}\\
z_{j,d} & := \arg\max\left\{ t \in (0,1); P_{j,d}(t)=0 \right\}~.~~&\text{(largest root)}
\end{align}

We have the following sufficient condition based on the spectrum $(\beta_{j,d})_j$:
[Spectral characterization of $\LPG$]propositionsufficientcond

Suppose there exist constants $K,C>0$ and $s \in \mathbb{N}$ such that we both have $\beta_{s,d} \geq C$ and $\sum_{j > s} \beta_{j,d} j (j+d-2)  \upsilon_{j-1,d+2} \leq K d^{(3-s)/2}$ .

Then, taking $s^*$ as the infimum of such $s$, $L$ has the property $\LPG(s^*-1,z_{s^*,d})$. In particular, whenever $s^*\ll d$, we have $z_{s^*,d} \leq 2\sqrt{s^*/d}$.

This proposition thus establishes that, modulo a mild regularity assumption expressed though the decay of the coefficients $\beta_{j,d}$, the spherically symmetric non-Gaussian setting has the same geometry as the Gaussian setting, for correlations slightly above the equator. Crucially, the required amount of correlation to `feel' the local polynomial growth is a $O(\sqrt{s})$ factor from the typical initialization, and can be thus obtained with probability $\simeq e^{-s}$ over the initialization, according to Lemma <ref>, a lower bound which is independent of $d$.

The sufficient condition for $\beta_{j,d}$ appearing in Proposition <ref> involves the minimum values $\upsilon_{j,d}$ of Gegenbauer polynomials $P_{j,d}$, as well as sums of the form $\sum_j j^2 \beta_{j,d}$.
In order to obtain a more user-friendly condition, we now provide an explicit control of $\upsilon_{j,d}$, and leverage mild regularity of $\phi$ to control $\sum_j j^2 \beta_{j,d}$. This motivates the following assumption on $\phi$ and $\nu$:

The link function $\phi$ satisfies $\phi \in L^2_\eta$ and $\phi' \in L^4_\eta$, and the radial distribution $\rho$ has finite fourth moment $\mathbb{E}_\rho [r^4] < \infty$.

Assume that $\phi$ and $\nu$ satisfy Assumption <ref>, and let $s^* = \inf\{j; \langle \phi, q_j \rangle_\eta \neq 0 \}$. Then $L$ has the property $\LPG(s^*-1,2\sqrt{s^*/d})$.

The proof is provided in Appendix <ref>. At a technical level, the main challenge in proving Theorem <ref> is to achieve a uniform control of $\upsilon_{j,d}$ in $j$, a result which may be of independent interest. We address it by combining state-of-the-art bounds on the roots of the Gegenbauer polynomials, allowing us to cover the regime where $j$ is small or comparable to $d$, together with integral representations via the Cauchy integral formula, providing control in the regime of large $j$. On the other hand, we relate the sum $\sum_j j^2 \beta_{j,d}$ to a norm of $\phi'$ using a Cauchy-Schwartz argument, where we leverage the fourth moments from Assumption <ref>.

Since we are in a setting where $\phi$ is known, an alternative to the original recovery problem from Eq (<ref>) is to consider a pure Gegenbauer `student' link function of the form $\tilde{\phi} = P_{s,d}$, where $s$ is the information exponent from Proposition <ref>. Indeed, the resulting population loss
$\tilde{L}(\theta) = \mathbb{E} [ (\tilde{\phi}(x \cdot \theta) - \tilde{\phi}( x \cdot \theta^*) )^2]$ satisfies the $\LPG$ property, as easily shown in Fact <ref>.

For the sake of completeness, we describe more precisely two concrete case studies below.

[Uniform Measure on the Sphere]
When $\nu = \text{Unif}(\sqrt{d} \S)$, we have $\rho = \delta_{\sqrt{d}}$, and therefore $\beta_{j,d} = \bar{\alpha}^2_{j, \sqrt{d},d}$.
In that case, the orthogonal polynomial basis $\{q_j(t)\}_j$ of $L^2_\eta$ coincides with the rescaled Gegenbauer polynomials, $q_j(t) = P_{j,d}(t/\sqrt{d}) $.
Consider now a link function
$\phi$ with $s-1$ vanishing moments with respect to $L^2_\eta$, i.e. such that
$\bar{\alpha}_{j,d}=\langle \phi, q_j \rangle_\eta = 0$ for $j < s$ and $\bar{\alpha}_{s,d}=\langle \phi, q_s \rangle_\eta \neq 0$; and with sufficient decay in the higher harmonics as to satisfy the bound on the sum presented in Proposition <ref> (for example, $\phi(t) = q_s(t)$ trivially satisfies this condition).
Then Proposition <ref> applies and we conclude that the resulting population landscape satisfies $\LPG(s-1, O(\sqrt{s/d}))$.

In [41, 40] it is shown that monotonically increasing link functions[or link functions where their monotonic behavior dominates; see [40].] lead to a benign population landscape, provided the data distribution $\nu$ satisfies mild anti-concentration properties.
We verify that in our framework.

[Non-decreasing $\phi$]
Indeed, Proposition <ref> is verified with $s=1$, provided $\phi' \in L^2_\eta$. Indeed, if $\phi\neq 0$ is monotonic, then we have

$\beta_1 = \langle \phi, \mathcal{K}_1 \phi \rangle_\eta = C_d \left(\mathbb{E}_{\eta} [t \phi(t)] \right)^2 \neq 0~,$
since we can assume without loss of generality that $\phi(t)\geq 0$ for $t\geq 0$ and $\phi(t) \leq 0$ for $t \leq 0$.

We emphasize that the results of [41, 40] extend beyond the spherically symmetric setting, which is precisely the focus of next section.

§.§ Non-Spherically Symmetric Case

We now turn to the setting where $\nu$ is no longer assumed to have spherical symmetry.
By making further regularity assumptions on $\phi$, our main insight is that distributions that are approximately symmetric (defined in an appropriate sense) still benefit from a well-behaved optimization landscape.

*Two-dimensional Wasserstein Distance.
When $\nu$ is not spherically symmetric, the machinery of spherical harmonics does not apply, and we thus need to rely on another structural property. Consider a centered and  isometric data distribution $\nu \in \mathcal{P}_2(\mathbb{R}^d)$, i.e. such that $\mathbb{E}_\nu x=0$ and $\Sigma_\nu=\mathbb{E}_\nu [x x^\top] = I_d$.
We consider the two-dimensional $1$-Wassertein distance <cit.> –see also [30]– between a pair of distributions $\nu_a, \nu_b \in \mathbb{P}(\mathbb{R}^d)$, defined as

\begin{equation}
\widetilde{W}_{1,2}(\nu_a, \nu_b) := \sup_{P \in \mathrm{Gr}(2, d)} W_1( P_\# \nu_a, P_\# \nu_b)~,
\end{equation}

where the supremum runs for any two-dimensional subspace $P \in \mathrm{Gr}(2, d)$, and $P_\# \nu \in \mathcal{P}(\mathbb{R}^2)$ is the projection (or marginal) of $\nu$ onto the span of $P$.
$\widetilde{W}_{1,2}$ is a distance (<cit.>) and measures the largest $1$-Wasserstein distance between any two-dimensional marginals.

We are in particular interested in the setting where $\nu_a = \nu$ is our data distribution, and $\nu_b$ is a reference symmetric measure –  for instance the standard Gaussian measure $\gamma_d$. Consider the fluctuations

\begin{align}
\Delta_L(\theta) & := |L(\theta) - \bar{\ell}(m_\theta)| \text{ and }~,\\
\Delta_{\nabla L}(\theta) &:= | \nabla^{\mathcal{S}}_\theta L(\theta) \cdot \theta^* - \bar{\ell}'(m_\theta) (1 - m_\theta^2) |~,
\end{align}

where $\bar{\ell}(m)$ is the Gaussian loss defined in Section <ref> and $L(\theta) = \mathbb{E}_\nu [ | \phi(  x \cdot \theta ) - \phi(  x \cdot \theta^* )|^2]$.

$\Delta_L$ and $\Delta_{\nabla L}$ thus measure respectively the fluctuations of the population loss and the relevant (spherical) gradient direction.
By making additional mild regularity assumptions on the link function $\phi$, we can obtain a uniform control of the population loss geometry using the dual representation of the $1$-Wasserstein distance.

[Regularity of link function]assumptionasslipreg

We assume that $\phi, \phi'$ are both $B$-Lipschitz,
and that $\phi''(t) = O(1/t)$.

The data distribution $\nu$ is $M$-subgaussian: for any $v \in \S$, we have $\|  x\cdot v  \|_{\psi_2} \leq M $, where $\|z\|_{\psi_2}:= \inf\{ t>0; \, \mathbb{E}[\exp(z^2/t^2)\leq 2\}$ is the Orlitz-2 norm.

[Uniform gradient approximation]propositiongradwassers

Under Assumptions <ref> and <ref>, for all $\theta \in \S$,

\begin{equation}
%\sup_{\theta \in \S} | \langle \nabla^{\mathbb{S}}_\theta L(\theta), \theta^* \rangle - \bar{\ell}'(m) (1 - m^2) | \leq 2(L'B + L^2) \widetilde{W}_{1,2}(\nu, \gamma)~,
%    \sup_{\theta \in \S} | \langle \nabla^{\mathbb{S}}_\theta L(\theta), \theta^* \rangle - \bar{\ell}'(m) (1 - m^2) | \leq \tilde{C} \widetilde{W}_{1,2}(\nu, \gamma) \left(\log(\widetilde{W}_{1,2}(\nu, \gamma)^{-1})\right)^2~,
%      | \langle \nabla^{\mathbb{S}}_\theta L(\theta), \theta^* \rangle - \bar{\ell}'(m) (1 - m^2) | \leq \sqrt{1-m^2} \tilde{C} \widetilde{W}_{1,2}(\nu, \gamma) \left(\log(\widetilde{W}_{1,2}(\nu, \gamma)^{-1})\right)^2~,
\Delta_{\nabla L}(\theta) = (1-m^2) O\left( \widetilde{W}_{1,2}(\nu, \gamma) \log ( \widetilde{W}_{1,2}(\nu, \gamma)^{-1} ) \right)
%(1-m^2)^{3/2} \widetilde{W}_{1,2}(\nu, \gamma) \left|\log\left(\left((1-m^2) \widetilde{W}_{1,2}(\nu, \gamma)\right)^{-1/C} \right)\right|~,
\end{equation}

where the $O(\cdot)$ notation only hides constants appearing in Assumptions <ref> and <ref>.

In words, the population gradient under $\nu$ is viewed as a perturbation of the population
gradient under $\gamma$, which has the well-behaved geometry already described in Section <ref>. These perturbations
can be uniformly controlled by the projected 1-Wasserstein distance, thanks to the subgaussian tails of $\nu$.

Our focus will be in situations where $\widetilde{W}_{1,2}(\nu, \gamma)=O(1/\sqrt{d})$. This happens to be the `natural' optimistic scale for this metric in the class of isotropic distributions $\Sigma_\nu = I_d$, as can be seen for instance when $\nu = \mathrm{Unif}(\sqrt{d} \S)$. Under such conditions, it turns out that link functions with information exponent $s\leq 2$ can be recovered with simple gradient-based methods, by paying an additional polynomial (in $d$) cost in time complexity.

The Gaussian information exponent of $\phi$, $s:=\arg\min\{ j; \langle \phi, H_j \rangle  \neq 0\}$ satisfies $s \leq 2$.

The projected Wasserstein distance satisfies
$\widetilde{W}_{1,2}(\nu, \gamma) \leq M'/\sqrt{d}$.

[$\LPG$, non-symmetric setting]propositionpropgrad_expo2

Under Assumptions <ref>, <ref>, <ref> and <ref>, $L$ verifies $\LPG\left(1,O\left(\sqrt{\frac{(\log d^\kappa)}{d}}\right)\right)$, where $\kappa$ depends only on $B,M,M'$.

This proposition illustrates the cost of breaking spherical symmetry in two aspects: (i) it requires additional regularity on $\phi$, and notably restricts its (Gaussian) information exponent to $s=2$, and (ii) the scale to reach $\LPG$ is now no longer dimension-free, but has a polynomial dependency on dimension, since from Lemma <ref>, picking $\delta$ any positive constant we have

$$\P_{\theta_0}\left(m_{\theta_0} \geq \sqrt{\log d^\kappa}/\sqrt{d}\right) \geq \frac{\delta}{4}  e^ {- (\sqrt{\log d^\kappa} + \delta)^2} \geq \Omega\left( d^{-(1+o(1))\kappa}\right)~. $$

At present, we are not able to rule this out as a limitation of our proof;  establishing whether this polynomial dependency on dimension is an inherent cost of the symmetry breaking is an interesting question for future work.

While assumptions <ref>, <ref> and <ref> are transparent and impose only mild conditions on the link function and tails of $\nu$, the `real' assumption of Proposition <ref> is the concentration of $\widetilde{W}_{1,2}(\nu, \gamma)$ (Assumption <ref>). The ball $\{ \nu; \widetilde{W}_{1,2}(\nu, \gamma) = O(1/\sqrt{d}) \}$ contains many non-symmetric measures, for instance empirical measures sampled from $\gamma$ with $n = \omega(d^2)$ <cit.>, and we suspect it contains many other examples, such as convolutions of the form $\nu \ast \gamma_{\sigma}$ arising for instance in diffusion models. That said, one should not expect the distance $\widetilde{W}_{1,2}(\nu, \gamma)$ to be of order $1/\sqrt{d}$ for generic `nice' distributions $\nu$; for instance, log-concave distributions are expected to satisfy $W_1( P_\# \nu, \gamma) \simeq 1/\sqrt{d}$ for most subspaces $P$, as captured in Klartag's CLT for convex bodies [27]. In summary, many situations of interest fall outside this regime, which motivates us to relax the uniform Wasserstein criterion.

Localized Growth via Stein's method.
To illustrate the mileage of the previous techniques beyond this `quasi-symmetric' case, we consider now an idealised setting where the data is drawn from a product measure $\nu = \eta^{\otimes d}$, with $\eta \in \mathcal{P}(\mathbb{R})$ and $W_1( \eta, \gamma_1) = \Theta(1)$. In other words, $x=(x_1, \ldots, x_d) \sim \nu$ if $x_i \sim \eta$ are i.i.d. In this setting, the distances $W_1( P_\# \nu, \gamma)$ reflect a CLT phenomena, which requires the subspace $P$ to `mix' across independent variables. Consequently, one may expect the expression of the hidden direction $\theta^*$ in the canonical basis to play a certain role.
For that purpose, we make the following additional regularity assumption on the tails of $\phi$ to simplify the quantitative bounds:
[Additional Regularity in third derivatives] assumptionassthird

$\phi$ admits four derivatives bounded by $L$, with $|\phi^{(3)}(t)| = O(1/t)$ and $|\phi^{(4)}(t)| = O(1/t^2)$. Moreover, the third moment of the data distribution is finite:
$\tau_3=\mathbb{E}_{t \sim \eta}[t^3] < \infty$.

Stein's method provides a powerful control on $\Delta_L(\theta)$ and
$\Delta_{\nabla L}(\theta)$, as shown by the following result:

[Stein's method for product measure]propositionpropstein

Let $\chi(\theta, \theta^*):= \| \theta\|_4^2 + \| \theta^*\|_4^2$.
Under Assumptions <ref>, <ref> and <ref>, there exists a universal constant $C=C(M, B, \tau_3)$ such that

\begin{equation}
\label{eq:loss_conc}
\Delta_L(\theta) \leq C \chi(\theta, \theta^*)~,\text{and }~\Delta_{\nabla L}(\theta) \leq C \sqrt{1-m^2} \chi(\theta, \theta^*)~.
%        \max(\Delta_{L}(\theta), \Delta_{\nabla L}(\theta)) \leq C (1-m_\theta^2) O\left(\chi(\theta, \theta^*) \log( \chi(\theta, \theta^*)^{-1} ) \right)~.
%   \left(\log \chi(\theta, \theta^*)(1-m_\theta^2) \right)^2~.
%        \left( \|\theta\|_3^3 + \| \theta\|_4^2 + \|\theta^*\|_3^3 + \| \theta^*\|_4^2 \right) \left(\log\left( \|\theta\|_3^3 + \| \theta\|_4^2 + \|\theta^*\|_3^3 + \| \theta^*\|_4^2 \right)\right)^2~.%\left( \|\theta\|_3^3 + \| \theta\|_4^2 + \|\theta^*\|_3^3 + \| \theta^*\|_4^2 \right) ~,
\end{equation}

The proof is based on the Stein's method for multivariate variables <cit.> with independent entries, which provides a quantitative CLT bound. Contrary to the quasi-symmetric case, here the concentration is not uniform over the sphere, but crucially depends on the sparsity of both $\theta$ and $\theta^*$, measured via the $\ell_4$ norms $\| \theta \|_4$, $\| \theta^* \|_4$: for incoherent, non-sparse directions, we have $\|\theta \|_4^2 \simeq 1/\sqrt{d}$, recovering the concentration rate that led to Proposition <ref>,
while for sparse directions we have $\| \theta\|_4^2 = \Theta(1)$, indicating an absence of concentration to the Gaussian landscape.

Therefore, the natural conclusion is to assume a planted model where $\theta^*$ is incoherent with the data distribution, i.e. $\| \theta^* \|_4 = O(d^{1/4})$. While the $\LPG$ property does not directly apply in this setting, we outline an argument that suggests that the single-index model can still be efficiently solved using gradient-based methods.
For that purpose, we assume that $\theta^*$ is drawn uniformly in $\S$, which implies that its squared-$L^4$ norm $\| \theta^*\|_4^2$ is of order $d^{-1/2}$ with high probability:

Assume $\theta^* \sim \mathrm{Unif}(\S)$. Then $\mathbb{P}( \| \theta^*\|_4^2 \leq C /\sqrt{d}) \geq 1 - C'\exp(-C)$.

Because $\theta_0$ is also drawn uniformly on the sphere, the typical value of $\chi(\theta, \theta^*)$ is of order $d^{-1/2}$. For Gaussian information exponent $s=2$, the population gradient under Gaussian data satisfies $-\nabla_\theta^{\S}L_\gamma(\theta) \cdot \theta^* \geq C m_\theta$. As a consequence, by Proposition <ref>, whenever $|\theta_0 \cdot \theta^*| > c \sqrt{d}$ (which happens with constant probability lower bounded by $e^{-c}$), we enter a `local' $\LPG$ region where $-\nabla_\theta^{\S}L(\theta) \cdot \theta^* \geq C (m - c/\sqrt{d}) > 0$.
While this condition is sufficient in the quasi-symmetric setting to start accumulating correlation (Theorem <ref>), now this event is conditional on $\theta$ being dense, ie so that $\chi(\theta, \theta^*) = O(1/\sqrt{d})$.

Since the typical value of $\chi(\theta, \theta^*)$ is of scale $1/\sqrt{d}$, one would expect that SGD will rarely visit sparse points where $\chi(\theta, \theta^*)\gg O(1/\sqrt{d})$, and thus that the local $\LPG$ property will be valid for most times during the entropic phase of weak recovery —
and therefore that the correlation $m_\theta$ will pile-up as in the quasi-symmetric setting.

We summarise this property in the following conjecture:

Assume $\theta^*, \theta_0$ are drawn from the uniform measure, and let $\theta_t$ be the $t$-th iterate of SGD with $\delta \simeq 1/(d \log d)$.
There exists a universal constant $C$ such that for any $\xi>0$, we have

\begin{equation}
\mathbb{P}\left( \sup_{t \leq T} \| \theta_t\|_4^2 \geq \sqrt{\frac{{\xi} \log T }{{d}}} \right) \leq C \exp(-\xi^2 d)~.
\end{equation}

Since the time to escape mediocrity in the case $s=2$ is $T \simeq d \log(d)^2$, this conjecture would imply that SGD does not effectively `see' any sparse points, and thus escapes mediocrity. If one assumed that in this phase the dynamics is purely noisy, now pretending that $\theta_i$ were drawn independently from the uniform measure, and that $\|\theta\|_4^4$ is approximately Gaussian with mean $d^{-1}$ and variance $d^{-3}$, the result follows by simple concentration. The challenging aspect of Conjecture <ref> is precisely to handle the dependencies across iterates, as well as the spherical projection steps.

§ EXPERIMENTS

In order to validate our theory, and inspect the degree to which our bounds may be pessimistic, we consider empirical evaluation of the training process in our two primary settings.  Specifically, we consider random initialization on the half-sphere (with the sign chosen to induce positive correlation as in [5]), and investigate how often strong recovery occurs relative to the information exponent of the link function.

Correlation with the true signal throughout training, under different choices of link function and input distribution.

The 4th degree Gegenbauer polynomial with dimension 50, 75, 100.  Each is equivalent (up to rescaling) to the loss landscape of learning the 4th degree Gegenbauer in the appropriate dimension.  Points indicate largest zeros.

*Symmetric Case

For the spherically symmetric setting, we experiment with the input distribution that is uniform on the sphere.  We are primarily interested in verifying that, unlike the Gaussian case, strong recovery depends on whether the initial correlation is sufficiently high to avoid local minima and benefit from the $\LPG$ guarantee.  This is not evident in the 2nd degree Gegenbauer case, which is monotonic and quickly reaches strong recovery, but it is clear from the 4th degree Gegenbauer link function.

In the infinite sample setting, Figure <ref> exactly characterizes the loss landscape when learning the 4th degree Gegenbauer under inputs uniform on $\S$ for different values of $d$.  Note that the largest zero for $d = 50$ occurs at $\approx \pm 0.31$, and the loss is monotonic for $m$ values initialized outside that region.  This phenomenon persists for higher dimensions, and one may observe that $d$ increases, the critical points become smaller in magnitude, according to the scaling $\simeq \sqrt{1/d}$.

The bottom right subplot in Figure <ref> indicates training runs in this setting, where red lines are initialized uniformly on the sphere, and blue lines are initialized uniformly conditioned on $m = 0.4$, which is slightly past the last zero of the polynomial.  We observe that random initialization infrequently exceeds the threshold necessary for strong recovery, but planting the initialization above this threshold gives a high probability of recovery.

*Non-Symmetric Case

For the non-spherically symmetric setting, we compare the performance of Gaussian inputs with inputs that are approximately Gaussian under a two-dimensional projection.  For simplicity, we loosen our assumptions slightly, and consider the input distribution as the $d$ dimensional product distribution of uniform random variables (rescaled to have unit variance), and allow for a non-Lipschitz link function.  Here, we are primarily interested in whether Assumption <ref> is tight and $s \leq 2$ is necessary for recovery, as well as whether Conjecture <ref> holds in practice.

To evaluate, we compare strong recovery rates when training on a "tricky" function with $s = 2$ (chosen to be $\frac{1}{2} \left(h_2 - h_3 - h_4 + h_5 \right)$) versus a function with $s = 3$ (simply the degree three hermite polynomial $h_3$).  We make this choice for the $s = 2$ function in order to produce a function which is not monotonic, for which learning is easy under many distributions, as discussed in [41].

In Figure <ref> we observe that strong recovery reliably occurs for both the Gaussian and hypercube input distributions when $s = 2$.  There is more variance in the Gaussian runs, likely because the magnitude of the gradients will be larger due to the inclusion of high degree terms.  But for the $s = 3$ case, the Gaussian distribution converges quickly while the hypercube distribution frequently cannot escape the equator.

§ CONCLUSION AND PERSPECTIVES

In this work, we have asked whether the remarkable properties of high-dimensional Gaussian SGD regression of single-index models are preserved as one loses some key aspects that make Gaussian distributions so special (and so appealing for theorists).
Our results are mostly positive, indicating a robustness of the Gaussian theory, especially within the class of spherically symmetric distributions, where a rich spherical harmonic structure is still available.  As one loses spherical symmetry, the situation becomes more dire, motivating a perturbative analysis that we have shown is effective via projected Wasserstein and Stein couplings.

That said, there are several open and relevant avenues that our work has barely touched upon, such as  understanding whether the robustness can be transferred to other algorithms beyond SGD, or addressing the semi-parametric problem when the link function is unknown, along the lines of [12, 2, 14, 11]. A particularly interesting direction of future work is to extend the analysis of product measures to `weakly dependent' distributions, motivated by natural images where locality in pixels captures most (but not all) of the statistical dependencies. Stein's method appears to be a powerful framework that can accommodate such weak dependencies, and deserves future investigation.

[1]
Abbe, E., Boix-Adsera, E., and Misiakiewicz, T. (2022).
The merged-staircase property: a necessary and nearly sufficient
condition for sgd learning of sparse functions on two-layer neural networks.
arXiv preprint arXiv:2202.08658.

[2]
Abbe, E., Boix-Adsera, E., and Misiakiewicz, T. (2023).
Sgd learning on neural networks: leap complexity and saddle-to-saddle
arXiv preprint arXiv:2302.11055.

[3]
Area, I., Dimitrov, D. K., Godoy, E., and Ronveaux, A. (2004).
Zeros of gegenbauer and hermite polynomials and connection
Math. Comput., 73:1937–1951.

[4]
Arnaboldi, L., Stephan, L., Krzakala, F., and Loureiro, B. (2023).
From high-dimensional & mean-field dynamics to dimensionless odes: A
unifying approach to sgd in two-layers networks.
arXiv preprint arXiv:2302.05882.

[5]
Arous, G. B., Gheissari, R., and Jagannath, A. (2021).
Online stochastic gradient descent on non-convex losses from
high-dimensional inference.
The Journal of Machine Learning Research, 22(1):4788–4838.

[6]
Bach, F. (2017).
Breaking the curse of dimensionality with convex neural networks.
The Journal of Machine Learning Research, 18(1):629–681.

[7]
Barak, B., Edelman, B., Goel, S., Kakade, S., Malach, E., and Zhang, C. (2022).
Hidden progress in deep learning: Sgd learns parities near the
computational limit.
Advances in Neural Information Processing Systems,

[8]
Barron, A. R. (1993).
Universal approximation bounds for superpositions of a sigmoidal
IEEE Transactions on Information theory, 39(3):930–945.

[9]
Ben Arous, G., Gheissari, R., and Jagannath, A. (2021).
Online stochastic gradient descent on non-convex losses from
high-dimensional inference.
Journal of Machine Learning Research (JMLR), 22:106–1.

[10]
Ben Arous, G., Gheissari, R., and Jagannath, A. (2022).
High-dimensional limit theorems for sgd: Effective dynamics and
critical scaling.
arXiv preprint arXiv:2206.04030.

[11]
Berthier, R., Montanari, A., and Zhou, K. (2023).
Learning time-scales in two-layers neural networks.
arXiv preprint arXiv:2303.00055.

[12]
Bietti, A., Bruna, J., Sanford, C., and Song, M. J. (2022).
Learning single-index models with shallow neural networks.
In Advances in Neural Information Processing Systems.

[13]
Chen, S., Dou, Z., Goel, S., Klivans, A. R., and Meka, R. (2023).
Learning narrow one-hidden-layer relu networks.
arXiv preprint arXiv:2304.10524.

[14]
Damian, A., Lee, J., and Soltanolkotabi, M. (2022).
Neural networks can learn representations with gradient descent.
In Conference on Learning Theory.

[15]
De Carli, L. (2008).
Local lp inequalities for gegenbauer polynomials.
In Topics in classical analysis and applications in honor of
Daniel Waterman, pages 73–87. World Scientific.

[16]
Diakonikolas, I., Kane, D. M., and Stewart, A. (2017).
Statistical query lower bounds for robust estimation of
high-dimensional gaussians and gaussian mixtures.
In 2017 IEEE 58th Annual Symposium on Foundations of Computer
Science (FOCS), pages 73–84. IEEE.

[17]
Dimitrov, D. K. and Nikolov, G. P. (2010).
Sharp bounds for the extreme zeros of classical orthogonal
Journal of Approximation Theory, 162(10):1793–1804.
Special Issue dedicated to the memory of Borislav Bojanov.

[18]
DLMF (2022).
NIST Digital Library of Mathematical Functions.
<https://dlmf.nist.gov/>, Release 1.1.10 of 2023-06-15.
F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider,
R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and
M. A. McClain, eds.

[19]
Driver, K. and Jordaan, K. (2012).
Bounds for extreme zeros of some classical orthogonal polynomials.
Journal of Approximation Theory, 164(9):1200–1204.

[20]
Dudeja, R. and Hsu, D. (2018).
Learning single-index models in gaussian space.
In Bubeck, S., Perchet, V., and Rigollet, P., editors, Proceedings of the 31st Conference On Learning Theory, volume 75 of Proceedings of Machine Learning Research, pages 1887–1930. PMLR.

[21]
Frei, S., Cao, Y., and Gu, Q. (2020).
Agnostic learning of a single neuron with gradient descent.
Advances in Neural Information Processing Systems,

[22]
Frye, C. and Efthimiou, C. J. (2012).
Spherical harmonics in p dimensions.
arXiv preprint arXiv:1205.3548.

[23]
Gautschi, W. (1959).
Some elementary inequalities relating to the gamma and incomplete
gamma function.
J. Math. Phys, 38(1):77–81.

[24]
Goel, S., Gollakota, A., Jin, Z., Karmalkar, S., and Klivans, A. (2020).
Superpolynomial lower bounds for learning one-layer neural networks
using gradient descent.
In International Conference on Machine Learning, pages
3587–3596. PMLR.

[25]
Kakade, S. M., Kanade, V., Shamir, O., and Kalai, A. (2011).
Efficient learning of generalized linear and single index models with
isotonic regression.
Advances in Neural Information Processing Systems, 24.

[26]
Kalai, A. T. and Sastry, R. (2009).
The isotron algorithm: High-dimensional isotonic regression.
In COLT.

[27]
Klartag, B. (2007).
A central limit theorem for convex sets.
Inventiones mathematicae, 168(1):91–131.

[28]
Laforgia, A. and Natalini, P. (2013).
On some inequalities for the gamma function.
Advances in Dynamical Systems and Applications, 8(2):261–267.

[29]
Niles-Weed, J. and Rigollet, P. (2022).
Estimation of wasserstein distances in the spiked transport model.
Bernoulli, 28(4):2663–2688.

[30]
Paty, F.-P. and Cuturi, M. (2019).
Subspace robust wasserstein distances.
In International conference on machine learning, pages
5072–5081. PMLR.

[31]
Röllin, A. (2013).
Stein's method in high dimensions with applications.
In Annales de l'IHP Probabilités et statistiques,
volume 49, pages 529–549.

[32]
Shalev-Shwartz, S., Shamir, O., and Sridharan, K. (2010).
Learning kernel-based halfspaces with the zero-one loss.
arXiv preprint arXiv:1005.3681.

[33]
Shamir, O. (2018).
Distribution-specific hardness of learning neural networks.
The Journal of Machine Learning Research, 19(1):1135–1163.

[34]
Soltanolkotabi, M. (2017).
Learning relus via gradient descent.
Advances in neural information processing systems, 30.

[35]
Song, M. J., Zadik, I., and Bruna, J. (2021).
On the cryptographic hardness of learning single periodic neurons.
Advances in Neural Processing Systems (NeurIPS).

[36]
Stam, A. J. (1982).
Limit theorems for uniform distributions on spheres in
high-dimensional euclidean spaces.
Journal of Applied probability, 19(1):221–228.

[37]
Szego, G. (1939).
Orthogonal polynomials, volume 23.
American Mathematical Soc.

[38]
Ursell, F. (2007).
Integrals with nearly coincident branch points: Gegenbauer
polynomials of large degree.
Proceedings of the Royal Society A: Mathematical, Physical and
Engineering Sciences, 463(2079):697–710.

[39]
Watson, G. N. (1922).
A treatise on the theory of Bessel functions, volume 3.
The University Press.

[40]
Wu, L. (2022).
Learning a single neuron for non-monotonic activation functions.
In International Conference on Artificial Intelligence and
Statistics, pages 4178–4197. PMLR.

[41]
Yehudai, G. and Shamir, O. (2020).
Learning a single neuron with gradient methods.
In Conference on Learning Theory, pages 3756–3786. PMLR.

$ $

We gather in the appendix the proofs of the theorems, propositions and lemmas stated in the main text. In Section <ref>, the reader will find a short proof of Lemma <ref>. In Section <ref>, we prove Theorems <ref> and <ref> on the SGD dynamics. Sections <ref> and <ref> are respectively devoted to prove that the $\LPG$ property holds in some spherical symmetric case and under some perturbative regime.

§ PROOF OF LEMMA <REF>

Let us first recall the Lemma before writing a proof of it.

For all $a > 0$, we have $ \P_{\theta_0}(m_{\theta_0} \geq a/\sqrt{d}) \leq a^{-1}e^{-a^2/4}  $. Additionally, for any $\delta>0$ such that $ \max\{a, \delta\} \leq \sqrt{d}/4$, we have the lower bound:
$\P_{\theta_0}(m_{\theta_0} \geq a/\sqrt{d}) \geq \frac{\delta}{4} e^ {- (a + \delta)^2} $.

By rotation invariance of the uniform distribution of the sphere, $m_{\theta_0}$ is distributed according to $\theta_0[1]$, the first coordinate of the vector $\theta_0 \in \S$. By a particular case of Stam’s formula <cit.>, we know that for $d \geq 3$, both are distributed according to the probability of density, $\forall  t \in \R$,

\begin{align*}
\tau(t):=\frac{\Gamma(d/2)}{\sqrt{\pi} \Gamma((d-1)/2)} \left(1 - t^2\right)^{(d-3)/2} \mathds{1}_{[-1, 1]}.
\end{align*}

First, note that we can upper and lower bound the constant by the following:

\begin{align*}
\sqrt{\frac{d}{3}} \leq \frac{\Gamma(d/2)}{\Gamma((d-1)/2)} \leq \sqrt{\frac{d}{2}},
\end{align*}

for $d \geq 6$ by <cit.>, which was already proved in [23].

Hence, in terms of the upper bound, we have:

\begin{align*}
\P_{\theta_0}(m_{\theta_0} \geq a/\sqrt{d}) &\leq \sqrt{\frac{d}{2\pi}} \int_{a/\sqrt{d}}^1  \left(1 - t^2\right)^{(d-3)/2} dt \\
&\leq \sqrt{\frac{d}{2 \pi}} \int_{a/\sqrt{d}}^1  e^{-\frac{d-3}{2} t^2} dt \\
&\leq \frac{1}{\sqrt{2\pi}}\frac{d}{a} \int_{a/\sqrt{d}}^1 t  e^{-\frac{d-3}{2} t^2} dt \\
&\leq \frac{1}{2a}   e^{-a^2/4}~,
\end{align*}

which concludes the first part of the result.

Second, let $a \leq \sqrt{d}/4$ and take any $0<\delta< \sqrt{d}/4$. We have,

\begin{align*}
\P_{\theta_0}(m_{\theta_0} \geq a/\sqrt{d}) &\geq \sqrt{\frac{d}{3\pi}} \int_{a/\sqrt{d}}^1  \left(1 - t^2\right)^{(d-3)/2} dt \\
&\geq  \sqrt{\frac{d}{3\pi}} \int_{a/\sqrt{d}}^{(a + \delta)/\sqrt{d}}  \left(1 - t^2\right)^{(d-3)/2} dt \\
&\geq \sqrt{\frac{d}{3\pi}} \frac{\delta}{\sqrt{d}} \left(1 - \frac{(a + \delta)^2}{d} \right)^{(d-3)/2},
\end{align*}

where the last inequality simply comes from the fact that $t \to (1-t^2)^{(d-3)/2}$ is non-increasing. Going further, if we lower bound the term with the negative $-3/2$ power by $1$, we have

\begin{align*}
\P_{\theta_0}(m_{\theta_0} \geq a/\sqrt{d}) &\geq \frac{\delta}{4} \exp\left( \frac{d}{2} \log\left(1 - \frac{(a + \delta)^2}{d} \right)\right) \\
&\geq  \frac{\delta}{4} \exp\left( - \frac{(a + \delta)^2}{2 \left( 1 - (a + \delta)^2/d \right)} \right)~,
\end{align*}

where the last inequality come from the classical bound $\log (1 + x) \geq x/(1+x)$, for $x > - 1$. Furthermore, as, $(a+\delta)^2 \leq d/2$, we have finally

\begin{align*}
\P_{\theta_0}(m_{\theta_0} \geq a/\sqrt{d}) &\geq  \frac{\delta}{4} e^ {- (a + \delta)^2},
\end{align*}

which finalizes the proof of the Lemma.

§ THE LOL PROPERTY IN THE SYMMETRIC CASE: PROOFS OF SECTION <REF>

§.§ Useful Facts about Gegenbauer Polynomials

We recall known facts on Gegenbauer Polynomials.

Recall that $P_{j,d}$ denotes the Gegenbauer polynomial of degree $j$ and dimension $d$, normalized so that $P_{j,d}(1) = 1$ for all $j,d$.
We denote also $\bar{P}_{j,\lambda}$ the Gegenbauer polynomials normalized so that $\| \bar{P}_{j,\lambda} \|^2_{L^2(\mathbb{R}, u_{2\lambda+2})} = \pi 2^{1-2\lambda} \frac{\Gamma(j + 2\lambda)}{(j+\lambda)\Gamma^2(\lambda)\Gamma(j+1)}$. Throughout the proof, we will use either $d$, and from time to time the mute symbol $\lambda$ to denote the dimension variable of Gegenbauer polynomials. They satisfy the following recurrence:

\begin{align}
\label{eq:recurrgegen}
(j+1)\bar{P}_{j+1, \lambda}(t) & = 2(j+\lambda) t \bar{P}_{j,\lambda}(t)-(j+2\lambda-1)\bar{P}_{j-1,\lambda}(t)~,
\end{align}

with first terms: $ \bar{P}_{0, \lambda}(t) = 1$ and $\bar{P}_{1, \lambda}(t)  = 2\lambda t$.

Rodrigues Formula for Gegenbauer Polynomials. The Gegenbauer polynomials can be represented as repeated derivatives of a simple polynome.

We have the formula

\begin{equation}
P_{j,d}(t) = \frac{(-1)^j}{2^j (j + (d-3)/2)_j }(1-t^2)^{(3-d)/2} \left(\frac{d}{dt}\right)^j (1-t^2)^{j+ (d-3)/2}~,
\end{equation}

where $(x)_j=\prod_{k=0}^{j-1} (x-k)$ is the falling factorial.

Hecke-Funk Formula.  Recall that we use the notation $\tau_d$ to denote the uniform distribution on the sphere and $u_d$ the distribution of, e.g., its first coordinate: $u_d \propto (1-t^2)^{(d-3)/2} \mathds{1}_{[-1,1]}$.

For $\theta, \theta' \in \S$, $f \in L^2_{u_d}(\mathbb{R})$ and $j \in \mathbb{N}$,

\begin{align}
\langle f_\theta, (P_{j,d})_{\theta'} \rangle_{\tau_{d}} &= \Omega_{d-2} P_{j,d}( \theta \cdot \theta' ) \langle f, P_{j,d} \rangle_{u_d} \nonumber \\
&= \Omega_{d-2} P_{j,d}( \theta \cdot \theta' ) \int_{-1}^1 f(t) P_{j,d}(t) (1-t^2)^{(d-3)/2} dt ~.
\end{align}

[Derivative Representation] We have the following derivation property for all $j,d$:

\begin{align}
%   P_{j,d}'(m)= \frac{-j(j + d -2)}{(d-1)/2} P_{j-1,d+2}(m)
P_{j,d}'  &= \frac{j (j+d-2)}{(d-1)} P_{j-1,d+2}~.
\end{align}

Recall the normalization relationships   $\lambda = \frac{d}{2}-1$,
$\bar{P}_{j,\lambda}(1) = \frac{\Gamma(j+2\lambda)}{\Gamma(j+1)\Gamma(2\lambda)}$ ,
$P_{j,d} = \frac{\Gamma(j+1)\Gamma(d-2)}{\Gamma(j+d-2)}\bar{P}_{j, \frac{d}{2}-1}$, as well as the identity
$\bar{P}'_{j,\lambda} = 2 \lambda \bar{P}_{j-1,\lambda+1}$. Thus,

\begin{align}
%   P_{j,d}'(m)= \frac{-j(j + d -2)}{(d-1)/2} P_{j-1,d+2}(m)
P_{j,d}' &=  \frac{\Gamma(j+1)\Gamma(d-2)}{\Gamma(j+d-2)}\bar{P}'_{j, \frac{d}{2}-1} \nonumber \\
&=  2\frac{\Gamma(j+1)\Gamma(d-2)}{\Gamma(j+d-2)} (\frac{d}{2}-1) \bar{P}_{j-1,\frac{d}{2}} \nonumber \\
& = (d-2)  \frac{\Gamma(j+1)\Gamma(d-2)}{\Gamma(j+d-2)}  \frac{\Gamma(j-1+d)}{\Gamma(j)\Gamma(d)} P_{j-1,d+2} \nonumber \\
&= \frac{j (j+d-2)}{(d-1)} P_{j-1,d+2}
\end{align}

We have the following bound of the location of the largest root $z_{j,d}$ of $P_{j,d}$:
[Bound on the Largest Root,<cit.>]

\begin{equation}
z_{j,d} \leq \sqrt{\frac{(j-1)(j+d-4)}{(j+d/2-3)(j+d/2-2)}} \cos(\pi/(j+1)) ~.
\end{equation}

And we have the following bound on the Taylor expansion of the Gegenbauer polynomials:
[Taylor Upper bound beyond largest root]

$$P_{j,d}(t) \geq ( t - z_{j,d})^{j}~,~\text{for } t \geq z_{j,d}~,$$

Note that all families of orthogonal polynomials have exclusively real, simple roots.  Therefore, by Rolle's theorem, the $j-1$ critical points of  $P_{j,d}$ must be interlaced with the $j$ zeroes.  So all zeroes of $P_{j,d}'$ are upper bounded by $z_{j,d}$.  Futhermore, by Fact <ref>, $P_{j,d}'$ is itself an orthogonal polynomial.  So applying this argument recursively, we see the zeros of $P_{j,d}^{(k)}$ for $k \leq j$ are all upper bounded by $z_{j,d}$.

Note also by Fact <ref> that, because $P_{j,d}(1) = 1$ for any choice of $j$ and $d$, it follows that $P_{j,d}^{(k)}(1) > 0$.  This implies $P_{j,d}^{(k)}(z_{j,d}) > 0$, as in order to flip signs there would need to be a zero in the range $[z_{j,d}, 1]$ which we've confirmed above cannot exist.

Now, consider a Taylor expansion

\begin{align}
P_{j,d}(t) = \sum_{i=0}^j c_i (t - z_{j,d})^i
\end{align}

Observe that $P_{j,d}^{(k)}(z_{j,d}) = k! c_k$, and therefore by the above argument we have $c_k > 0$.  So it remains to show that $c_j \geq 1$.

Consider applying Fact <ref> repeatedly, then we have:

\begin{align}
P_{j,d}^{(j)}(1) &= \frac{j! \prod_{l=1}^j (j+d-3 + l)}{\prod_{l=1}^j (d-3 + 2l)} \\
&= j! \prod_{l=1}^j \frac{j+d-3 + l}{d-3 + 2l}\\
&\geq j!
\end{align}

And from the fact that $P_{j,d}^{(j)}(1) = j!c_j$, we conclude $c_j \geq 1$.

§.§ Proof of Proposition <ref>

The $\beta_{j,d}$ defined in (<ref>) have the integral representation

\begin{equation}
\beta_{j,d} = \langle \phi, \mathcal{K}_j \phi \rangle_{L^2(\mathbb{R},\eta)} ~,
\end{equation}

where $\mathcal{K}_j$ is a positive semi-definite integral operator of $L^2_\eta$ that depend solely on $\rho$ and $\phi$,
with kernel

\begin{align}
\mathcal{K}_j(t,t') &= \frac{\Omega_{d-2}N(j,d)}{\Omega_{d-1}}  \int_0^\infty  P_j(r^{-1} t) P_j(r^{-1} t') \bar{u}_d(r^{-1} t) \bar{u}_d(r^{-1} t') \rho(dr)~,
\end{align}

where we defined the conditional density

$$\bar{u}_d(r^{-1}t) = \frac{r^{-1} {u}_d(r^{-1} t)}{\int_0^\infty (r')^{-1} {u}_d((r')^{-1} t) \rho(dr')}~. $$

Moreover, we have

\begin{equation}
\mathbb{E}_\eta [\phi^2] = \frac{\Omega_{d-2}}{\Omega_{d-1}}\sum_j \beta_{j,d} = \frac{\Gamma((d-2)/2)}{\sqrt{\pi}\Gamma((d-1)/2)} \sum_j \beta_{j,d}~.
\end{equation}

The marginal conditioned on $\|x\|=r$ is precisely given by
$\eta(x_1=t ~|~\|x\|=r) = r^{-1} u_d(r^{-1} t)$, so

$$\eta(t) = \int_0^\infty r^{-1} u_d(r^{-1} t) \rho(dr)~.$$

We have

\begin{align}
\alpha_{j,r} &= \| P_j\|^{-2} \int_{-1}^1 P_j(t) \phi^{(r)}(t) \tau_d(dt) =  \| P_j\|^{-2} \int_{-1}^1 P_j(t) \phi(r t) u_d(t) dt \nonumber\\
%&= r^{-1}\int_{-r}^r \phi(t) P_j(r^{-1} t) (1-r^{-2}t^2)^{(d-2)/3} dt \\
&= \frac{\Omega_{d-2}N(j,d)}{\Omega_{d-1}} r^{-1} \int_{-\infty}^\infty \phi(t) P_j(r^{-1} t) (1-r^{-2}t^2)^{(d-3)/2}_+ dt~,
\end{align}

\begin{align}
\beta_{j,d} &= \frac{\Omega_{d-2}N(j,d)}{\Omega_{d-1}}  \int_0^\infty r^{-2} \iint_{-\infty}^\infty \phi(t) P_j(r^{-1} t) (1-r^{-2}t^2)^{(d-3)/2}_+ \phi(t') P_j(r^{-1} t') (1-r^{-2}(t')^2)^{(d-3)/2}_+ dt dt' \rho(dr) \nonumber\\
&= \langle \phi, \mathcal{K}_j \phi \rangle_{L^2(\mathbb{R}, \eta)}~,
\end{align}

with the $L^2(\mathbb{R},\eta)$ positive semi-definite integral kernel operator

\begin{align}
\mathcal{K}_j(t,t') &= \frac{\Omega_{d-2}N(j,d)}{\Omega_{d-1}}  \eta(t)^{-1} \eta(t')^{-1} \int_0^\infty r^{-2} P_j(r^{-1} t) (1-r^{-2}t^2)^{(d-3)/2}_+ P_j(r^{-1} t') (1-r^{-2}(t')^2)^{(d-3)/2}_+ \rho(dr) \nonumber \\
&= \frac{\Omega_{d-2}N(j,d)}{\Omega_{d-1}}  \int_0^\infty  P_j(r^{-1} t) P_j(r^{-1} t') \bar{u}_d(r^{-1} t) \bar{u}_d(r^{-1} t') \rho(dr)~,
\end{align}

where we defined the conditional density

$$\bar{u}_d(r^{-1}t) = \frac{r^{-1} {u}_d(r^{-1} t)}{\int_0^\infty (r')^{-1} {u}_d((r')^{-1} t) \rho(dr')}~. $$

Finally, let us establish (<ref>).
We have

\begin{align}
\label{eq:enercons1}
\mathbb{E}_\eta \phi^2 &= \mathbb{E}_\rho [\mathbb{E}_{x_1|\|x\|=r} \mathbb{E} \phi(x_1)^2 ] \nonumber \\
& = \mathbb{E}_\rho [ \mathbb{E}_{u_d} (\phi^{(r)})^2] \nonumber \\
& = \mathbb{E}_\rho \sum_j \alpha_{j,d,r}^2 \| P_j \|^2 \nonumber \\
& = \mathbb{E}_\rho \sum_j \alpha_{j,d,r}^2 \frac{\Omega_{d-1}}{\Omega_{d-2} N(j,d)} \nonumber \\
&= \frac{\Omega_{d-1}}{\Omega_{d-2}} \mathbb{E}_\rho \sum_j \bar{\alpha}_{j,r,d}^2 = \frac{\Omega_{d-1}}{\Omega_{d-2}} \sum_j \beta_{j,d}~.
\end{align}

§.§ Proof of Proposition <ref>

If $\beta_{j,d}=0$ for $j < s$, then $\alpha_{j,r,d}=0$ for $j<s$ and $\rho$-ae $r$. We want to show that for any polynomial $Q$ of degree $j'<s$, we must have
$\langle \phi, Q \rangle_\eta  = 0$.

For each $r$, consider $Q^{(r)}(t) = Q(rt)$, which is also a polynomial of degree $j'<s$, and its decomposition as $Q^{(r)} = \sum_{j=0}^{j'} b_{j,j',r} P_{j,d}$, which only involves terms of degree $j'<s$ since Gegenbauer polynomials of degree up to $r$ span all polynomials of degree up to $r$.
We have

\begin{align}
\langle \phi, Q \rangle_\eta &= \mathbb{E}_\eta [ \phi(x) Q(x)] \nonumber \\
&= \mathbb{E}_\rho \mathbb{E}_{x_1 | \|x\|=r} [\phi(x) Q(x)] \nonumber \\
&= \mathbb{E}_\rho \mathbb{E}_{u_d} [\phi^{(r)}(x) Q^{(r)}(x)] \nonumber \\
&= \mathbb{E}_\rho [\sum_{j \leq j'} b_{j,j',r} \alpha_{j,r,d} ] = 0~.
\end{align}

§.§ Proof of Proposition <ref>

Assume first that there are $\bar{C}, \bar{\zeta}$ such that

\begin{equation}
\label{eq:bi0}
P'_{s,d}(t) \geq \bar{C}( t - \bar{\zeta})^{s-1}~,~\text{ for } t \geq \bar{\zeta}~.
\end{equation}

Now, let

\begin{equation}
\label{eq:bi1}
B = \frac{1}{d-1} \sum_{j \geq s} \beta_{j,d} j(j+d-2)  \upsilon_{j-1,d+2} < 1~,
\end{equation}

and define

\begin{equation}
\label{eq:bi2}
\zeta^* := \left( \frac{B}{\beta_{s,d}\bar{C}}\right)^{1/(s-1)} + \bar{\zeta}~.
\end{equation}

From (<ref>), (<ref>) and (<ref>) we verify that $\ell'(m) = \sum_j \beta_{j,d} P'_{j,d}(m)$ satisfies, for $m \geq \zeta^*$,

$$\ell'(m) \geq \beta_{s,d} \bar{C} \left((m-\bar{\zeta})^{s-1} - (m-\zeta^*)^{s-1} \right) \geq  \beta_{s,d} \bar{C} \left[\left(\frac{1-\bar{\zeta}}{1-\zeta^*} \right)^{s-1} -1 \right] (m-\zeta^*)^{s-1}~.$$

Finally, we have that for any $j,d$, the largest root $z_{j,d}$ satisfies $z_{j,d} \leq \sqrt{\frac{(j-1)(j+2d-2)}{(j+d-2)(j+d-1)}} \simeq j/\sqrt{d}$ and

$$P_{j,d}(t) \geq \frac12 ( t - z_{j,d})^{s}~,~\text{for } t \geq z_{j,d}~,$$

which implies that

\begin{equation}
P'_{s,d}(t) \geq \frac{s(s+d-2)}{2(d-1)} ( t - z_{s-1,d+2})^{s-1} ~,~\text{ for } t \geq z_{s-1,d+2}~.
\end{equation}

We thus have $\bar{C} = \frac{s(s+d-2)}{2(d-1)}$ with $\bar{\zeta} = z_{s-1,d+2}$.

Finally, we verify that

\begin{align}
\label{eq:bi3}
\frac{d^{(s-1)/2}}{s(s+d-2)} \sum_{j > s} \beta_{j,d} j (j+d-2) \upsilon_{j-1,d+2} &\leq K
\end{align}

ensures a local polynomial growth of order $s-1$ at scale $O(1/\sqrt{d})$.
Indeed, plugging (<ref>) into (<ref>), together with $\beta_{s,d} \geq C$ yields

\begin{align}
\left( \frac{B}{\beta_{s,d}\bar{C}}\right)^{1/(s-1)} &\leq (CK)^{1/(s-1)}d^{-1/2} ~,
\end{align}

which shows that $\zeta^* = O(1/\sqrt{d})$.
Finally, we observe that $\bar{C} \geq s=\Theta(1)$ if $s<d$.

§.§ Proof of Theorem <ref>

To prove the theorem, we will establish the sufficient conditions of Proposition <ref> under our mild assumptions.
The key technical results we need are explicit bounds for $\upsilon_{j,d}$ and for the sum $\sum_j j^2 \beta_{j,d}$, established in the following two lemmas. Since the parameter $\lambda = d/2-1$ is more convenient to express many relationships in Gegenbauer polynomials, we will adopt it in this proof instead of $d$, without loss of generality.

We have

\begin{equation}
\upsilon_{j,\lambda} \lesssim  \begin{cases}
\left[1 - \left(\frac{\lambda}{j+\lambda} \right)^2 \right]^{j/2}
%O\left( \left(\frac{j}{\lambda}\right)^{j/2} \right)
& \text{ if } j=\Omega(1) ~,\\
\lambda^{\frac{(\alpha-1)\lambda^\alpha}{2}} & \text{ if } j = \Theta(\lambda^\alpha), \text{ with } 0< \alpha < 1~, \\
e^{-\frac12 \lambda^{2-\alpha}} & \text{ if } j=\Theta(\lambda^\alpha) \text{ with } 1\leq \alpha < 3/2 ~, \\
e^{-\lambda} & \text{ if } j = \Omega(\lambda^{3/2}) ~.
\end{cases}
\end{equation}

If $\phi \in L^2(\mathbb{R}, \mu)$ is such that $\phi' \in L^4(\mathbb{R}, \eta)$ and $\mathbb{E}_\rho [r^4] < \infty$, then
$\beta_j = \langle \phi, \mathcal{K}_j \phi \rangle$ satisfies

\begin{equation}
\sum_j j^2 \beta_{j,d} \leq \frac{\Omega_{d-2}}{\Omega_{d-1}} \mathbb{E}_\rho [r^4]^{1/2} \|\phi'\|_{L^4(\eta)}^2 = O(1/d)~.
\end{equation}

Let $s = \inf\{ j ; \beta_{j,d} \neq 0\}$.
We need to verify that
there exists a constant $K>0$ such that

\begin{equation}
\label{eq:flower}
\sum_{j > s} \beta_{j,d} j (j+d-2)  \upsilon_{j-1,d+2} \leq K d^{(3-s)/2}~.
\end{equation}

We will control the LHS by splitting it into appropriate regions, determined by $J_i$, $i\in \{1,2,3\}$.
Let $\alpha = \frac{4}{1+s}$ and $J_1 = \frac{\lambda^\alpha}{2}$.
From Lemma <ref>, part (i) we have that
$\upsilon_{j,\lambda} \leq C \left( \frac{j(j+2\lambda)}{(j+\lambda)^2}\right)^{j/2}$,
and in particular $\upsilon_{j,\lambda} \leq C \lambda^{(\alpha-1)j/2}$ for $j \leq J_1$.
As a result, using Lemma <ref>,

\begin{align}
\sum_{j = s+1}^{J_1} \beta_j j (j+\lambda)  \upsilon_{j-1,\lambda+1} & \leq \lambda^{(\alpha-1)(s+1)/2} \sum_{j=s+1}^{J_1} \beta_j j (j+\lambda) \nonumber\\
&\leq  \lambda^{(\alpha-1)(s+1)/2} (C_1 \lambda^{-1} + \lambda \sum_{j=s+1} ^{J_1} \beta_j j )\nonumber\\
&\leq  \lambda^{(\alpha-1)(s+1)/2} (C_1 \lambda^{-1} + \lambda \sum_{j=s+1} ^{J_1} \beta_j j^2 )\nonumber\\
& \leq \lambda^{(\alpha-1)(s+1)/2} C_2 \nonumber\\
&\leq C_2 \lambda^{(3-s)/2}~.
\end{align}

Let $J_2 = \lambda$.
We have

\begin{align}
\sum_{j=J_1+1}^{J_2} \beta_j j (j+\lambda) \upsilon_{j-1,\lambda+1} &\leq \lambda^{\frac{(\alpha-1)\lambda^\alpha}{2}} C_3  \nonumber \\
& \leq C_3 \lambda^{(3-s)/2}~.
\end{align}

Let $J_3 = \lambda^{3/2}$.
We have

\begin{align}
\sum_{j=J_2+1}^{J_3} \beta_j j (j+\lambda) \upsilon_{j-1,\lambda+1} &\leq e^{-\frac12 \sqrt{\lambda}} C_4  \nonumber \\
& \leq C_4 \lambda^{(3-s)/2}~.
\end{align}

Finally, the remainder satisfies

\begin{align}
\sum_{j>J_3} \beta_j j (j+\lambda) \upsilon_{j-1,\lambda+1} &\leq C_5 (e/2)^{-{\lambda}} \nonumber \\
& \leq C_5 \lambda^{(3-s)/2}~,
\end{align}

which proves (<ref>).

We prove this result by analysing different regimes for $j$ and $\lambda$. Concretely, we claim the following:

We have the following regimes:

* For $j = \Omega(1)$, we have

\begin{equation}
\label{eq:upsi_mainclaim1}
\upsilon_{j,\lambda} \lesssim \left[1 - \left(\frac{\lambda}{j+\lambda} \right)^2 \right]^{j/2}~.
%  \left(\frac{2j}{\lambda} \right)^{j/2}~.
\end{equation}

* For $j = \Theta(\lambda^\alpha)$, with $0< \alpha < 1$, we have

\begin{equation}
\label{eq:upsi_mainclaim2}
\upsilon_{j,\lambda} \lesssim \lambda^{\frac{(\alpha-1)\lambda^\alpha}{2}}~.
\end{equation}

* For $j = \Theta(\lambda^\alpha)$, with $1 \leq \alpha < 2$, we have

\begin{equation}
\label{eq:upsi_mainclaim2bis}
\upsilon_{j,\lambda} \lesssim e^{-\frac12 \lambda^{2-\alpha}}~.
\end{equation}

* For $j = \Omega(\lambda^\alpha)$, with $\alpha > 3/2$, we have

\begin{equation}
\label{eq:upsi_mainclaim3}
\upsilon_{j,\lambda} \lesssim e^{-\lambda}~.
\end{equation}

To prove the first three regimes of Claim <ref>, we control $\upsilon_{j,\lambda}$ based on the distribution of the roots of $P_{j,\lambda}$.
We recall that $(z_{k,j,\lambda})_{k \leq j}$ denotes the roots of $P_{j,\lambda}$ in increasing order, and $z_{j,\lambda} = z_{j,j,\lambda}$ its largest root.

We have

\begin{equation}
P_{j,\lambda}(t) = \begin{cases}
\prod_{k=j/2}^j \frac{t^2 - z_{k,j,\lambda}^2}{1 - z_{k,j,\lambda}^2} & \text{ if } j \text{ even}, \\
t \prod_{k=(j+1)/2}^j \frac{t^2 - z_{k,j,\lambda}^2}{1 - z_{k,j,\lambda}^2} & \text{ if } j \text{ odd }~.
\end{cases}
\end{equation}

From this representation, we deduce that $\upsilon_{j,\lambda}$ can be calculated explicitly. Indeed, as the local maxima of $|P_{j,\lambda}(t)|$ are increasing [37], <cit.>, we have the following equation:

\begin{equation}
\label{eq:upsi}
\upsilon_{j,\lambda} = - P_{j,\lambda}(z_{j-1,\lambda+1}) = - \begin{cases}
\prod_{k=j/2}^j \frac{z_{j-1,\lambda+1}^2 - z_{k,j,\lambda}^2}{1 - z_{k,j,\lambda}^2} & \text{ if } j \text{ even}, \\
z_{j-1,\lambda+1} \prod_{k=(j+1)/2}^j \frac{z_{j-1,\lambda+1}^2 - z_{k,j,\lambda}^2}{1 - z_{k,j,\lambda}^2} & \text{ if } j \text{ odd }~.
\end{cases}
\end{equation}

Let us focus first on the case $j$ even, for simplicity.
We can rewrite (<ref>) more conveniently as

$$\upsilon_{j,\lambda} = \frac{z_{j,\lambda}^2 - z_{j-1,\lambda+1}^2}{1-z_{j,\lambda}^2} \prod_{k=j/2}^{j-1} \frac{z_{j-1,\lambda+1}^2 - z_{k,j,\lambda}^2}{1 - z_{k,j,\lambda}^2}~.$$

For $\delta\in (0,z_{j-1,\lambda+1})$ let

$$m(\delta, j, \lambda) := | \{ k \in \{j/2,j\} ; z_{k,j,\lambda} \geq \delta \} |$$

denote the number of zeros of $P_{j,\lambda}$ in the interval $(\delta, 1)$.

Since the function $t \mapsto \frac{a^2 - t^2}{1-t^2}$ is decreasing in $t \in (0,a)$,
we have

We have the upper bound:

\begin{equation}
\upsilon_{j,\lambda} \leq  \frac{z_{j,\lambda}^2 - z_{j-1,\lambda+1}^2}{1-z_{j,\lambda}^2} \inf_{\delta} \left(\frac{z_{j-1,\lambda+1}^2 - \delta^2}{1-\delta^2} \right)^{m(\delta, j, \lambda)}~.
\end{equation}

Letting $\delta = z_{j/2,j,\lambda}$ the smallest positive root of $P_{j,\lambda}$ we have

\begin{equation}
\label{eq:bhi}
\upsilon_{j,\lambda} \leq \frac{z_{j,\lambda}^2 - z_{j-1,\lambda+1}^2}{1-z_{j,\lambda}^2} \left(\frac{z_{j-1,\lambda+1}^2 - z_{j/2,j,\lambda}^2}{1-z_{j/2,j,\lambda}^2} \right)^{j/2}~.
\end{equation}

We can thus obtain an explicit control on $\upsilon_{j,\lambda}$ from bounds on the zeros of the Gegenbauer polynomials. We complement the upper bound on the largest root (Fact <ref>) with lower bounds for all positive roots, as well as a sharp lower bound for its largest root [17]:

$$b_{j,\lambda} = j^3 + 2(\lambda - 1)j^2 - (3\lambda -5)j + 4(\lambda - 1)~,$$

$$a_{j,\lambda} = 2(j + \lambda-1)(j^2 + j(\lambda-1) + 4(\lambda+1))~\text{and}$$

$$c_{j,\lambda} = j^2(j+2\lambda)^2 + (2\lambda + 1)(j^2 + 2(\lambda+3)j + 8(\lambda - 1))~.$$
Then for every $k,j,\lambda$ we have

\begin{equation}
\frac{b_{j,\lambda} - (j-2)\sqrt{c_{j,\lambda}}}{a_{j,\lambda}} \leq z_{k,j,\lambda}^2 \leq \frac{b_{j,\lambda} + (j-2)\sqrt{c_{j,\lambda}}}{a_{j,\lambda}}~.
\end{equation}

\begin{equation}
z_{j,\lambda}^2 > 1 - \frac{(2\lambda+1)(2\lambda+3)}{(j-1)(j+2\lambda+1)+(2\lambda+1)(2\lambda+3)}:= 1 - \frac{g_{j,\lambda}}{h_{j,\lambda}}~.
\end{equation}

Rewriting Fact <ref> as
$z_{j,\lambda}^2 \leq \frac{e_{j,\lambda}}{f_{j,\lambda}}~,$

$$e_{j,\lambda} = (j-1)(j+2\lambda-2)~,~f_{j,\lambda} = (j+\lambda-2)(j+\lambda-1)~,$$

and using again the monotonocity of $t \mapsto \frac{t-p}{1-t}$
we can bound the first term in the RHS of (<ref>) as

\begin{equation}
\label{eq:bhi2}
\frac{z_{j,\lambda}^2 - z_{j-1,\lambda+1}^2}{1-z_{j,\lambda}^2} \leq \frac{e_{j,\lambda}/f_{j,\lambda} + g_{j-1,\lambda+1}/h_{j-1,\lambda+1} -1}{1-e_{j,\lambda}/f_{j,\lambda}}~.
\end{equation}

For $j, \lambda = \omega(1)$,
we have

\begin{align*}
a_{j,\lambda} &\simeq 2j (j+\lambda)^2,& ~b_{j,\lambda} &\simeq j^2(j+2\lambda)~,~\sqrt{c_{j,\lambda}}\simeq j(j+2\lambda)~, \\
e_{j,\lambda} &\simeq j(j+2\lambda),&~f_{j,\lambda}&\simeq (j+\lambda)^2~,\\
g_{j,\lambda} &\simeq 4\lambda^2 ,&~ h_{j,\lambda} &\simeq j(j+2\lambda) + 4\lambda^2~,
\end{align*}

and thus

\begin{equation}
\label{eq:bhi2}
\frac{z_{j,\lambda}^2 - z_{j-1,\lambda+1}^2}{1-z_{j,\lambda}^2} \lesssim \frac{3 j (j+2\lambda)}{j (j+2\lambda) + 4\lambda^2} \leq 3~.
\end{equation}

\begin{align}
\label{eq:upsi_smalljregime}
\upsilon_{j,\lambda} & \leq 3 \left(\frac{\frac{e_{j-1,\lambda+1}}{f_{j-1,\lambda+1}} - \frac{b_{j,\lambda} - (j-2)\sqrt{c_{j,\lambda}}}{a_{j,\lambda}}}{1-\frac{b_{j,\lambda} - (j-2)\sqrt{c_{j,\lambda}}}{a_{j,\lambda}}} \right)^{j/2}~ \nonumber \\
%&= \left( 1 + \frac{d_{j-1,\lambda+1} f_{j,\lambda}}{2  c_{j-%1,\lambda+1}  (e_{j,\lambda}- f_{j,\lambda})}\right)
%\left(\frac{a_{j,\lambda} e_{j-1,\lambda+1} - f_{j-1,\lambda+1}%(b_{j,\lambda} - (j-2)\sqrt{c_{j,\lambda}})}{f_{j-1,\lambda+1}(a_{j,\lambda} - b_{j,\lambda} + (j-2)\sqrt{c_{j,\lambda}})}\right)^{j/2} \nonumber \\
& \leq 3\left(\frac{a_{j,\lambda} e_{j-1,\lambda+1} - f_{j-1,\lambda+1}(b_{j,\lambda} - (j-2)\sqrt{c_{j,\lambda}})}{f_{j-1,\lambda+1}(a_{j,\lambda} - b_{j,\lambda} + (j-2)\sqrt{c_{j,\lambda}})}\right)^{j/2}~ \nonumber \\
& = 3\left(\frac{2j(j+\lambda)^2 j(j+2\lambda) - (j+\lambda)^2 (j^2(j+2\lambda) - j^2(j+2\lambda)) }{(j+\lambda)^2(2j(j+\lambda)^2 - j^2(j+2\lambda) + j^2(j+2\lambda))}  \cdot (1 + o_{j,\lambda}(1)) \right)^{j/2} \nonumber \\
& \lesssim \left(\frac{j(j+2\lambda)}{(j+\lambda)^2} \right)^{j/2} \nonumber \\
&= \left[1 - \left(\frac{\lambda}{j+\lambda} \right)^2 \right]^{j/2}~.
% & \leq \frac{z_{j,d}^2 - 1 + d_{j-1,\lambda+1}/(2 c_{j-1,\lambda+1}) }{1-z_{j,d}^2} \left(\frac{z_{j-1,d+2}^2 - (b_{j,\lambda} - (j-2)\sqrt{\delta_{j,\lambda}})/a_{j,\lambda}}{1-(b_{j,\lambda} - (j-2)\sqrt{\delta_{j,\lambda}})/a_{j,\lambda}} \right)^{j/2}
\end{align}

As a direct consequence of (<ref>), we immediately obtain Eqs (<ref>), (<ref>) and (<ref>). The case where $j$ is odd is treated analogously.

Let us now study the regime $j = \omega(\lambda^{3/2})$.

Given $z \in \mathbb{C}$ with $|z|<1$, Gegenbauer polynomials admit the following generating function <cit.>:

\begin{equation}
\frac{1}{(1 - 2z \cos \theta + z^2)^\lambda } = \sum_{j \geq 0} \bar{P}_{j, \lambda}(\cos \theta) z^j ~.
\end{equation}

From this generating function, the Cauchy integral formula leads to
the following integral representation:
For any $0 < \rho < 1$, we have

\begin{equation}
\bar{P}_{j, \lambda}(\cos \theta) = \frac{1}{2\pi i} \oint_{|z|=\rho}  \frac{dz}{(1-2z \cos \theta + z^2)^{\lambda} z^{j+1}}~.
\end{equation}

Assume $j = \Theta(\lambda^\alpha)$, with $\alpha > 3/2$.
We are interested in the above representation for  $\bar{\theta} = \arccos(z_{j-1,\lambda+1})$.
From Theorem <ref>, we have
$z_{j-1,\lambda+1}^2 \geq 1 - d_{j,\lambda}/(2c_{j,\lambda})$, and thus

$$\bar{\theta}^2 \lesssim \frac{d_{j,\lambda}}{2c_{j,\lambda}} \simeq \frac{32\lambda^2 j^4 }{16 j^6} = \frac{2\lambda^2}{j^2}~,$$

so $\bar{\theta} = O(\lambda/j)$.
Combining this upper bound with the lower bound obtained from Fact <ref> we have $\bar{\theta} = \Theta(\lambda/j)$.

Using $1-\cos \theta \simeq  \theta^2/2 \simeq  \lambda^2/j^2$

\begin{align}
|1 - 2 z \cos \theta + z^2| & = \left|(1-z)^2 + 2z(1-\cos \theta) \right| \nonumber \\
&\geq | 1- z|^2 - 2|z|(1-\cos \theta) \nonumber \\
& \geq 1 - \rho \left(2 + \Theta\left(\frac{\lambda^2}{j^2}\right)\right) + \rho^2~,
\end{align}

we have

\begin{align}
| \bar{P}_{j, \lambda}(\cos \theta) | \leq \inf_{0<\rho < 1}|\rho|^{-(j+1)} \left(1 - \rho(2 + c\lambda^2/j^2) + \rho^2 \right)^{-\lambda}:= g(\rho)~.
\end{align}

Optimizing the RHS over $\rho$ we obtain
$\rho^* = \frac{j - (\sqrt{2}-1)\lambda}{j + 2\lambda}$; substituting, we obtain

\begin{align}
g(\rho^*) \simeq e^{-(1+\sqrt{2}) \lambda} \left(\frac{j+2\lambda}{\lambda(1+\sqrt{2})} \right)^{2\lambda}~.
\end{align}

As a result, it follows that

\begin{align}
P_{j, \lambda}(\cos \theta) &= \bar{P}_{j, \lambda}(\cos \theta) \frac{j! (2\lambda-1)!}{(2\lambda+j-1)!}
\end{align}

satisfies, for $\theta = \Theta(\lambda/j)$ and $j = \omega(\lambda^{3/2})$,

\begin{align}
\log |P_{j, \lambda}(\cos \theta) | \simeq & j \log j - j + 2\lambda \log(2\lambda) - 2\lambda - (j +2\lambda) \log(2\lambda + j) + 2\lambda+j \nonumber \\
&~ -(1+\sqrt{2}) \lambda + 2\lambda \log (j + 2\lambda) - 2\lambda \log(\lambda (1+\sqrt{2})) \nonumber \\
\simeq & -(1+\sqrt{2}) \lambda~,
\end{align}

where we have used Stirling's approximation.
This proves Eq (<ref>) and completes the proof of  Lemma <ref>.

We have, using Fact <ref>, that

\begin{align}
\label{eq:bound0}
\mathbb{E}_\rho [ \mathbb{E}_{u_d} {(\phi^{(r)}}')^2]
&= \frac{\Omega_{d-1}}{(d-1)^2\Omega_{d-2}} \sum_j \mathbb{E}_\rho \left[  \bar{\alpha}_{j,d,r}^2 \left(j (j+d-2)\right)^2  \right] \nonumber\\
&\geq \frac{\Omega_{d-1}}{\Omega_{d-2}} \sum_j j^2 \mathbb{E}_\rho \left[  \bar{\alpha}_{j,d,r}^2 \right]~.
\end{align}

And we can upper bound via

\begin{align}
\mathbb{E}_\rho [ \mathbb{E}_{u_d} {(\phi^{(r)}}')^2] & = \mathbb{E}_\rho [ r^2 \mathbb{E}_{u_d} ({(\phi')^{(r)}})^2] \nonumber \\
&= \mathbb{E}_\rho \mathbb{E}_{x_1|\|x\|=r} [ r^2 (\phi(x_1)')^2 ] \nonumber \\
& \leq \sqrt{\mathbb{E}_\rho [r^4] \mathbb{E}_\eta (\phi')^4} ~,
\end{align}

where this last line is finite by our assumptions on $\phi$ and $\rho$,

so from (<ref>) we conclude that

\begin{equation}
\sum_j j^2 \beta_{j,d} \leq \frac{\Omega_{d-2}}{\Omega_{d-1}} \sqrt{\mathbb{E}_\rho [r^4] \mathbb{E}_\eta (\phi')^4}~.
\end{equation}

§ THE LOL PROPERTY IN THE NON-SYMMETRIC CASE: PROOFS OF SECTION <REF>

§.§ Proof of Proposition <ref>

Recall the notation $\phi_\theta(x) = \phi(\langle x, \theta\rangle)$. Let $v = \theta^* - m \theta$.
From the definition, we have that

\begin{align}
\langle \nabla_\theta^{\mathbb{S}} L(\theta), \theta^* \rangle &= 2 \mathbb{E}_\nu\left[ \phi'_\theta(\phi_\theta - \phi_{\theta^*}) \left(x \cdot v\right)\right] \nonumber \\
&:= \mathbb{E}_\nu [g_{\theta, \theta^*}]~.
\end{align}

Since $\mathbb{E}_\gamma [g_{\theta, \theta^*}]$ is precisely $\bar{\ell}'(m)(1-m^2)$, we need
to establish that

\begin{equation}
\sup_\theta \left|\mathbb{E}_\nu g_{\theta, \theta^*} - \mathbb{E}_\gamma g_{\theta, \theta*} \right| \leq C \sqrt{1-m^2} \widetilde{W}_{1,2}(\nu, \gamma) \left(\log  \widetilde{W}_{1,2}(\nu, \gamma) \right)^2~.
\end{equation}

Fix $\theta$ and let $P_{\theta, \theta^*}$ be the orthogonal projection onto the subspace spanned by $\theta, \theta^*$.  For $R>0$ we consider $A_R = \{ x \in \mathbb{R}^d; \| P_{\theta, \theta^*} x \| \leq R \}$.

\begin{align}
\left|\mathbb{E}_\nu g_{\theta, \theta^*} - \mathbb{E}_\gamma g_{\theta, \theta*} \right| &= \left| \int g_{\theta, \theta^*}(x) (\nu(dx) - \gamma(dx) ) \right| \nonumber \\
&\leq \underbrace{\left| \int_{x \in A_R} g_{\theta, \theta^*}(x) (\nu(dx) - \gamma(dx) ) \right|}_{T_a} + \underbrace{\left| \int_{x \notin A_R} g_{\theta, \theta^*}(x) (\nu(dx) - \gamma(dx) ) \right|}_{T_b}~.
\end{align}

Let us first bound $T_a$. Denote by $v = \theta^* - m \theta$, with $\|v\|^2 = 1-m^2$ Since $\phi$ and $\phi'$ are Lipschitz and $|\phi''| \leq O((1+t)^{-1})$ by Assumption <ref>, we have that

\begin{align}
\nabla_x g_{\theta, \theta^*}(x) &= \phi''_{\theta} (\phi_\theta - \phi_{\theta^*}) x^\top v \theta + \phi'_\theta ( \phi'_\theta \theta - \phi'_{\theta^*} \theta^*) x^\top v + \phi'_\theta (\phi_\theta - \phi_{\theta^*}) v
\end{align}

\begin{align}
\| \nabla_x g_{\theta, \theta^*}(x) \| &\leq  2\|v\| C \mathrm{Lip}(\phi) R +  4 \|v\|\mathrm{Lip}(\phi)^2 R \nonumber \\
&\leq C \|v\| R ~,
\end{align}

and as a result we have that $g_{\theta, \theta^*}$ is $C \|v\| R$-Lipschitz when restricted to $A_R$, and thus

\begin{equation}
T_a \leq C R \|v\| \widetilde{W}_{1,2}(\nu, \mu)~.
\end{equation}

Let us now control the tail $T_b$. Since $x^\top v$ is $\sqrt{2} M \|v\|$-subgaussian and $\phi$ is Lipschitz, we have that $z = |g_{\theta, \theta^*}(x)|$ is $\tilde{M} \|v\|$-subexponential where $\tilde{M}$ only depends on $M$ and $L$.
It follows that

\begin{align}
T_b &\leq R (\mathbb{P}_\nu(z \geq R) + \mathbb{P}_\gamma( z \geq R) ) \nonumber \\
&\leq R \exp\left( - \frac{\beta}{\|v\|} R \right)~,
\end{align}

where $\beta$ is a constant that depends only on $\tilde{M}$.
As a result, we have

\begin{align}
\left|\mathbb{E}_\nu g_{\theta, \theta^*} - \mathbb{E}_\gamma g_{\theta, \theta*} \right| \leq \inf_{R>0} \left(C R \|v\| \widetilde{W}_{1,2}(\nu, \mu) + R \exp\left( -\frac{\beta}{\|v\|} R \right) \right)~.
\end{align}

$$R =  {-\|v\| \beta^{-1}} \log((C \|v\| \widetilde{W}_{1,2}(\nu, \gamma))) $$

we obtain

\begin{align}
\label{eq:blink}
\left|\mathbb{E}_\nu g_{\theta, \theta^*} - \mathbb{E}_\gamma g_{\theta, \theta*} \right| &\leq  2 C (1-m^2)  \beta^{-1} \left|\log(C \|v\| \widetilde{W}_{1,2}(\nu, \gamma))\right|  \widetilde{W}_{1,2}(\nu, \gamma) \nonumber \\
& \leq (1-m^2) O\left( \widetilde{W}_{1,2}(\nu, \gamma) \log ( \widetilde{W}_{1,2}(\nu, \gamma)^{-1} ) \right)~,
% \left|\mathbb{E}_\nu g_{\theta, \theta^*} - \mathbb{E}_\gamma g_{\theta, \theta*} \right| &\leq \beta^{-2} \|v\|^3 C \widetilde{W}_{1,2}(\nu, \gamma) \log\left(\frac{\beta}{C\|v\|^2 \widetilde{W}_{1,2}(\nu, \gamma)} \right)^2 \nonumber \\
% &= \beta^{-2} (1-m^2)^{3/2} C \widetilde{W}_{1,2}(\nu, \gamma) \log\left(\frac{\beta}{C(1-m^2) \widetilde{W}_{1,2}(\nu, \gamma)} \right)^2~,
\end{align}

as claimed.

§.§ Proof of Proposition <ref>

We leverage Proposition <ref> and the fact that if
$\phi$ has information exponent $s=2$, then $\bar{\ell}'(m) \simeq m$ for small $m$.

We need to show that for $b = \Theta(\log d)$ we have

\begin{equation}
\label{eq:lpgcond0}
\langle \nabla_\theta L(\theta), \theta^*\rangle \geq C \left( m - \frac{b}{\sqrt{d}}\right)~,~\text{for }  \frac{b}{\sqrt{d}} \leq m \leq \frac12~,
\end{equation}

as well as

\begin{equation}
\label{eq:contractcond0}
\langle \nabla_\theta L(\theta), \theta^*\rangle \geq C' (1-m^2)
\end{equation}

for $m\geq \frac12$.

From (<ref>) and $\widetilde{W}_{1,2}(\nu, \gamma) \leq C /\sqrt{d}$,
we obtain

\begin{align}
\langle \nabla_\theta L(\theta), \theta^*\rangle &= \bar{\ell}'(m)(1-m^2) + \langle \nabla_\theta L(\theta), \theta^*\rangle - \bar{\ell}'(m)(1-m^2)   \nonumber \\
&\geq 2 \alpha_2^2 m (1-m^2) - (1-m^2)\tilde{C} \widetilde{W}_{1,2}(\nu, \gamma) \log(\widetilde{W}_{1,2}(\nu, \gamma)^{-1}) \nonumber \\
& \geq  \left(\alpha_2^2 m  - \tilde{C} \frac{C}{\sqrt{d}} \log(\sqrt{d}/C)\right)(1-m^2) \nonumber \\
& \geq  \alpha_2^2 \left( m  - \tilde{C} \frac{C}{\alpha_2^2 \sqrt{d}} \log(\sqrt{d}/C)\right)(1-m^2) \nonumber \\
& \geq \alpha_2^2 \left( m - \frac{\log d^{C'}}{\sqrt{d}} \right)(1-m^2)~,
\end{align}

which proves (<ref>) and (<ref>).

§.§ Proof of Proposition <ref>

Recall the notation $\phi_\theta(x) = \phi( \langle x, \theta\rangle)$, and,  using $v = \theta^* - m \theta$,

\begin{align}
h_{\theta, \theta^*}(x) &:= \phi_\theta^2 - 2 \phi_\theta \phi_{\theta^*}~,\\
g_{\theta, \theta^*}(x) &:= 2 \phi'_\theta ( \phi_\theta - \phi_{\theta^*}) ( x \cdot v )~,
\end{align}

so that

\begin{align}
\Delta_L(\theta) &= \mathbb{E}_\nu[ h_{\theta, \theta^*}(x)] - \mathbb{E}_\gamma[h_{\theta, \theta^*}(x)]~, \\
\Delta_{\nabla L}(\theta) &= \mathbb{E}_\nu[ g_{\theta, \theta^*}(x)] - \mathbb{E}_\gamma[g_{\theta, \theta^*}(x)]~.
\end{align}

The result is obtained via the following Stein coupling method for product measures:

Let $X$ be a $d$-dimensional random vector of independent coordinates, such that $\mathbb{E} X = 0$, $\mathbb{E} [XX^\top] = I_d$ and $\mathbb{E}|X_i|^3 = \tau_i^3 < \infty$. If $Z$ is a standard Gaussian random vector, and $h: \mathbb{R}^d \to \mathbb{R}$ is three-times differentiable, then

\begin{equation}
| \mathbb{E} h(X) - \mathbb{E} h(Z) | \leq \frac{5}{6} \sum_{i=1}^d \tau_i^3 \| \partial_{x_i}^3 h\|_\infty~.
\end{equation}

We verify that, thanks to the decay assumptions in Assumption <ref>, we have

\begin{align}
\partial_{x_i}^3 g_{\theta, \theta^*}(x) &= \lambda_1(x) \theta_i^3 + \lambda_2(x) \theta_i^2 \theta^*_i + \lambda_3(x) \theta_i (\theta^*_i)^2 + \lambda_4(x) (\theta^*_i)^3~, \\
\partial_{x_i}^3 h_{\theta, \theta^*}(x) &= \lambda_5(x) \theta_i^3 + \lambda_6(x) \theta_i^2 \theta^*_i + \lambda_7(x) \theta_i (\theta^*_i)^2 + \lambda_8(x) (\theta^*_i)^3~,
\end{align}

\begin{equation}
\sup_{k \in \{1,2,3,4\}} |\lambda_k(x)| \leq C \|v\|~,~ \sup_{k \in \{5,6,7,8\}} |\lambda_k(x)| \leq \tilde{C}~.
\end{equation}

Observing by Cauchy-Schwartz that

\begin{align*}
\max\left\{\sum_i |\theta_i|^2 |\theta_i^*| , \sum_i |\theta_i|^3 \right\} &\leq \|\theta\|_4^2 ~,\\
\max\left\{\sum_i |\theta^*_i|^2 |\theta_i| , \sum_i |\theta^*_i|^3 \right\} &\leq \|\theta^*\|_4^2~,
\end{align*}

we obtain from Theorem <ref> that

$$\Delta_L(\theta) \leq C \chi(\theta, \theta^*)~,~\Delta_{\nabla L}(\theta) \leq C' \|v\| \chi(\theta, \theta^*)~,$$

as claimed.