1 Introduction↩︎

Single-Index Models (SIMs) are a broad family of semi-parametric models subsuming linear regression, logistic regression, phase retrieval, and generalized linear models as special cases. They model the response variable \(Y\in\mathbb{R}\) as a nonlinear function of a one-dimensional projection of the covariates \(X\in\mathbb{R}^d\): \[\label{eq:SIM} Y = f(X^\top\beta^\star) + \zeta,\tag{1}\] where \(f:\mathbb{R}\to\mathbb{R}\) is a known link function, \(\zeta\) is stochastic noise, and \(\beta^\star\in\mathbb{R}^d\) is an unknown index vector to be recovered. The recovery of \(\beta^\star\) is a fundamental problem in semi-parametric statistics [1]–[6] and machine learning [7]. In this paper, we focus on Gaussian designs \(X\sim\mathcal{N}(0,\mathbf{I}_d)\), the canonical setting for studying the geometry of non-convex SIM loss landscapes [7]–[14].

While classical SIM recovery techniques relied on clean-data assumptions [15], [16], in practice, data is invariably subject to noise and corruption. The field of robust statistics [17]–[19] developed estimators with well-understood robustness guarantees such as high breakdown points, primarily in low-dimensional settings. However, in high dimensions, achieving strong robustness often leads to estimators that are computationally intractable, with many classical formulations known to be NP-hard [20], [21]. Recent breakthroughs in algorithmic robust statistics [22], [23] overcame this computational barrier by providing efficient high-dimensional subroutines for robust mean estimation and robust PCA that tolerate both heavy-tailed noise and strong adversarial corruption. These subroutines serve as algorithmic building blocks for efficient robust recovery in structured models such as linear regression [24], [25], logistic regression [26], [27], and phase retrieval [28]–[30].

For monotonic link functions, [31] gave the first efficient recovery algorithm for monotone Lipschitz SIMs in the clean setting via the Isotron algorithm, later extended to broader classes of SIMs by [32] and [33]. Under heavy-tailed noise and strong adversarial contamination, near-linear time algorithms with optimal sample complexity \(\tilde{O}(d)\) were obtained for linear regression, i.e., \(f(x)=x\),  [24], [25] with [34] later obtaining information-theoretically optimal error rates under the same guarantees. Beyond linear regression, monotonic link functions such as logistic regression have been studied in the context of Generalized Linear Models (GLMs) [35]; for robust recovery under heavy-tailed noise and strong adversarial contamination, [26] proposed a polynomial-time, optimal-sample-complexity algorithm, while [27] later obtained near-linear time at the cost of polynomial sample complexity in the streaming setup. For non-monotonic link functions, phase retrieval, i.e., \(f(z) = z^2\), stands out as the canonical example, with applications in optics, crystallography, X-ray imaging, and astrophysics. [36] and [37] gave the first efficient recovery algorithms under Gaussian covariates in the absence of noise and corruption. Under strong adversarial contamination without noise (\(\zeta = 0\)), [28] proposed a near-linear time algorithm with optimal sample complexity \(\tilde{O}(d)\). [29] extended these results to heavy-tailed noise with an exponential-time algorithm, which was recently subsumed by [30], who gave a polynomial-time optimal sample complexity robust recovery algorithm for the same setting.

We note that the lack of robust recovery results for generic non-monotonic SIMs beyond phase retrieval is telling: standard first-order proof techniques  [11], [38], [39] break down under strong adversarial contamination in high dimensions, since the adversary destroys the statistical structure that gradient-based convergence arguments rely upon. The special case of phase retrieval enjoys two structural properties that are conducive for robust recovery: (i) its squared-loss landscape admits a convex basin of constant radius around the true parameter, enabling second-order convergence guarantees within the basin, and (ii) its non-vanishing second Hermite coefficient¹ ensures the signal direction is reachable via spectral methods  [36], [37]. It is a priori unclear whether these properties extend to generic non-monotonic SIMs, which often possess a high information exponent² (IE). Even in the clean setting, assuming one existed, accessing a convex basin for generic SIMs would require computationally demanding methods such as Tensor PCA [40], in stark contrast to the simple spectral methods that suffice for phase retrieval. Indeed, it remains unclear for which classes of non-monotonic link functions efficient provable robust recovery is achievable. This uncertainty naturally brings us to the following question:

Our Contributions: We make significant progress on the above question by identifying two general structural conditions on the link function \(f\) outlined in 1 2, under which efficient provable robust recovery is achievable. These conditions are remarkably general, capturing a large class of SIMs with generative exponent³ at most \(2\). This class includes phase retrieval, Tanh, Probit, and Logistic, as well as modern activation functions such as GeLU [41], Swish [42] which arise naturally as scalar primitives in gated neural architectures [43] that serve as the fundamental building blocks of Transformer architectures (e.g., GPT and BERT [44], [45]). We now state an informal characterization of our main result below.

1.1 Technical Overview↩︎

A natural strategy for robust recovery in generic non-monotonic SIMs is to leverage the geometry of the squared-loss landscape. Global convexity of the loss landscape would guarantee that any optimizer avoids spurious local minima, but for generic non-monotonic SIMs, global convexity of the squared loss holds if and only if \(f\) is affine, reducing to linear regression. For non-affine link functions, a straightforward approach is therefore to establish the existence of a convex neighborhood around the true parameter \(\beta^\star\). Once an optimizer enters this region, standard convex optimization guarantees ensure reliable recovery. However, for this strategy to be computationally viable in high dimensions, the convex basin must have radius \(R = O(1)\) strictly independent of the ambient dimension \(d\), since a vanishing basin of attraction offers no algorithmic guarantee in high-dimensional non-convex optimization. Prior to this work, such a dimension-independent convex basin was known to exist only for phase retrieval and monotonic link functions. This motivates our first contribution.

