A regularization-based framework combining a Huberized data fidelity with generalized folded-concave penalties (SCAD, MCP), and a two-block proximal alternating algorithm with backtracking (NLD-PALM) whose whole iterate sequence provably converges to critical points under the Kurdyka--\L{}ojasiewicz property, with local linear rates.
Abstract
We consider the recovery of a pair of sparse vectors from a limited number of nonlinear observations of their superposition: $y_i=g(\inner{\ba_i}{\bPhi\bw^\ast+\bPsi\bz^\ast})+e_i$, $i=1,\dots,m$, with $m\ll n$, incoherent orthonormal bases $\bPhi,\bPsi$, a scalar link $g$, and noise $e_i$ that may be heavy-tailed or contaminated. We propose a regularization-based framework combining a Huberized data fidelity with generalized folded-concave penalties (SCAD, MCP), and a two-block proximal alternating algorithm with backtracking (NLD-PALM) whose whole iterate sequence provably converges to critical points under the Kurdyka--\L{}ojasiewicz property, with local linear rates. On the statistical side we establish restricted strong convexity of the Huberized nonlinear loss through an exact sign-definite decomposition, and derive estimation error bounds of order $\sigma\sqrt{s\log(n)/m}$ that hold at \emph{every} localized stationary point, an oracle rate $\sigma\sqrt{s/m}$ free of $\log n$ and shrinkage bias under a beta-min condition, and a co-equal recovery theorem for \emph{unknown} monotone links via a linear surrogate and a clipped Plan--Vershynin decoupling. The estimator requires no knowledge of the sparsity levels, and its guarantees hold under symmetric noise with only finite variance. Experiments at $n=512$ under a frozen data-driven regularization rule show an earlier phase transition than convex $\ell_1$ demixing and greedy hard-thresholding baselines, a $35\times$ accuracy advantage over squared-loss estimation under $5\%$ gross outliers, and successful demixing of spike-plus-background signals observed through a saturating amplifier.
It is proved that membership in the approximation space $k_t$ is equivalent to polynomial decay of the best $n$-term approximation error, which is equivalent to polynomial decay of the best $n$-term approximation error.
Abhishake Rastogi, T. Bubba, T. Helin et al.· 0 citations
The problem of recovering an (approximately) low-rank Hermitian matrix $\pmb{M}_0 \in \mathbb{C}^{n \times n}$ of rank $r$ from quadratic sampling matrices of the form $\{\pmb{a}_k \pmb{a}_k^*\}_{k=1}^m$ arises in a variety of applications, including phase retrieval. To obtain rigorous recovery guarantees, the sampling vectors $\{\pmb{a}_k\}_{k=1}^m$ are typically modeled probabilistically. However, most existing theoretical results rely on Gaussian or sub-Gaussian assumptions, which may not accurately capture practical data models. In many applications, sampling vectors exhibit heavier tails, while theoretical understanding in such regimes remains scarce. In this paper, we bridge this gap. We show that two widely used convex approaches, nuclear norm minimization and semidefinite-constrained empirical risk minimization, achieve uniform, stable, and robust recovery under the mild assumption that the entries of the sampling vectors have only finite $4+\delta$ moments, with the optimal sample complexity $m = \mathcal{O}(rn)$ up to moment-dependent constants. The two main ingredients of our analysis are moment estimates for quadratic forms established via decoupling, together with recent advances in covariance estimation in heavy-tailed settings. As byproducts, we also establish the optimal sample complexity for low-rank matrix recovery under complex projective $4$-design sampling, thereby improving upon previous results, and obtain stability guarantees for phase retrieval under similarly weak moment assumptions.
A classical problem in sparse Fourier transforms, which dates back to the work by Prony in 1795 at least, is to learn a $k$-Fourier-sparse signal $x(t):=\sum_{j=1}^k \alpha_j e^{2 \pi \mathbf{i} f_j t}$ with arbitrary frequencies $f_1,\ldots,f_k$. We study this problem of learning $x(t)$ in a fixed time window $[-T,T]$ under adversarial noise with bounded $\ell_2$ norm, where the frequencies $f_1,\ldots,f_k$ may be"off-grid"-- arbitrarily located in a given bandlimit $[-F,F]$. In particular, our goal is to output a sparse interpolation $\tilde{x}$ such that $\tilde{x}(t) \approx x(t)$ in the time window $[-T,T]$. 1. Our first result shows that the sample complexity of interpolation is $k^2 \cdot O(\log \frac{k FT}{\epsilon})^2$. While its running time is $(\frac{k FT}{\epsilon})^{O(k)}$, this improves the previous upper bound $k^{4} \cdot (\log FT)^{O(1)}$ on the sample complexity substantially and leaves a gap of about $k$ to the lower bound $\Omega(k \log FT)$. 2. Our second result provides efficient algorithms to interpolate $x(t)$. The first algorithm takes $m=k^{3.75} \cdot (\log FT)^{O(1)}$ samples and $m^{\omega+o(1)}$ time ($\omega$ is the matrix multiplication exponent). Assuming that the growth of any $k$-Fourier-sparse signal cannot be significantly larger than the growth of the degree-$(k-1)$ Chebyshev polynomial -- specifically, $x(t) \le e^{k \cdot O\big( \sqrt{\frac{|t|}{T}-1} \big)} \cdot \underset{s \in [-1,1]}{\max} |x(s)|$ for any $t \notin [-T,T]$, the second algorithm further improves the sample complexity to $m'=k^{3} \cdot (\log FT)^{O(1)}$ and the time complexity to $(m')^{\omega+o(1)}$.
Dongrun Cai, Xue Chen, Xiaowei Shao et al.· 0 citations
Recovering a low-dimensional latent subspace from nonlinear observations of Gaussian covariates in high dimensions is a fundamental problem in feature learning. Here, we consider Gaussian multi-index models in which the covariates $\boldsymbol{x}_i \stackrel{\mathrm{i.i.d.}}{\sim} \mathcal{N}(0,\boldsymbol{I}_d)$ and the responses $\boldsymbol{y}_i$ depend on $\boldsymbol{x}_i$ only through its projection onto an unknown $r$-dimensional subspace. Earlier work based on approximate message passing (AMP) identified a sharp threshold for weak recovery [Troiani et al., 2025], raising the question of whether it can be attained, without side information, by a spectral method. We answer this affirmatively and develop a general random matrix theory for matrix-valued spectral estimators of the form \[\boldsymbol{D}_n=\frac{1}{n}\sum_{i=1}^n\boldsymbol{T}(\boldsymbol{y}_i)\otimes\boldsymbol{x}_i\boldsymbol{x}_i^\top,\] where $\boldsymbol{T}$ is an arbitrary bounded symmetric matrix-valued preprocessing map of fixed dimension. As $n,d \to \infty$ with $n/d\to\alpha$, we prove that the empirical spectral measure of $\boldsymbol{D}_n$ converges almost surely to a deterministic compactly supported distribution characterized by a matrix-valued self-consistent equation. We then establish a spectral phase transition for the largest eigenvalue: below threshold it sticks to the bulk edge, while above threshold an outlier emerges. We characterize the outlier location through a finite-dimensional deterministic equation and show that the associated spectral estimator achieves weak recovery of the latent subspace. Finally, we prove that the AMP-derived preprocessing of [Defilippis et al., 2025] is optimal among all bounded matrix-valued preprocessing maps of any fixed dimension. Its transition coincides with the AMP weak-recovery threshold, proving the general spectral conjecture of [Defilippis et al., 2025].
Florent Krzakala, Pierre Mergny, Vanessa Piccolo· 0 citations
A mathematical theory of superposition in neural networks using tools from frame theory and compressed sensing and a novel characterization of the distribution of signs in the Gram matrix is developed.
Michael I. Ivanitskiy, J. Jasper, Emily J. King et al.· 0 citations
We study online discrepancy minimization: vectors $v_1,\ldots,v_T\in\mathbb{R}^n$ arrive sequentially, and each must immediately be assigned a sign $x_t\in\{\pm1\}$, with the aim of minimizing $\|\sum_{t=1}^T x_t v_t\|_\infty$. We give a polynomial-time potential-based algorithm combining a regularization of the $\ell_\infty$-norm with restriction to an adaptively chosen coordinate set. For i.i.d. inputs with independent, symmetric, centered, unit-variance sub-Gaussian coordinates of sub-Gaussian norm at most $\sigma$, the algorithm achieves terminal discrepancy $O(\sigma^8\sqrt{n})$ with probability at least $1-\exp(-\Omega(\sigma^3\sqrt{n}))$. If the coordinates are independently masked by Bernoulli variables with mean $k/n$, where $k\gtrsim(\log n)^2$, the bound improves to $O(\sigma^8\sqrt{k})$, with failure probability $\exp(-\Omega(\sigma^3\sqrt{k}))$. Both guarantees hold for every prescribed finite horizon $T$, with no dependence on $T$. The dense result substantially generalizes a theorem of Bansal and Spencer (2020) for Rademacher inputs and gives an efficient $O(\sqrt{n})$ bound for Gaussian inputs, as conjectured by Gamarnik et al. (2022). When $T$ is polynomially larger than $n$, this is conditionally close to optimal: under worst-case hardness assumptions for standard approximate lattice problems, Vafa and Vaikuntanathan (2025) showed that no polynomial-time algorithm, even offline, can improve the $\sqrt{n}$ scale by a fixed polynomial factor in $T/n$.