Establishing the existence of a convex basin does not by itself yield a computational guarantee. One must also certify that the basin can be reached efficiently from an initialization, even under heavy-tailed noise and strong adversarial contamination. The recent works of [29] and [30] make progress in this direction for phase retrieval, but their approaches suffer from two limitations. First, the robust PCA subroutines they employ run in time polynomial in the dimension. Second, and more fundamentally, their reachability arguments exploit the symmetric structure of the quadratic link and do not extend to asymmetric non-monotonic link functions such as GeLU and Swish. Beyond these computational limitations, certifying reachability under strong adversarial contamination in high dimensions poses a deeper analytical challenge. Prior proof techniques [11], [38], [39] that implicitly leverage convex basin structure rely on martingale-drift decompositions that require the stochastic deviations to be mean-zero, a property the adversary destroys by corrupting a constant fraction of samples in high dimensions. Our second contribution resolves both the computational and analytical obstacles.

2 demonstrates that the shortcomings of the previous approaches [29], [30] are not intrinsic. By coupling second-order Stein identities with a refined analysis of robust spectral estimation, we generalize the reachability insights of phase retrieval to a broad class of non-monotonic link functions while achieving near-linear initialization time. Under the Gaussian design, second-order Stein’s identity allows us to decompose the population moment matrix into a rank-one component aligned with \(\beta^\star\) and an isotropic component. Crucially, we show that the coefficient of the signal direction in this decomposition is governed precisely by the ESC condition, which allows us to identify the moment matrix whose leading eigenvector is \(\beta^\star\). Using this structural insight, we apply the off-the-shelf robust PCA subroutine of  [46], which explicitly leverages the hypercontractivity of the Gaussian design to extract the leading eigenvector of the sample moment matrices efficiently under adversarial contamination. This contrasts with earlier analyses  [29], [30], [47], which either fell short of this general structural characterization or did not exploit such concentration properties to ensure both robustness and efficiency.

However, attaining better error rates requires going beyond spectral initialization. In classical phase retrieval algorithms like Wirtinger Flow [36], spectral methods are used primarily to initialize the optimization within a basin of attraction, after which Gradient Descent (GD) is employed to converge to the exact solution. In the robust setting, this two-stage architecture is equally critical but for a fundamental statistical reason: we observe that relying solely on robust spectral estimators encounters a fundamental error floor of \(O(\epsilon^{1/4})\). To break this barrier, we use an off-the-shelf robust Gradient Descent (GD) subroutine that refines the initial estimate significantly while maintaining near-linear time complexity.

1.1.1 Related Work↩︎

To put our contributions in context, prior to this work, it was not known if efficient robust recovery in under both heavy-tailed noise and strong adversarial contamination was possible for generic link functions. We now give a brief overview of known results for robust recovery in SIMs beyond linear regression, and defer a detailed literature review to the appendix (see 6).

[26] proposed a polynomial-time \(\tilde{O}(d)\) sample recovery algorithm (SEVER) for logistic regression under strong adversarial contamination with respect to the hinge loss and the logistic loss. With respect to the squared loss, their sample complexity blows up to \(\tilde{O}(d^5)\). [27] obtain a near-linear time robust recovery streaming algorithm for logistic regression that requires \(\tilde{O}(d^2)\) samples, with respect to the squared loss. We note that both [26] and [27] obtain an error of \(\sigma\sqrt{\epsilon}\), similar to us. [35] obtain optimal sample complexity robust recovery guarantees under strong adversarial contamination for GLMs (which rely on only monotonic link functions) without any guarantees on the running time of their approach, and only under Gaussian noise, in contrast to our near-linear time robust recovery algorithm under heavy-tailed noise and strong adversarial contamination. However, our error rate of \(O(\sigma\sqrt{\epsilon})\) is worse than their \(O(\sigma\epsilon\log\frac{1}{\epsilon})\) error rate, which is optimal. The caveat, however, is that our error rate also holds for a large class of non-monotonic link functions, which are not addressed by [35]. We also remark that our error rate matches the best known error rate for existing robust recovery algorithms for non-monotonic link functions [26], [27], [29], [30] under heavy-tailed noise and strong adversarial contamination.

1.2 Organization of the paper↩︎

We detail our problem setup, model, and technical tools in 2. In 3, we prove the existence of a convex basin. Next, we obtain a linear-time optimal-sample complexity robust recovery algorithm in 4. Finally in 5, we discuss our contributions and outline multiple avenues for future directions.

2 Preliminaries↩︎

Notation For any positive integer \(n\), let \([n]\) denote the set \(\{1, 2, \ldots, n\}\). For a vector \(v \in \mathbb{R}^d\), \(\|v\|_2\) denotes its Euclidean norm. For two vectors \(u, v \in \mathbb{R}^d\), \(\langle u, v \rangle\) denotes their inner product. For a matrix \(M\), \(\|M\|_{\mathrm{op}}\) denotes its operator norm, \(\operatorname{Tr}(M)\) its trace, and we write \(M \succeq 0\) to indicate that \(M\) is PSD. We denote the identity matrix in \(d\) dimensions by \(\mathbf{I}_d\). The indicator function of an event \(E\) is denoted by \(\mathbb{I}(E)\). \(\mathcal{N}(\mu, \Sigma)\) denotes a multivariate Gaussian with mean \(\mu\) and covariance \(\Sigma\). \(\mathbb{E}[\cdot]\) denotes expectation. \(\widetilde{O}(\cdot)\) hides logarithmic factors in the dimension \(d\), sample size \(n\), corruption level \(\epsilon\), and failure probability \(\delta\). Finally, \(f'\) and \(f''\) denote the first and second derivatives of the link function \(f\), respectively.

Problem Setup We consider the single index model (SIM) defined in 1 , with covariates \(X \in \mathbb{R}^d\) and a known⁴ link function \(f\). In this section, we formally introduce the problem our paper tackles. We first define the notion of strong adversarial contamination, and then introduce our model.

We now state the natural structural assumptions we require on the link function. The first assumption controls the smoothness of the population loss landscape through one-dimensional expectations of \(f\) and its derivatives under the Gaussian measure.

In 1, the local Lipschitz complexity measure \(C_{\mathrm{lip}}(R)\) measures how quickly the curvature of the population loss can deteriorate as \(\beta\) moves away from \(\beta^\star\). A simple upper bound on \(C_{\mathrm{lip}}(R)\) (via Cauchy–Schwarz) is \(5(M_2 M_3 + M_1 M_4)\), where \[M_k := \sup_{\|\beta-\beta^\star\|\leq R} \left(\underset{Z\,\sim\,\mathcal{N}(0,\|\beta\|_2^2)}{\mathbb{E}} \left[f^{(k-1)}(Z)^4\right]\right)^{\frac{1}{4}}, k \in [4],\] where \(f^{(0)} := f\), \(f^{(1)} := f'\), \(f^{(2)} := f''\), \(f^{(3)} := f'''\). Crucially, since \(g\) is evaluated pointwise on \(Z \sim \mathcal{N}(0, \|\beta\|_2^2)\), the quantities \(M_1, M_2, M_3, M_4\) (and hence \(C_{\mathrm{lip}}(R)\)) are determined entirely by one-dimensional Gaussian integrals of \(f\) and its derivatives. Therefore, any link function satisfying 1 has a dimension-independent basin radius \(R\).

The second assumption complements the first by addressing a different question: while 1 guarantees that the loss landscape admits a dimension-independent convex basin near \(\beta^\star\), it does not guarantee (efficient) reachability to this convex basin from a random initialization. 2, termed Expected Squared Convexity (ESC), precisely characterizes when the signal direction \(\beta^\star\) is identifiable from the second-order moments of the data, enabling a robust spectral initialization that lands within the convex basin.

2.1 Tools↩︎

In this subsection, we present the main tools used in our work.

The univariate version of [lem:steins] states that for a Gaussian variable \(z \sim \mathcal{N}(0,1)\), \(\mathbb{E}[g(z)(z^2 - 1)] = \mathbb{E}[g''(z)]\). We now present the main results from algorithmic robust statistics that our paper builds upon.

We use a robust PCA subroutine of [46], for which we recall two definitions. For \(M \in \mathbb{S}^{d \times d}_{\succeq 0}\) and \(\rho \in [0,1]\), a unit vector \(u \in \mathbb{R}^d\) is a \(\rho\)-approximate energy \(1\)-PCA of \(M\) [46] if \(\langle u u^\top, M \rangle \ge (1-\rho)\,\|M\|_1\), where \(\|M\|_1 := \max_{\|v\|_2=1} \langle vv^\top, M\rangle\). A random vector \(X \in \mathbb{R}^d\) is \((p, C_p)\)-hypercontractive [46] if for all \(u \in \mathbb{R}^d\), \(\mathbb{E}[\langle u, X\rangle^p]^{1/p} \le C_p(\mathbb{E}[\langle u, X\rangle^2])^{1/2}\), for some constant \(C_p\) and even integer \(p\). We now state the linear-time robust PCA guarantee [46].

3 Existence of Convex Basins in the Loss Landscape of SIMs↩︎

Consider the population loss \(\mathcal{L}(\beta):= \frac{1}{2}\mathbb{E}[(f(X^{\top}\beta)-Y)^{2}]\) with the corresponding Hessian \(H(\beta):=\mathbb{E}[\left(f'(X_i^{\top}\beta))^{2} + (f(X_i^{\top}\beta) - Y_i) f''(X_i^{\top}\beta) \right) X_i X_i^{\top}].\) We now give a sufficient condition to characterize a class of link functions that admit a dimension-independent constant-sized convex basin around \(\beta^\star\) in the loss landscape, in the following theorem:

[lem:local95geometry] establishes that whenever the link function satisfies 1, the population loss landscape is \(\frac{\mu}{2}\)-strongly convex and \(\frac{\mu + 2\mu_{1}}{2}\)-smooth throughout \(\mathcal{B}(\beta^\star,R)\). Crucially, both the strong convexity constant \(\mu/2\) and the basin radius \(R\) are dimension-independent: they are determined entirely by one-dimensional integrals of \(f\) and its derivatives against the standard Gaussian. This dimension-independence is the key structural fact that enables our robust recovery guarantees in 4, since it allows the curvature of the loss to dominate the adversarial bias uniformly over the entire ball \(\mathcal{B}(\beta^\star,R)\), regardless of the ambient dimension \(d\). We present an extended proof sketch below and defer the full proof to 7 in the appendix.

4 Linear-Time Robust Recovery of SIMs↩︎

In this section, we present our linear-sample and time algorithm for robustly recovering the true signal \(\beta^\star\) when the link function admits a convex basin in the loss landscape under heavy-tailed noise and strong adversarial contamination. We now formally state our main result.

In the above 2, \(m\) is the sample complexity for the robust spectral initialization subroutine (see 2) and \(\tilde{m}\) is the sample complexity for the robust gradient descent subroutine (see 3). We begin by giving the pseudocode of [alg:lineartimealgo] in 2 along with a high-level overview.

First, we begin by recalling that if our link function \(f\) satisfies 1 2, the population loss-landscape admits a convex basin in a neighborhood of \(\beta^\star\) (see [lem:local95geometry]). [alg:lineartimealgo] exploits this structure via the LRSI subroutine (see 2) which combines generalized higher-order Stein’s identities (see [lem:steins]) together with the linear-time robust hypercontractive \(1\)-ePCA algorithm as described in 2 to construct an estimate \(\beta_0\) of the true signal \(\beta^\star\) that lies within the convex basin. Finally, using \(\beta_0\) as a warm start, the LRGD algorithm (see 3) converges to \(\beta^\star\). We first carefully detail the robust spectral initialization step, then the robust gradient descent step, and finally the proof of 2.

4.1 Warm Start via Linear time Spectral Initialization↩︎

[lem:spectral] and 2 state that for the random variable \(\widetilde{Y} = YX\), where \(Y,X\) are defined in 2 s.t. the link function \(f\) satisfies 2, the vector \(\beta^\star\) is the leading eigenvector of the matrix \(\mathbb{E}[\widetilde{Y}\widetilde{Y}^\top]\), with eigenvalue \(\lambda_{\max}\). For a warm start, we robustly estimate the leading eigenvector of the second-moment matrix of \(\widetilde{Y}\). To this end, we employ the robust hypercontractive \(1\)-ePCA algorithm. The applicability of this method requires the distribution of \(\widetilde{Y}\) to be hypercontractive, a property we establish in [lem:hypercontracitivity32of32YX].

We now state the guarantees for the LRSI algorithm (see 2) that uses the PCA algorithm of [46] as a subroutine.

4.2 Linear Robust Gradient Descent↩︎

The second key step, following phase retrieval, is to perform gradient descent on the population risk, \(\beta_{t+1} = \beta_t - \eta \nabla \mathcal{L}(\beta_t).\) By Theorem 3.12 of [51] together with Lemma [lem:local95geometry], the iteration stated above converges linearly to the global minimizer, provided that all iterates remain within the ball \(\mathcal{B}(\beta^\star,R)\) and the step size is chosen as \(\eta = \frac{2}{\alpha+\gamma}\), where \(\alpha\) and \(\gamma\) denotes the smoothness and the strong convexity parameters, respectively. [lem:local95geometry] implies that, for the population loss \(\mathcal{L}(\beta)\), the strong convexity parameter is \(\gamma =\frac{\mu}{2}\) and the smoothness parameter is \(\alpha = \frac{\mu}{2} + \mu_1\). Since the learning algorithm only has access to the dataset and not the population gradient \(\nabla \mathcal{L}(\beta)\), the update stated above cannot be implemented directly. To address this, we follow the standard approach of expressing the gradient as an expectation [25], [29]. In particular, \(\nabla \mathcal{L}(\beta) = \mathbb{E}\!\left[(f(X^\top \beta)-Y)f'(X^\top \beta)X\right].\) Substituting the value of \(\nabla \mathcal{L}(\beta)\) into above update yields \[\beta_{t+1} = \beta_t - \eta\,\mathbb{E}\!\left[(f(X^\top \beta_t)-Y)f'(X^\top \beta_t)X\right].\] We then replace the expectation with a robust estimator of the gradient. We state the definition of a robust gradient estimator (see Definition 2.1.1 in [29]) below.

Following [29], we use the notion of a robust gradient estimator, stated formally in Definition 4. Let \(g_t := g(\beta_t; T, \delta, \epsilon)\) denote such an estimator computed from the dataset \(T\). The resulting robust gradient descent update is \(\beta_{t+1} = \beta_t - \eta g_t\). We summarize the resulting robust gradient descent procedure in Algorithm 3, and then state its main guarantees.

Proof of 2. Under the assumptions on \(m\) and \(\epsilon\), the output \(\beta_{0}\) of the LRSI algorithm (see 2) satisfies \(\|\beta_{0} - \beta^{\ast}\| = O\!\left( \frac{C_{4}\bigl(\sigma^{2} + \mathbb{E}[f^{2}] + c\bigr)^{1/2}\epsilon^{1/4}}{\sqrt{c}} \, \right).\) Moreover, under the additional assumption on the corruption level \(\epsilon \leq \frac{R^{4} c^{2}}{C_{4}^{4}\bigl(\sigma^{2} + \mathbb{E}[f^{2}] + c\bigr)^{2}},\) we have \(\|\beta_{0} - \beta^{\ast}\| \leq R\). The proof now follows directly from [theorem:32li-time-grads-descent].

4.3 Applications↩︎

In this section, we now demonstrate near-linear time and sample robust recovery for SIMs with \(6\) different link functions under heavy-tailed noise and strong adversarial contamination. Our representative link functions can be broadly classified into three categories: (i) Monotonic Links: Logistic/Sigmoid (\(\sigma(z)\)), Tanh (\(\tanh(z)\)), and Probit (\(\Phi(z)\)). Here, \(\phi(z) = \frac{e^{-z^2/2}}{\sqrt{2\pi}}\) and \(\sigma(z) = \frac{1}{1+e^{-z}}\), (ii) Phase Retrieval (\(f(z) = z^2\)), and (iii) Asymmetric Non-monotonic Links: GeLU (\(z\Phi(z)\)), and Swish (\(z\sigma(z)\)).

5 Discussion and Future Work↩︎

In this work, we established the first framework achieving near-linear time and optimal sample complexity for the robust recovery of Single-Index Models with generic, non-monotonic link functions (e.g., GeLU, Swish) under heavy-tailed noise and strong adversarial contamination. Below we outline a few interesting future directions.

Optimal Error Rates. Our estimator achieves an \(\ell_2\) error rate of \(O(\sigma\sqrt{\epsilon})\) under Gaussian covariates in the presence of heavy-tailed noise and \(\epsilon\)-fraction adversarial contamination. In comparison, for the same setting, [24], [34] established the information-theoretically optimal rate \(\tilde{O}(\sigma\epsilon)\). However, for non-linear models with non-convex population loss, existing provably robust algorithms (e.g., phase retrieval), are only known to achieve a \(O(\sigma\sqrt{\epsilon})\) rate [29], [30]. Our result matches this best-known rate for non-linear single-index models while accommodating a significantly broader class of link functions, including non-monotonic ones. Closing the gap between the achievable rate and the information-theoretically optimal rate for general SIMs remains an important open problem and we leave it for future work.

Non-Gaussian Covariates. Our theoretical guarantees heavily leverage the Gaussianity of the design matrix to derive the ESC condition, and in our basin radius analysis by exploiting rotational invariance. We leave extending our proofs to even sub-Gaussian designs as an open question.

Alternative Loss Landscapes and Adversary Models. While we focused on the squared loss, investigating if convex basins persist under Huber loss or general \(M\)-estimators remains an important open question. Additionally, since the existence of a convex-basin does not depend on the corruption model, adapting our framework to Agnostic Learning or Differential Privacy settings is a promising future direction.

Multi-Index Models (MIMs) As noted in our introduction, functions like GeLU and Swish are scalar primitives for GLUs, which are inherently Multi-Index Models (MIMs) defined by interactions between multiple projections, as \(y = \langle \beta_1, x \rangle \cdot f(\langle \beta_2, x \rangle)\). Extending the guarantees of our work to MIMs (particularly 1 2) requires disentangling the interaction terms between multiple weight vectors. This likely necessitates robust tensor decomposition techniques of order significantly higher than those required for SIMs, which presents a distinct set of algebraic and algorithmic challenges.

Robust Recovery for Links with Information Exponent \((k \geq 3)\). Our framework primarily targets link functions where the signal is detectable via low-order derivatives (specifically, where ESC is non-trivial). However, for link functions with an information exponent \(k \geq 3\), the signal is entirely suppressed in lower-order moments. Robustly recovering \(\beta^\star\) in this regime would require working with higher-order moment tensors (of order at least \(k\)). This is an interesting avenue for future work and one concrete line of investigation is discussed next.

Identifying Label Transforms. Suppose there exists a map \(\tau:\mathbb{R}\to\mathbb{R}\) such that \(\tilde{f}=\tau\circ f\) has IE \(k^\star\leq 2\) and satisfies 1 2. Then our framework applies directly to \(\tilde{f}\), yielding efficient robust recovery for the original link function \(f\). Identifying such transforms is highly non-trivial, since \(k^\star\) is the infimum of the IE over all square-integrable label transforms [13], and characterizing which link functions admit a transform \(\tau\) such that the resulting \(\tilde{f}\) has IE \(\leq 2\) and satisfies 1 2 is an open problem. We view this as a promising direction for future work, and note that our results provide the first motivation for investigating such regularity conditions on label transforms.

Empirical Verification. While our focus is theoretical, empirical evaluation is a vital next step. Implementing high-order robust spectral estimators involves practical engineering challenges, particularly regarding numerical stability and hyperparameter tuning for the filtering subroutines. We leave the extensive experimental benchmarking of these algorithms on real-world datasets, and the potential development of practically optimized heuristics based on our theory, for future work.

Acknowledgements↩︎

The authors are grateful to Ankit Pensia, for many valuable discussions and in particular for pointing us to the robust PCA and robust mean estimation subroutines used in 4 which was instrumental in obtaining the near-linear time guarantee of 2. The authors also thank the anonymous reviewers of ICML 2026 for their thorough and constructive feedback, which helped improve the presentation of this work. This work was supported by the Department of Atomic Energy, Government of India, under project no. RTI4014.

6 Literature Review↩︎

Monotone link functions without corruption and heavy-tailed noise: [31] study the recovery of monotone Lipschitz single-index models (SIMs) with an unknown link function \(f\) via the Isotron algorithm, which operates in \(\mathrm{poly}(d)\) time with \(O(1)\) sample complexity. For arbitrary link functions, [52] provide a linear-time estimator requiring \(O(d)\) samples for consistent recovery. Further extending this scope, [12] analyze link functions with polynomial-tailed derivatives using online stochastic gradient descent on a smoothed loss function. Their approach achieves a sample complexity of \(n = O(d + d^{k^*/2})\) and time complexity of \(O(nd)\), where \(k^*\) denotes the information exponent-the smallest nonzero coefficient of the link function’s Hermite expansion under the Gaussian measure.

Monotone link function under adversarial corruption and heavy-tailed: [55] study linear and polynomial regression with Gaussian covariates, and propose a Sum-of-Squares–based algorithm achieving polynomial sample and computational complexity in \(d\). [26] study robust learning in two settings: for SIMs with a monotone (logistic) link, they give a polynomial-time algorithm with sample complexity \(\tilde{O}(d)\) achieving generalization error \(O(\epsilon^{1/4})\) under logistic loss; for linear regression with heavy-tailed noise, they propose a polynomial-time algorithm with sample complexity \(O(d^5)\) achieving estimation error \(O(\sqrt{\epsilon})\), where \(\epsilon\) denotes the contamination level. [23] study linear regression under Gaussian noise, and propose a polynomial-time algorithm achieving sample complexity \(\tilde{O}(d)\) and estimation error optimal up to constant factors. [35] study monotone SIMs under adversarial contamination and propose a algorithm based on trimmed MLE, achieving sample complexity \(\tilde{O}(d)\) with heuristic time complexity.

Phase retrieval: Phase retrieval is a classical inverse problem with applications in optics, crystallography, and X-ray imaging. In the non-contaminated, light-tailed setting, the first provable algorithm was due to [37], followed by a line of work including [54], [57] and others. [37] propose an alternating minimization approach, while [36] use a gradient-based method; both achieve sample complexity \(\widetilde{O}(d)\) and runtime \(\widetilde{O}(nd)\). Notably, the method of [37] requires fresh samples at each iteration, whereas [36] does not. Whereas, [54] study sparse phase retrieval (\(\beta^\star\) is \(s\) sparse) and propose an algorithm, achieves sample complexity \(\tilde{O}(s^2)\) and computational complexity \(\tilde{O}(s^2n)\). [53] study phase retrieval under dense bounded noise and label corruption using a median-truncated Wirtinger Flow algorithm with sample complexity \(\widetilde{O}(d)\) and runtime \(\widetilde{O}(nd)\), while [56] consider same settings as [53] and propose an SDP-based algorithm achieving sample complexity \(O(d)\). Phase retrieval under strong adversarial contamination of both labels and covariates has been studied in [28]–[30]. Specifically, [28] consider the noiseless setting and propose a linear-time algorithm with sample complexity \(\widetilde{O}(d)\), while [29] study heavy-tailed label noise and give an exponential-time algorithm achieving sample complexity \(\widetilde{O}(d)\). Subsequently, [30] obtain the same optimal sample complexity \(n=\tilde{d}\) with runtime \(\widetilde{O}(n^{2}d)\).

7 Details and Omitted Proofs of 3↩︎

7.1 Gaussian symmetry↩︎

In this section, we state and prove an important property of the Gaussian distribution, which will be used later in the proof of [lem:local95geometry].

7.2 Proof of [lem:local95geometry]↩︎

8 Sample splitting in [alg:lineartimealgo]↩︎

9 Omitted Proofs of 4↩︎

First, we present two key lemmas needed for the proof of [theorem:32li-time-grads-descent], followed by the proof itself. The first key lemma expresses the distance of the iterate at time \(t+1\) from \(\beta^\star\) in terms of the distance of the iterate at time \(t\) from \(\beta^\star\). The relation is as follows:

The above lemma and its proof are standard in robust statistics literature; for example, see Theorem 1 in [25] and Lemma 2.2 from [29]. Look at any of the above-mentioned references for the proof of the 4. We are omitting the proof here. Our second key lemma computes the trace and operator norm of the covariance matrix of the gradient.

Now we go for the proof for [theorem:32li-time-grads-descent]. The proof is an adaptation of the proof of the guarantees of the robust gradient descent algorithm (particularly theorem 3.3) from [29]. The calculations will vary, but the overall flow of arguments remains the same.

10 Numericals↩︎

10.1 Minimum and maximum eigenvalue of the expected Hessian at true signal↩︎

From the proof of [lem:local95geometry], we know that \[\mu = \lambda_{\min}(H(\beta^\star)) = \lambda_{\min}\left( \mathbb{E}_{x}\left[ (f'(x^\top\beta^\star))^2 xx^\top \right] \right),\quad \mu_{1} = \lambda_{\max}(H(\beta^\star)) = \lambda_{\max}\left( \mathbb{E}_{x}\left[ (f'(x^\top\beta^\star))^2 xx^\top \right] \right)\] Let \(z^* = x^\top\beta^\star\). Since \(x \sim \mathcal{N}(0, \mathbf{I}_d)\), then \(z^*\sim \mathcal{N}(0, 1)\) because \(\|\beta^\star\|=1\). Due to the rotational symmetry of the Gaussian distribution, the matrix \(H(\beta^\star)\) has only two distinct eigenvalues:

1. \(\lambda_{\parallel}\): Associated with the eigenvector aligned with \(\beta^\star\).

2. \(\lambda_{\perp}\): Associated with eigenvectors orthogonal to \(\beta^\star\) (with multiplicity \(d-1\)).

These are computed as: \[\begin{align} \lambda_{\perp} &= \mathbb{E}_{z^*}\left[ (f'(z^*))^2 \right] \\ \lambda_{\parallel} &= \mathbb{E}_{z^*}\left[ (f'(z^*))^2 (z^*)^2 \right] \end{align}\] We compute these explicitly for each link function. Thus, \(\mu=\min\{\lambda_{\perp},\lambda_{\parallel}\}\), and \(\mu_{1}=\max\{\lambda_{\perp},\lambda_{\parallel}\}\).

10.2 Derivation of \(C_{lip}\) and \(R\)↩︎

The condition for convexity relies on the bound: \[8.43 \cdot C_{lip} \cdot R < \mu\] where \(C_{lip}\) is defined as: \[C_{\mathrm{lip}}(R) := \sup_{\|\beta-\beta^\star\| \le R} \sqrt{\mathbb{E}_{z\sim \mathcal{N}(0,\|\beta\|^2)}[g(z)]},\] and \(g(z) := 18 f'(z)^2f''(z)^2 + 2f'''(z)^2f(z)^2\).

1. For Bounded Derivatives (Logistic, Tanh, Probit): From the proof of [lem:local95geometry], we have \[\sqrt{\mathbb{E}\!\left[A(z,z^*)^2\right]} \le C_{\mathrm{lip}}(R),\] where \[A(z,z^*) = 3 f'(z')f''(z') + 2 f'''(z')f(z'),\] and \(z'\) is a point on the line segment joining \(z\) and \(z^*\). For the case of bounded derivatives, we bound \(C_{\mathrm{lip}}(R)\) using a uniform global bound on \(|A(z,z^*)|\), namely, \[C_{\mathrm{lip}}(R) := \sup_{z,\,z^*} |A(z,z^*)|.\] Then \(R\leq \frac{\mu}{8.43 C_{\mathrm{lip}}(R)}.\)

2. For Unbounded Derivatives (Phase Retrieval, GeLU, Swish, SwiGLU): For the case of unbounded derivatives, we proceed differently. We first compute a function \(g(z)\) such that and then express \[\mathbb{E}_{z \sim \mathcal{N}(0,\|\beta\|^2)}[g(z)]\] explicitly as a function of \(\|\beta\|\). Substituting \(\|\beta\| \le R+1\) yields an upper bound of the form \(\chi(R)\), where \(\chi(\cdot)\) is a deterministic function of \(R\).

   We then choose \(R\) so as to satisfy \[8.43\, R\, \chi(R) \le \mu,\] using any suitable numerical or analytical method. With this choice, \(\chi(R)\) serves as \(C_{\mathrm{lip}}(R)\), and the corresponding \(R\) defines the radius of the convex basin.

10.3 Derivations of \(c,\phi_{1}\) and \(\phi_{2}\).↩︎

We know that \(c=\mathrm{ESC}(\beta; f) :=\mathbb{E}_{Z\sim\mathcal{N}(0,1)}\!\left[ \bigl(f'(Z)\bigr)^2 + f(Z)\, f''(Z) \right]\). From [theorem:32li-time-grads-descent], we know that \(\phi_{1}=\sup_{\beta\in \mathcal{B}(\beta^\star, R)}\mathbb{E}\!\left[f'(X^{\top}\beta)^{16})\right]^{1/4}\) and \(\phi_{2}=\sup_{\beta\in \mathcal{B}(\beta^\star, R)}\mathbb{E}\!\left[f'(X^{\top}\beta)^{4}\right]^{1/2}\).

1. For Bounded Derivatives (Logistic, Tanh, Probit): For the case of bounded derivatives, we obtain a uniform upper bound on \(f'\), which can then be used to compute the constants \(\phi_{1}\) and \(\phi_{2}\).

2. For Unbounded Derivatives (Phase Retrieval, GeLU, Swish, SwiGLU): Note that \[\begin{align} \phi_{1}&=\sup_{\beta\in \mathcal{B}(\beta^\star, R)}\mathbb{E}\!\left[f'(X^{\top}\beta)^{16})\right]^{1/4}=\sup_{\beta\in \mathcal{B}(\beta^\star, R)}\mathbb{E}_{Z\sim\mathcal{N}(0,1)}\!\left[f'(\|\beta\|Z)^{16})\right]^{1/4}\\ &=\sup_{\|\beta\| \in [1-R, 1+R]}\left( \mathbb{E}\left[f'(\|\beta\|Z)^{16}\right] \right)^{1/4}=\sup_{s \in [1-R, 1+R]}\left( \mathbb{E}\left[f'(sZ)^{16}\right] \right)^{1/4}. \end{align}\] Similarly \(\phi_{2}=\sup_{s \in [1-R, 1+R]}\left( \mathbb{E}\left[f'(sZ)^{4}\right] \right)^{1/2}\). So, once \(R\) is determined, we can use the above forms for \(\phi_{1}\) and \(\phi_{2}\) to compute their values, either via numerical methods or with the aid of standard software packages. In certain cases, such as phase retrieval, these quantities can also be computed in closed form.

10.4 Code for Numerical Calculations↩︎

# -*- coding: utf-8 -*-
"""Calculations
"""

import numpy as np
# Enable JAX 64-bit precision BEFORE importing other JAX modules
from jax import config
config.update("jax_enable_x64", True)

import jax.numpy as jnp
from jax import grad, jit, vmap
from jax.scipy.stats import norm as jax_norm
from scipy.integrate import quad
from scipy.optimize import minimize_scalar, brentq
import pandas as pd

# ==========================================
# 1. Definitions of the Functions
# ==========================================

def sigmoid(z):
    return 1.0 / (1.0 + jnp.exp(-z))

def probit(z):
    return jax_norm.cdf(z)

def phi_pdf(z):
    return jax_norm.pdf(z)

# Function dictionary mapping names to their JAX implementations
functions = {
    "Logistic/Sigmoid": lambda z: sigmoid(z),
    "Tanh": lambda z: jnp.tanh(z),
    "Probit": lambda z: probit(z),
    "Phase Retrieval": lambda z: z**2,
    "GeLU": lambda z: z * probit(z),
    "Swish": lambda z: z * sigmoid(z),
    "GeGLU": lambda z: z**2 * probit(z),
    "SwiGLU": lambda z: z**2 * sigmoid(z)
}

# ==========================================
# 2. Derivative & Expectation Helpers
# ==========================================

class FunctionAnalyzer:
    def __init__(self, func_name, func_jax):
        self.name = func_name
        self.f = func_jax

        # JIT compile derivatives for speed
        self.d1 = jit(grad(self.f))
        self.d2 = jit(grad(self.d1))
        self.d3 = jit(grad(self.d2))

    def get_derivatives(self, z):
        """Returns f(z), f'(z), f''(z), f'''(z)"""
        return self.f(z), self.d1(z), self.d2(z), self.d3(z)

    def expected_value(self, integrand_fn, sigma=1.0):
        """
        Computes E[integrand(Z)] where Z ~ N(0, sigma^2).
        We transform variables: Z = sigma * x where x ~ N(0, 1).
        """
        # pdf of standard normal
        norm_pdf = lambda x: (1 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * x**2)

        def wrapper(x):
            z_val = sigma * x
            # Get function terms at z_val
            f, df, d2f, d3f = self.get_derivatives(z_val)
            # Calculate specific integrand term
            val = integrand_fn(z_val, f, df, d2f, d3f)
            # Cast to float to ensure compatibility with scipy.quad
            return float(val * norm_pdf(x))

        # Integrate from -10 to 10.
        # Mass outside [-10, 10] is < 1e-23, which is effectively 0 for integration.
        # points=[0.0] ensures the integrator splits intervals at the peak.
        res, _ = quad(wrapper, -10.0, 10.0, points=[0.0], limit=100)
        return res

# ==========================================
# 3. Calculation Logic for Each Constant
# ==========================================

def calculate_all_constants(analyzer):
    print(f"--- Processing: {analyzer.name} ---")

    # --- 0. Calculate ESC (at sigma=1) ---
    # Definition: E[ f'(Z)^2 + f(Z)f''(Z) ]
    def term_ESC_fn(z, f, df, d2f, d3f):
        return df**2 + f * d2f
    val_ESC = analyzer.expected_value(term_ESC_fn, sigma=1.0)

    # --- 1. Calculate mu and mu1 (at sigma=1) ---
    # Term A: E[f'(Z)^2]
    def term_A_fn(z, f, df, d2f, d3f): return df**2
    val_A = analyzer.expected_value(term_A_fn, sigma=1.0)

    # Term B: E[Z^2 * f'(Z)^2]
    def term_B_fn(z, f, df, d2f, d3f): return (z**2) * (df**2)
    val_B = analyzer.expected_value(term_B_fn, sigma=1.0)

    mu = min(val_A, val_B)
    mu1 = max(val_A, val_B)

    # --- 2. Define C_lip(R) Calculator ---
    # C_lip is the sup over beta ball.
    # Ball ||beta - beta*|| <= R implies ||beta|| is in [max(0, 1-R), 1+R].
    # We simplify this to finding max expectation over sigma in this range.

    def get_Clip_integrand_at_sigma(sigma):
        # E[18 f'(Z)^2 f''(Z)^2 + 2 f'''(Z)^2 f(Z)^2]
        def inner(z, f, df, d2f, d3f):
            return 18 * (df**2) * (d2f**2) + 2 * (d3f**2) * (f**2)
        return analyzer.expected_value(inner, sigma=sigma)

    def solve_C_lip(R):
        # We want to maximize the expectation over sigma
        bounds = (max(0.001, 1 - R), 1 + R)

        # minimize_scalar finds minimum, so we minimize negative
        res = minimize_scalar(lambda s: -get_Clip_integrand_at_sigma(s),
                              bounds=bounds, method='bounded')
        return -res.fun # Return the maximum value found

    # --- 3. Solve for R ---
    # Equation: R = mu / (2 * (315^0.25) * C_lip(R))
    # Let g(R) = LHS - RHS. We want g(R) = 0.

    const_factor = 2 * (315**0.25)

    def equation_to_solve(R):
        if R <= 0: return -1.0 # R must be positive
        c_val = solve_C_lip(R)
        # Avoid division by zero if C_lip is 0
        if c_val < 1e-9: c_val = 1e-9
        return R - (mu / (const_factor * c_val))

    # Dynamic Bracketing to find R
    # We check a range. If signs are same, we expand the range.

    low, high = 1e-6, 5.0
    try:
        f_low = equation_to_solve(low)
        f_high = equation_to_solve(high)

        solved_R = 0.01 # Fallback

        if np.sign(f_low) == np.sign(f_high):
            if f_low > 0:
                # Both positive: R > RHS for low and high.
                # This Mean Value Theorems the root is very small (Left of low).
                # Try finding in [1e-9, low]
                try:
                    solved_R = brentq(equation_to_solve, 1e-9, low)
                except:
                    solved_R = 1e-9 # Cap at min
            else:
                # Both negative: R < RHS for low and high.
                # The root is to the right of high.
                # Try expanding high significantly.
                try:
                    solved_R = brentq(equation_to_solve, high, 100.0)
                except ValueError:
                    print(f"Warning: R > 100.0 for {analyzer.name}. Using 100.0")
                    solved_R = 100.0
        else:
            # Signs differ, standard solve
            solved_R = brentq(equation_to_solve, low, high)

    except Exception as e:
        print(f"Error solving R for {analyzer.name}: {e}. Defaulting to 0.01")
        solved_R = 0.01

    # Recalculate C_lip at the solved R for reporting
    final_Clip = solve_C_lip(solved_R)

    # --- 4. Solve for phi_1 and phi_2 ---
    # Both are supremums over the ball (sigma in [1-R, 1+R])

    # phi_1: sup (E[f'(Z)^16])^(1/4)
    def get_phi1_val(sigma):
        def inner(z, f, df, d2f, d3f): return df**16
        return analyzer.expected_value(inner, sigma)

    res_phi1 = minimize_scalar(lambda s: -get_phi1_val(s),
                               bounds=(max(0.001, 1-solved_R),
                               1+solved_R), method='bounded')
    phi1 = (-res_phi1.fun)**0.25

    # phi_2: sup (E[f'(Z)^4])^(1/2)
    def get_phi2_val(sigma):
        def inner(z, f, df, d2f, d3f): return df**4
        return analyzer.expected_value(inner, sigma)

    res_phi2 = minimize_scalar(lambda s: -get_phi2_val(s),
                               bounds=(max(0.001, 1-solved_R),
                               1+solved_R), method='bounded')
    phi2 = (-res_phi2.fun)**0.5

    return {
        "ESC": val_ESC,
        "mu": mu,
        "mu1": mu1,
        "R": solved_R,
        "C_lip(R)": final_Clip,
        "phi1": phi1,
        "phi2": phi2
    }

# ==========================================
# 4. Main Execution
# ==========================================

if __name__ == "__main__":
    results = []

    for name, func_jax in functions.items():
        analyzer = FunctionAnalyzer(name, func_jax)
        try:
            res = calculate_all_constants(analyzer)
            res["Function"] = name
            results.append(res)
        except Exception as e:
            import traceback
            traceback.print_exc()

    # Create DataFrame
    df = pd.DataFrame(results)

    # Reorder columns
    cols = ["Function", "ESC", "mu", "mu1", "R", "C_lip(R)", "phi1", "phi2"]
    if not df.empty:
        df = df[cols]

    print("\n\n=== Final Constants Table ===")
    # Format for nicer reading
    pd.set_option('display.float_format', lambda x: '    print(df.to_string(index=False))

References↩︎