A Recovery Guarantee for Sparse Neural Networks
Authors: Sara Fridovich-Keil
Mert Pilanci
Venue: ICLR 2026 Poster (2026)
Keywords: compressed sensing
neural networks
model pruning
sparse weight recovery
Primary area: learning theory
TL;DR: We prove the first identifiability and recovery results for sparse ReLU neural networks, and show from-scratch recovery of sparse networks without first training a dense network.

Abstract: We prove the first guarantees of sparse recovery for ReLU neural networks, where the sparse network weights constitute the signal to be recovered. Specifically, we study structural properties of the sparse network weights for two-layer, scalar-output networks under which a simple iterative hard thresholding algorithm recovers these weights exactly, using memory that grows linearly in the number of nonzero weights. We validate this theoretical result with simple experiments on recovery of sparse planted MLPs, MNIST classification, and implicit neural representations. Experimentally, we find performance that is competitive with, and often exceeds, a high-performing but memory-inefficient baseline based on iterative magnitude pruning. Code is available at https://github.com/voilalab/MLP-IHT.

Published as a conference paper at ICLR 2026
A RECOVERY GUARANTEE
FOR SPARSE NEURAL NETWORKS
Sara Fridovich-Keil
Mert Pilanci
School of Electrical and Computer Engineering Georgia Institute of Technology s fkk afkkgatech.ede
Department of Electrical Engineering Stanford University pi pilanci dllanci@stanford.com @stanford.edu
ABSTRACT
We prove the first guarantees of sparse recovery for ReLU neural networks, where the sparse network weights constitute the signal to be recovered. Specifically, we study structural properties of the sparse network weights for two-layer, scalar- output networks under which simple iterative hard thresholding algorithm re- covers these weights exactly, using memory that grows linearly in the number nonzero weights. We validate this theoretical result with simple experiments recovery of sparse planted MLPs, MNIST classification, and implicit neural representations. Experimentally, we find performance that is competitive with, and often exceeds, high terfhrerfoniing but meemory-inefficie baseline based on iterative magnitude pruning. Code is available at https: / /gi tuub com/ voi 1 ab/MLPITHT
INTRODUCTION
Consider the task of training a sparse multilayer perceptron (MLP). We view this task through the lens of sparse signal recovery, in which the signal to be recovered is the vectorized MLP weights, most of which are zero so exact recovery requires finding the indices and values of the few nonzero MLP weights. Are these weights uniquely identifiable from training data? Can they be recovered efficiently in both memory and iteration complexity? For scalar-output, two-layer MLPs we answer both questions in the affirmatvve, proving what is to our knowledge the first recovery guarantee for sparse MLP weights.
Large neural networks are widely used as universal function approximators (Hornik et al., 1989), but as model size grows networks require ever larger memory and compute time to train (Kaplan al., 2020). Although large networks tend to be trainable to the highest quality, trained network weights are often highly compressible, e.g. by pruning, allowing for dramatic savings in memory and computation at inference time (Cheng et al., 2024). While sparse and high high-performing networks known to exist, efficiently optimizing them is open challenge. Existing approaches often compromise either memory efficiency reffiinncyrreuurringg memory to first train a dense network (Frankle & Carbin, 2019; Saikumar & Varghese, 2024; Gharatappeh & Sekeh, 2025)-or quality, failing to match the performance of dense counterparts (Frankle et al., 2021; Saikumar & Varghese, 2024). While some strategies can empirically balance efficiency and quality (Parger al., 2022; Jin al., 2016; Damadi et al., 2024), all existing approaches to sparse network training are heuristic in nature and lack formal guarantees of weight recovery.
At the same time, the compressed sensing literature is rich with theoretically-jusiiee algorithms to leverage sparsity in large-scal scale optimization tasks (see e.g., Wright & Ma (2022) for an accessible overview). However, these results are typically designed for linear models and convex optimization, and do not directly apply to recovery of sparse MI.P weights (Tropp, 2004; Khanna & Kyrillidis, 2018 8; Aghazadeh et al., 2018).
Our work bridges this gap by leveraging the recent development of convex reformulation of MLPs (Pilanci & Ergen, 2020a; Ergen & Pilanci, 2024), which allows to apply strong results from sparse signal estimation (Jain et al., 2014) to the task of training a sparse MLP. In its convex reformulation,
"This work initiated postdoc Stanford, and completed an assistant professor Georgia Tech.

---

Published as a conference paper at ICLR 2026
sparse MLP optimization can be viewed highly structured linear sensing problem in which the network weights are the signal to be recovered. We show that, when the training data consists of net- work evaluations random Gaussian sample points, this highly structured sensing matrix satisfies (with high probability) the classic restricted strong convexity and restricted smoothness conditions that suffice enable efficient sparse recovery via simple projected gradient descent method known as Iterative Hard Thresholding (IHT). Concretely, we make the following contributions:
We prove the first sparse recovery result applicable to ReLU MLPs, focusing on the case of shallow scalar-output network and random Gaussian data. Our result includes both unique identifiability of sparse network weights as well as a high-probability guarantee of efficient recovery of these weights via IHT, building result from Jain et al. (2014).
We demonstrate in a suite of illustrative small-scale experiments that IHT indeed tends to outperform a strong but
baseline of iterative magnitude pruning (IMP) (Frankle & Carbin, 2019), recovering gegher-performing sparse networks while using less memory during optimization. Our experiments include both 2-layer and 3-layer MLPs with both scalar and vector valued outputs, extending beyond the regime of our theoretical results.
2 RELATED WORK
2.1 SPARSE NEURAL NETWORKS
Prior work has shown that, in diverse contexts, large neural network may be well approximated by a sparse subnetwork, for example with only 10% of the original parameters left nonzero (Frankle & Carbin, 2019; Nowak et al., 2023). Sparse networks are far cheaper and faster to evaluate and store, making them attractive for applications on edge and resource-constrainet platforms as well as for democratizing access to large foundation models. Moreover, in many cases a sparse subnetwork can even outperform the prediction accuracy (Frankle & Carbin, 2019) and nut-of-diitrrbution robustness (Diffenderfer et al., 2021; Wu et al., 2024) of the original dense network.
However, sparse networks notoriously difficult to optimize. Existing approaches to finding sparse networks fall into three categories: iterative pruning (Frankle & Carbin, 2019; Liu et al., 2024), prun- ing initialization (Wang et al., 2022; Frankle et al., 2021), and dynamic sparse training (Jin al., 2016; Ji et al., 2024; Nowak et al., 2023; Damadi et al., 2024; Kusupati et al., 2020). Respectively, these approaches tend be high-performing but require high memory during optimization, memory and computation efficient to optimize but with reduced final model performance, and efficient but heuristic optimize reasonable final performance. None of the existing sparse network optimiza- tion paradigms come with theoretical understanding or recovery guarantees.
Prior theoretical results for sparse neural networks are present in Boursier & Flammarion (2023) and Ergen & Pilanci (2021) (see Lemma 10 therein), which derive conditions under which the sparsest two-layer MLP may be recovered by minimizing the Euclidean norm of the weights (i.e., applying weight decay). However, Boursier & Flammarion (2023) focuses on univariate data and Ergen & Pi- lanci (2021) considers sparsity of the second (output) layer weights, whereas our analysis considers arbitrary data dimension with focus sparsity of the first (hidden) layer weights. The recovery result in Ergen & Pilanci (2021) also requires fewer data points than dimensions, while our result does not. Further, the conditions in Ergen & Pilanci (2021) are based on the KKT optimality condi- tions of a semi-infinite convex formulation (Hettich & Kortanek, 1993) and are not straightforward to verify, nor is tractable recovery algorithm presented in Ergen & Pilanci (2021). In contrast, guarantee of sparse weight recovery relies on verifiable and satisfiable conditions that we show hold with high probability under random training data, and we prove that an iterative algorithm (iterative hard thresholding) achieves successful recovery of sparse neuron weights.
2.2 CONVEX NEURAL NETWORKS
Recent work has revealed an equivalence between training shallow (Pilanci & Ergen, 2020a) or deep (Ergen& & Pilanci, 2024) neural networks and solving convex optimization problems defined by network architectures. The core idea involves enumerating or sampling neuron activation paths to form fixed dictionary, whose coefficients are optimized via convex programming.
2

---

Published as a conference paper at ICLR 2026
Specifically, two-layer ReLU network approximates labels y using the nonconvex form y x where U = [u1..... Up] and v are the network weights and X the data ma- trix. Instead, the convex formulation uses activation patterns D; = Diag(l[Xx 2 0)) enumerated over all u to express the same network as
P Σ DiX(W; - -w,, 1=1
y R
(I)
2 0 and 1w) 2 0 for all i. Optimal values of the nonconvex weights U and can be recovered from optimal values of the convex optimization parameters t and w. Note that we use the term activation pattern to refer to a binary pattern whose length matches the number of training examples, and whose values denote which training examples are attended to by particular neuron (each neuron has its own activation pattern). The total number of activation patterns P derived from all possible is bounded exponentially in the data rank r, typically requiring subsampling for computational tractability. However, assuming sparsity in weights dramatically reduces the number possible patterns, cnabling exact convex optimization for large-scale datasets. Section 3 describes how we adapt and specialize this convex MLP reformulation for sparse networks in our theory and experiments.
subject to (2D,-1,)
2.3 ITERATIVE HARD THRESHOLDING (IHT)
Iterative Hard Thresholding (IHT) is special case of projected gradient descent, in which the pro- jection is onto the nonconvex set of sparse vectors. For large-scale sparse recovery problems, IHT and additive algorithms such basis pursuit and matching pursuit (Tropp, 2004) are often the only feasible algorithms, due to their memory efficiency compared to convex relaxations such as LASSO. IHT is also well-studied theoretically and comes with convergence guarantees both in its classic implementation (Blumensath & Davies, 2009; 2010; Jain et al., 2014) and accelerated variants (Blu- mensath, 2012; Khanna& & Kyrillidis, 2018). Some results also exist for variant of IHT augmented with a count sketch data structure (Aghazadeh et al., 2018), which can expand the regimes of spar- sity under which IHT enjoys successful recovery. Of these theoretical results for sparse recovery by IHT, most require the measurement matrix to satisfy either the restricted isometry property (RIP) with a small enough RIP constant, or restricted strong convexity and restricted smoothness proper- ties with small enough condition number; these conditions are too strict for the sparse MLP weight recovery task we consider.
However, Jain et al. (2014) proved more general sparse recovery result for IHT, showing recovery under restricted strong convexity and restricted smoothness with an arbitrary finite condition number. Jain al. (2014)'s result holds for classic IHT with the relaxation that the hard thresholding step of IHT must project onto a larger sparsity level than that of the true signal, where the inflation factor grows with condition number. Our theoretical results build on this result show that the task of recovering sparse MLP weights can be reformulated so as to satisfy the restricted strong convexity and restricted smoothness properties in expectation over Gaussian data, allowing us to show that IHT is guaranteed to recover the weights of a planted sparse MLP.
3 PRELIMINARIES
Consider a ReLU neural network with vector-valued input, scalar output, and a single hidden layer. We use X Raxd to denote the (Gaussian) data matrix with data points and data (input) dimen- sion d. We denote the ground truth labels or values as y E R", and the neural network output as ÿ E R2, The hidden weights of the 1-hidden-layer MLP are denoted U E Rdd where p is the width of the hidden layer. The columns of this weight matrix are u, E RRd = 1,....p, P, and the second layer weights U11
We now describe convexifying the model by fusing the first and second layer weights of the non- convex ReLU model D/ vz. We can express this model follows:
UjU1
أ = (diag (I(Xu1 / 2 0})X diag (I[Xup 2 0})x
(2)
E R
Up
3

---

Published as a conference paper at ICLR 2026
where use I(x 2 0} to denote the elementwise indicator function, taking value 1 at indices where 24 2 0 and value otherwise.
Training (e.g. with MSE loss) the 2-layer MLP in Equation (2) presents nonconvex optimization problem, because the parameters uj appear in both the weight vector and the A matrix. We convexify by simply replacing the p weight vectors u, in the A matrix with p separate, fixed generator vectors h; E RR and fusing the weights via w; = u Vi. This parameterization was previously studied in Mishkin al. (2022), where it shown to yield the gated ReLU (GReLU) network class, which is equivalent in expressivity to standard ReLU networks. Here, we extend this approach and show that sparse ReLU networks can also be recovered using similar strategy. We obtain
W
-= diag (I(Xh) 2))X diag (I(Xh₂ 2 0))X]
Wp
in this formulation exact recovery amounts to finding the sparse vector w* RP whose values the weights of a ground truth, planted MLP. If we allow the effective hidden dimension to be very large (up to 2d( ) (Pilanci & Ergen, 2020b)), we can choose a set of vectors h2 such exactly the set of all possible distinct activation patterns achievable for dataset X. Recall that in our notation, the term activation paltern refers to a binary pattern whose length matches the number of training examples n. and whose values denote which train- ing examples are attended to by particular neuron (each neuron has its own activation pattern). Moreover, for sparse neural networks with most s' nonzero weights per hidden neuron, we have p I 2s' (m)* by a counting argument; this may far fewer total activation patterns than needed to model dense weights. Consider neuron whose weight vector has at most s' nonzero entries. First, the support of this weight vector must be selected, which corresponds to choosing s' input dimensions out of d, resulting in (d) possible choices. For each choice of these s' dimensions, the neuron computes a linear threshold function in an s' dimensional subspace of IR.. A classi- cal result in the theory of hyperplane arrangements (Stanley al., 2007) shows that such linear
that ((T(X); > 0})}
threshold function can generate most 2 "''0" ("T)) < ((x) distinct activation patterns n data points. Multiplying the number ways to select the support and the number of patterns per support, and incorporating factor of s' for indexing neurons, we arrive the stated bound:
p S 2s' (:) (#)". With this large but fixed set of generator vectors h; we can solve similarly large but convex program to recover hidden weights w, corresponding to the globally optimal 2-layer nonconvex MLP.
Alternatively, we can operate with an arbitrary hidden dimension m and select the generator vectors h; at random such that the activation patterns {(Il > 0}) }m) are a random subset (drawn without replacement) of all possible activation patterns. In our theoretical results (Section 4) we assume patterss are enumerated; in experiments (Section 5) we sample S p patterns using random generator vectors.
THEORETICAL RESULTS
Consider the sparse recovery problem defined by Equation (3) of the form y Aw" for some unknown vector w", with sensing matrix
A == [diag (TT(M), 0))
diag (T(Xh, 2 0)) gaxdy,
Our main result leverages connections between sparse recovery methods and convex formulations of ReL.U networks. For simplicity, we will assume that the data matrix E Raad has entries drawn i..d. (0, 1); similar effect may be achieved in practice by data whitening. We also assume that the columns of A are unit-normalized before optimization.
To recover the planted weights w", we consider the following simple variant of the classic Iterative Hard Thresholding (IHT) algorithm.
w k+ = Hy aaA'"A"= - v))-
(4)
Here 17 > is a step size parameter and the hard thresholding operation Hş is a projection onto the set of 8-sparse vectors, where following Jain et al. (2014). In Lemma I we show that A
4

---

Published as a conference paper at ICLR 2026
satisfies restricted strong convexity and restricted smoothness with high probability over the random data X, making the sparse MLP weights uniquely identifiable. In Theorem we show that IHT efficiently recovers these sparse MLP weights.
Au is the planted neural network model. Recall that the relation between the standard and fused form of the weights is w* [uivt, u525 where sign (X,) sign(* X ut) Vi as defined in (2) and (3), Assumption 1 gives conditions on planted network under which we can ensure exact recovery of its weights.
Suppose that = (x,((x
Assumption (Properties of the planted sparse network). Assume that either
(a) u (11100,13 llu - ki of E RVvi E
[p] and kp S.a
(b) u TR. ||u Ilu = S; E [Smin k), u E (- (-1, [p] and Ei= 8; S holds.
Both parts of Assumption have to do with what values the planted MLP weights can take, and both parts restrict the number of nonzero hidden weights. Assumption (() requires that the nonzero hidden weights take binary values, but allows the output layer weights to take any real values. As- sumption ((b) captures the more relaxed and common scenario in which the nonzero hidden weights take any real values, but the output layer weights are restricted 11, since the flexibility to model any real value is already captured by the hidden layer weights.
In Appendix E show that Assumption 2 follows from either of Assumption ((a) and 1(b) with high probability, and we give weight constructions that satisfy each option in Assumption 1. We note that only Assumption 2 used in our proof of convergence and sparse recovery, Assumption is sufficient for Assumption 2 but may not be necessary. Likewise, we show that Assumption 2 is sufficient for sparse recovery but we do not prove that is necessary.
Assumption 2 (Properties of activation patterns). Let D; diag (1, > 0}) E RTX" with {D.Yi as the set of all such distinct activation patterns possible with data Raxd, whose entries are drawn i.i.d. ~ (O, 1). We assume the following properties about this set of enumerated activation patterns:
1. D; 2 En for all E [p], some E (0,1).
2. For all + i, the diagonals of D; and Dyr differ at least "yn positions, for some yf E
(0. 1).
In the appendix we prove that both of these hold with high probability under Assumption 1. Specif- ically, Assumption 2.1 holds with probability at least - pe pe- - -H(E)))(()) as long as 2 4k. Here H denotes binary entropy. Assumption 2.2 follows from Assumption I(a) with probability at least 2e-ee =05nn as long as n 2 C6) and S .(7+5)- Here e and C are positive ab- solute constants, 5 > O, and w () is the normalized Gaussian mean width of a subset K C TRd,, where K represents the set of (normalized) newron weights that satisfy Assumption (a)) Assump- tion 2.2 likewise follows from Assumption 1(b) with probability at least 1 - 2e 22--2-ff -o----f-- - ē as long as n 2 C8 15-0w(KK ww)) where now K represents the set of (normalized) neuron weights that satisfy A Assumption 1(b). Note that ē and some odditional restrictions on 8min and k are described in the appendix proof.
0.69
(Sample complexity). Note that Assumption 2 requires the number of training examples n 2 max(4k; C8-6wc K)), where ls is the sparsity level of each neuron, is the set of (normal- ized) neuron weights that satisfy Assumption 1(a) or I(b), w(K K) is its normalized Gaussian mean width, and C is a positive absolute constant. This is modest requirement that grows with the num- ber of active (nonzero) neuron weights rather than the total number of neuron weights, enabling compressive sensing of sparse newron weights.
Remark
Below we show that Assumption 2 is sufficient to ensure recovery of sparse MLP weights, for a 2-layer scalar-output ReLU MLP. Intuitively, the first part of Assumption 2 requires that every neuron attends to at least an e fraction of the training data, rather than fitting or overfitting to a tiny number of examples. Since our data covariates are assumed Gaussian, this first part of Assumption enables a concentration argument. The second part of Assumption 2 requires that any two different neurons must attend to subsets of the training dataset that differ by at least y fraction. Without
5

---

Published as a conference paper at ICLR 2026
this requirement, neurons might be very similar to each other and thus more difficult to distinguish and recover correctly during optimization. This second portion of Assumption bears similarity in spirit with the incoherence property common in compressive sensing.
Lemma (Restricted strong convexity and restricted smoothness). Let A E Ruxdp be as defined in Equation (3), with the modification that all columns are normalized to have unit lz norm. Assume that entries of the data malrix X E "kx are drawn i..d. N(0,) 1) and Assumption holds. Consider n,, and the induced - - 2s(s - 1)exp
an index set S S [dp] with |S|
matrix A As- For any ô
(=c5n>/'))) n2/440
(0, I-7), with probability at least 1 - 20(s -1) 1)exp
(s(s exp
(ee5ee
146
al, 1 A5 As 3 BI,
with 2 1
1 5 5
n1/4(1 8) 35 1-
8
1 5
11
1 5
n n1/4(1 5)
Here and y are the same as in Assumption 2, and is a positive universal constant.
Lemma ensures that the condition number of A, restricted to any of S n columns, is finite and bounded above by vB/a. The condition number shrinks as Y grows, because this enforces greater separation (incoherence) between columns of A. The conditioning worsens with increasing s, as this increases the number of columns in As and thus the potential for a coherent pair of columns. Theorem ensures that IHT recovers planted sparse weights regardless of this condition number (as long as it is finite), though the rate of convergence slows with increasing condition number.
Theorem (IHT recovers sparse MLP weights). Suppose that Assumption holds, the data matrix X E Rnd has entries drawn i.idd NO0,, 1), the activation patterns D; = diag (I{Xh; 2 03) in the sensing matrix are enumerated to include all unique patterns that can result from Ihllo S s, the columns of A are pre-normalized in l norm, and the planted neural network weights satisfy ||ww S 8. Consider the following variant of Iterative Hard Thresholding (IHT) to minimize the MSE objective f(w*) = Aw y y122.
w k+1
H3& HH(w"
"AA - -y)),
(5)
where 3 22) 17 33. and a, 3 are the restricted strong convexity and restricted smoothness constants from Lemma corresponding to sparsity level 23 + With the same high probability as
2
in Lemma 1, after K = ( _ log ([("))) steps, IHT finds sparse weights w such that
f(w")- SE and - w0 S -
Remark 2. Theorem I is, to knowledge, the first sparse recovery result that applies to sparse neural network weights. It extends Lemma to show that sparse MLP weights are not only uniquely identifiable with high probability from network's 's behavior on random data, but that these sparse weights may be recovered efficiently by IHT with high probability. If the underlying function map- ping data points to values is indeed a planted sparse MLP, recovery of the weights of this sparse MLP also guarantees generalization, in the sense that the labels of fresh data following the same function will be perfectly predicted by the sparse MLP recovered by IHT.
Proofs of all theoretical results may be found in the appendix.
5 EXPERIMENTAL RESULTS
Our experiments compare the performance of IHT and a strong MLP-pruning baseline method, iterative magnitude pruning (IMP), the algorithm from the Lottery Ticket Hypothesis (Frankle & Carbin, 2019), at training sparse MIPs. While these experiments are intended to complement and validate sparse recovery theoretical results, they also extend beyond the setting of Theorem I in several respects, to demonstrate that IHT empirically recovers high-performing sparse MLPs even under more flexible settings than those for which we can prove sparse recovery succeeds.
Specifically, our range of experiments for IHT include (i) both full-batch (deterministic) and mini- batch (stochastic) gradients, (ii) both scalar and vector-valued MLP outputs, (iii) both single-hidden- layer and deeper MLPs, (iv) both vanilla and accelerated IHT, (v) randomized (rather than enumer- ated) initialization for the sensing matrix A, for computational efficiency, and (vi) sequential convex

---

Published as a conference paper at ICLR 2026
updates A during IHT, rather than keeping A fixed as do in theoretical analysis. These sequential convex updates interpolate between the fully convex formulation in our theory and the nonconvex training that is standard practice for MLPs, enabling empirically strong performance for compared to what is required in Theorem 1.
IHT even with a much smaller, randomly-initialion The extension of IHT to vector-oupput ML.s and deeper MIPP is enabled by employing count- sketch datastructure (following Aghazadeh et al. (2018)) for noisy but memory-efficient estimation of all weights, slightly relax the hard thresholding in vanilla IHT (which we do use for shallow, scalar-output MLPs closer to the setting of our theoretical guarantees). We also strengthen the IMP baseline by pruning only 10% of the weights in each iteration (rather than the default 20%), which allows IMP to spend extra time finding a higher-performin experimental settings for both IHT and IMP are provided Appendix and our open-sourced code. We also provide some ablations and variability experiments in Appendix B.
sparse network. The details of our
We present experimental results on three illustrative tasks: fitting a planted sparse MLP, classifying handwritten MNIST digits (Deng, 2012), and fitting implicit neural representation to MNIST and CIFAR-IO 10 images (Krizhevsky et al.. 2009). In each task, we compare the performance of IHT (ours) and IMP (Frankle & Carbin, 2019), implemented described in Appendix A. For all figures, we show heatmaps comparing model performance as function of the hidden dimension m (vertical axis) and sparsity level (horizontal axis). We emphasize that IMP requires first training dense MIP and then iteratively pruning it to achieve sparse weights, whereas IHT optimizes sparse weights directly and thus has far smaller memory requirements during training.
Results on scalar-output and vector-output planted sparse MLPs are presented in Figure 1 and Fig- 2, respectively. Within Figure and Figure the left subfigures compare IHT and IMP on 2-layer (1-hidden-layer) sparse MLPs while the right two subfigures compare IHT and IMP on 3-layer (2-hidden-layer) sparse MLPs. For the planted MLP fitting tasks, we optimize sparse MLP with hidden dimension (heatmap vertical axis) and a budget of s nonzero weights (heatmap horizontal axis) to match the input-output behavior of unknown planted model of the same ar- chitecture and sparsity. Specifically, we draw random sparse weights and use these to generate dataset (X R g500000 1/ E IR (2,0000.c), for c 10). For planted MLP fitting tasks, we report peak signal to noise ratio (PSNR) at fitting this input-output behavior of the planted model. PSNR defined as PSNR = 10 log ( /MSE), where mms the largest magnitude value in the ground truth signal and MSE is the mean squared error; higher PSNR is better. If certain setting of and s yields a planted model whose outputs y are all zero, we skip evaluation and report PSNR of zero.
30
43
20
20
10
100
10
10
200
200
9
-38
-10
-10
100
20 1000
-20
IHT
IMP
IHT
IMP
Figure 1. Average PSNR for fitting a planted one-hidden-layer (left) and two-hidden-laye (right) sparse scalar-output MLP of hidden dimension (vertical axis) and most nonzero parameters (horizontal axis). Colorbar shows average PSNR over random trials. exhibits more robust performance than a strong but memory-inee memory-inefficient iterative magnitude pruning (IMP) baseline (Frankle Carbin, 2019).
Figure 3 presents results on MNIST digit classification. In Figure 3, all MIPs have one hidden layer; the left two subfigures compare IHT and IMP on binary classification and the right two subfigures compare IHT and IMP 10-way classification. We consider both binary classification (digit 0 vs. 1) posed as a regression problem with MSE loss, as well as 10-way classification of all digits using cross-entropy loss and one-hot labels 10 For MNIST classification, report classification accuracy, the fraction of test digits correctly classified.
Although runtime is not key component of our analysis experiments, trend worth noting is that the runtime for IHT is increasing in s, whereas the opposite true for IMP. This is because IMP requires iterative retraining with gradual pruning, so more steps of retraining required to reach

---

Published as a conference paper at ICLR 2026
16
53
102
100
297
202
-0 122
-p
500
100
1391
3150
IHT
IMP
IHT
IMP
Figure Average PSNR for fitting a planted one- me-hidden-layee (left) and two-hidden-layy layer (right) sparse vector-output (10-dimensional output) MI.P of hidden dimension m (vertical axis) and at most & nonzero parameters (horizontal axis). Colorbar shows average PSNR over random trials. IHT is competitive with strong but meemory-inefficien iterative magnitude pruning (IMP) baseline (Frankle & Carbin, 2019).
30
20
18
20
38
100
50
26
100
300
200
500
300
36
1032
532
22
100
500
1A0
2000
A000
IHT
IMP
IHT
IMP
Figure Average binary (left) and 10-class (right) classification accuracy for handwritten MNIST digits with 2-layer (one-hidden-layer) MLP of hidden dimension m (vertical axis) and at most s nonzero parameters (horizontal axis). Colorbar shows average classification accuracy over 3 random trials. IHT exhibits more robust performance than strong but memory -inffri-nefficciet iterative magnitude pruning (IMP) baseline (Frankle & Carbin, 2019).
sparser network (with smaller s). In other words, IHT is fastest exactly where IMP is slowest. Runtime varies for both IHT and IMP as function problem parameters, so we provide a few illustrative examples, all evaluated on an NVIDIA A6000 GPU.
For binary MNIST classification using the smallest scalar-output model with m = (hidden layer has single neuron) and sparsity = (meaning that neuron can attend to single pixel only), and 15 full-batch gradient steps, IHT reaches 98.85% % test accuracy in 1.2 seconds, while IMP reaches 50% test accuracy (random chance) in 27.78 seconds. With = 10 and 100, IHT reaches 99.2% accuracy in 1.66 seconds; IMP reaches 50.15% in 20.56 seconds. With m = 100 and s = 1000, IHT retains 99.2% accuracy in 8.4 seconds, and IMP achieves 77.66% accuracy in 20.91 seconds. For small, scalar-output MLPs IHT is dominant in memory, runtime, and accuracy.
However, for vector-output MLPs, deeper MLPs, and settings with large s. IHT (in its current im- plementation) can run more slowly than IMP, especially when using minibatch gradient updates. For full (10-class) MNIST classification, with 50 epochs, batch size 5000, 88.73% 73% test accuracy in 219.2 seconds, whereas IMP gets 77.66% accuracy in 68.98 sec-
10, and 1000,
IHT gets 88. onds. For fitting planted MLP with 10-dimensional output. with 15 epochs, full-batch gradients, = 10, and = 10, IHT reaches 48.67dB PSNR in 3.02 seconds while IMP reaches 24.81 db PSNR 28.46 seconds.
Figure 4 and Figure present results on fitting MNIST and CIFAR-10 images, respectively, with an MLP-based implicit neural representation. Specifically, we use fixed Fourier feature embedding with Gaussian-distributed frequencies followed by ReLU MLP, following Tancik et al. (2020). We use the same Fourier features for embedding pixel coordinates for both IMP and IHT, and vary the optimization strategy for fitting the sparse MLP weights. Within Figure 4 and Figure the left two subfigures compare IHT and IMP on 2-layer (1-hidden-layer) sparse MLPs while the right two subfigures compare [HT and IMP on 3-layer (2-hidden-layer) sparse MI.s.
8

---

Published as a conference paper at ICLR 2026
120
102
IHT
IMP
IHT
IMP
Figure Average 1-hidden-layer (left) and 2-hidden-layer (right) PSNR for overfitting MNIST digit image with an MLP-based implicit neural representation (Tancik et al., 2020) hidden dimension (vertical axis) at most s nonzero parameters (horizontal axis). Colorbar shows average PSNR over 3 random trials. IHT exhibits more robust performance than strong hut memoretieefficient itcrative magnitude pruning (IMP) baseline (Frankle & Carbin, 2019). We highlight that IHT exhibits stable recovery independent of m. in line with our theoretical results (see Remark 1). In contrast, IMP shows improved recovery with increasing m, likely because IMP here is solving nonconvex optimization problem whose landscape is made more benign by increasing m.
190
IHT
IMP
IHT
IMP
Figure Average 1-hidden-layer (left) and 2-hidden-layer (right) PSNR for overfitting CIFAR-10 digit image with MLP-based implicit neural representation (Tancik et al., 2020) of hidden dimension m (vertical axis) and at most s nonzero parameters (horizontal axis). Colorbar shows average PSNR over random trials. IIT exhibits more robust performance than strong memory-inefficien itcrative magnitude pruning (IMP) baseline (Frankle & Carbin, 2019).
Across these experimental settings, we find that IHT almost always finds gigher-performing sparse network weights compared to IMP, a strong baseline for pruning nonconvex MLPs (Frankle et al., 2021). Moreover, IHT uses a fixed parameter budget that scales with the sparsity level throughout optimization, whereas IMP requires initial training of a dense network whose parameter count grows with data dimension d and hidden dimension m.
DISCUSSION
This work presents, to our knowledge, the first sparse recovery result applicable the weights of ReLU MLP. For nonnegative scalar-output MLPs, we show that sparse weights are uniquely iden- tifiable and efficiently recoverable from measurements on random Gaussian data, with high proba- bility. We complement this theoretical result with an empirical demonstration that a simple iterative hard thresholding algorithm can find sparser and higher-performing network weights compared to strong network pruning baseline, while using far less memory during training.
Limitations and future work. Our results subject several limitations that expect future work may address. Our theoretical results are restricted to shallow, scalar-output MLPs, and are shown to hold with high probability over Gaussian data rather than more general data distributions. As these are the first recovery results for sparse MLPs, we are optimistic that future work may extend our results to deeper, vector-output networks with more diverse architectures and data distributions. Our experiments also suggest that sequential convex optimization from random initialization can find high-performing sparse MLPs; extending our convex- monvex-formulation recovery result to sequential convex programming is of interest. Our IHT recovery result also inherits an inflated sparsity level from Jain et al. (2014); tightening this result a compelling direction for further study. Finally, encourage future work to refine and scale up implementation of memory- IHT training of sparse MLPs to enable both memory-efficien and fast training of high performenn sparse MI.Ps.

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
This work supported in part by the NSF Mathematical Sciences Postdoctoral Research Fel- lowship under award number 2303178, in part by the National Science Foundation (NSF) CA- REER Award under Grant CCF-2236829, in part by the National Institutes of Health under Grant 1R01 AG08950901AL, in part by the Office of Naval Research under Grant N00014-24-15 1-14--11144,, and in part by the Defense Advanced Research Projects Agency under Grant HR001 12490441. We are grateful to Gordon Wetzstein for helpful discussions on applications of sparse optimization.
10

---

Published as a conference paper at ICLR 2026
REFERENCES
Amirali Aghazadeh, Ryan Spring, Daniel LeJeune, Gautam Dasarathy, Anshumali Shrivastava, et al. Mission: Ultra large-scale feature selection using count-sketches. In International conference on machine learning, pp. 80-88. PMLR, 2018.
Dennis Amelunxen, Martin Lotz, Michael B McCoy, and Joel Tropp. Living on the edge: Phase transitions in convex programs with random data. Information and Inference: A Journal of the IMA, 3(3):224-294, 2014.
Yoshua Bengio, Pascal Lamblin, Dan Popovici, and Hugo Larochelle. Greedy layer-wise training of deep networks. Advances in neural information processing systems, 19, 2006.
Thomas Blumensath. Accelerated iterative hard thresholding. Signal Processing, 92(3):752-756, 2012.
Thomas Blumensath and Mike Davies. Iterative hard thresholding for compressed sensing. Ap- plied and computational harmonic analysis, 27(3):265-274, 2009.
Thomas Blumensath and Mike E Davies. Normalized iterative hard thresholding: Guaranteed stabil-
ity and performance. IEEE Journal of selected topics in signal processing, 4(2):298 -009,,,,-,, 2010.
Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. Concentration Inqqualities: A Nonasymp- totic Theory of Independence. Oxford University Press, 02 2013. ISBN 9780199535255. .0001. URL https: /tt.///iiii///// http://iiorrgrr/////e/ee 1093/
doi: 10. 00.0999//ppmtta000/7771999135333
acprof osso 97801 99535255. 10001000011.
Etienne Boursier and Nicolas Flammarion. Penalising the biases in norm regularisation enforces sparsity. Advances in Neural Information Processing Systems, 36:57795-57824, 2023.
Hongrong Cheng, Miao Zhang, and Javen Qinfeng Shi. A survey on deep neural network pruning: Taxonomy, comparison, analysis, and recommendations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2024.
Saced Damadi, Soroush Zolfaghari, Mahdi Rezaie, and Jinglai Shen. Learning a sparse neural network using IHT. In 2024 International Joint Conference on Neural Networks (IJCNN), pp- 1-8. IEFE, 2024.
Li Deng. The mnist database of handwritten digit images for machine learning research [best of the web]. IEEE signal processing magazine, 29(6): 142-142, 2012.
James Diffenderferr Brian Bartoldson, Shreya Chaganti, Jize Zhang, and Bhavya Kailkhura. A winning hand: Compressing deep networks can improve duiti--distiibution robustness. Advances in neural information processing systems, 34:664-676 676, 2021.
Tolga Ergen and Mert Pilanci. Convex geometry and duality of over-parameterizee neural networks. Journal of machine learning research, 22(212):1 2021.
Tolga Ergen and Mert Pilanci. Path regularization: A convexity and sparsity inducing regularization parallel relu networks. Advances Neurul Information Processing Systems, 36, 2024.
Tolga Ergen, E Candès, and M Pilanci. Random projections for learning non-convex models. In 33rd Conference on Neural Information Processing Systems, Workshop on Beyond First Order Methods in Machine Learning, 2019.
Jonathan Frankle and Michael Carbin. The lottery ticket hypothesis: Finding sparse, trainable neu- ral networks. In 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019. OpenReview.net, 2019. URL httph:ttt///...... httpt:://ppenreeeeeeee /owp://veeerr...eeeemmmmm. net net/Eorum?? 20-oummii==00--33c---- D3RCF7.
Jonathan Frankle, Gintare Karolina Dziugaite, Daniel Roy, and Michael Carbin. Pruning neural net- works at initialization: Why are we missing the mark? In International Conference on Learning Representations, 2021. URL https pttttpss//opereeeeee /pwp//pppeeeee.eoooe/e.e.. net
wtw.eeiwwwwweeeeeue/e/e forum? -og.rruunnttttttt-----------

---

Published as a conference paper at ICLR 2026
Soheil Gharatappeh and Salimeh Yasaei Sekeh. Information consistent pruning: How to efficiently search for sparse networks? In International Workshop on Al for Transportation, pp. 270-284. Springer, 2025.
Rainer Hettich and Kenneth O Kortanek. Semi-infinite programming: theory, methods, and appli- cations. SIAM review, 35(3):380-422 429. 1993.
Kurt Hornik, Maxwell Stinchcombe, and Halbert White. Multilayer feedforward networks are uni- versal approximators. Neural networks, 2(5):359-366, 1989.
Prateek Jain, Ambuj Tewari, and Purushottam Kar. On iterative hard thresholding methods for high- dimensional m-estimation. Advances in neural information processing systems, 27, 2014.
Jie Ji, Gen Li, Lu Yin, Minghai Qin, Geng Yuan, Linke Guo, Shiwei Liu, and Xiaolong Ma. Ad- vancing dynamic sparse training by exploring optimization opportunities. In Ruslan Salakhut- dinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp (eds.), Proceedings of the 41st International Conference Machine Learning, vol- ume 235 of Proceedings of Machine Learning Research, pp. 21606-21613 PMLR, 21-27 Jul 2024. URL https::/ htttp::/prroceedngs .ml.. wwwppppsssssseeee j12/ html.
Xiaojie Jin, Xiaotong Yuan, Jiashi Feng, and Shuicheng Yan. Training skinny deep neural networks with iterative hard thresholding methods. arXiv preprint arXiv:1607.05423 2016.
Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei. Scaling laws for neural language models. arXiv preprint arXhv:2001.0836/ 2020.
Ali Karimi, Ahmad Kalhor, and Melika Sadeghi Tabrizi. Forward layer-wise learning of convo- lutional neural networks through separation index maximizing. Scientific Reports, 14(1):8576, 2024.
Rajiv Khanna and Anastasios Kyrillidis. Iht dies hard: Provable accelerated iterative hard thresh- olding. In International Conference on Artificial Intelligence and Statistics, pp. 188-198, PMLR, 2018.
Sungyoon Kim and Mert Pilanei. Convex relaxations of relu neural networks approximate global optima in polynomial time. In International Conference on Machine Learning, pp. 24458-24485. PMLR, 2024.
Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009.
Aditya Kusupati, Vivek Ramanujan, Raghav Somani, Mitchell Wortsman, Prateek Jain, Sham Kakade, and Ali Farhadi. Soft threshold weight reparameterization for learnable sparsity. In International Conference Machine Learning, pp. 5554-5445555 5555. PMLR, 2020.
Bohan Liu, Zijie Zhang, Peixiong He, Zhensen Wang, Yang Xiao, Ruimeng Ye, Yang Zhou, Wei- Shinn Ku, and Bo Hui. A survey of lottery ticket hypothesis. arXiv preprint vrXiv:2403.0486/ 2024.
Aaron Mishkin, Arda Sahiner, and Mert Pilanci. Fast convex optimization for two-layer relu net- works: Equivalent model classes and cone decompositions. In International Conference Ma- chine Learning, pp. 15770-15816. PMLR, 2022.
Aleksandra Nowak, Bram Grooten, Decebal Constantin Mocanu, and Jacek Tabor. Fantastic weights and how to find them: Where to prune in dynamic sparse training. Advances in Neural Information Processing Systems, 36:55 160-5599, 2023.
Mathias Parger, Alexander Ertl, Paul Eibensteiner, Joerg H Mueller, Martin Winter, and Markus Steinberger. Gradient- based weight density balancing for robust dynamic sparse training. arXiv preprint arXiv:2210.1400 2022.
12

---

Published as a conference paper at ICLR 2026
Mert Pilanci and Tolga Ergen. Neural networks convex regularizers: Exact polynomial-time convex optimization formulations for two-layer networks. In Hal Daumé III and Aarti Singh (eds.), Proceedings of the 37th International Conference
Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 7695-7705. PMLR, 13-18 Jul 2020a. URL. htttps://prosee
ngs ml press / v1 9/ p http/pplaacciii.eee.ee.eee..eaae lanci20a html..
Mert Pilanci and Tolga Ergen. Neural networks are convex regularizers: Exact polynomial-time con- vex optimization formulations for two-layer networks. In International Conference on Machine Learning, pp. 7695-7705. PMIR, 2020b.
Yaniv Plan and Roman Vershynin. Dimension reduction by random hyperplane tessellations. Dis- crete & Computational Geometry, 51(2):438-46 4618, 2014.
Dhananjay Saikumar and Blesson Varghese. Drive: Dual gradient-basee based rapid iterative pruning.
arXiv preprint arXiv:2404.03687, 2024.
Richard P Stanley et al. An introduction to hyperplane arrangements. Geometric combinatorics, 13: 389-496 496, 2007.
Matthew Tancik, Pratul Srinivasan, Ben Mildenhall, Sara Fridovich-Keil. Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan T. Barron, and Ren Ng. Fourier features let net- works learn high frequency functions in low dimensional domains. NeurIPS, 2020.
Joel A Tropp. Greed is good: Algorithmic results for sparse approximation. IEEE Transactions on Information theory, 50(10):2231-2242, 2004.
Huan Wang, Can Qin, Yue Bai, Yulun Zhang, and Yun Fu. Recent advances on neural network pruning at initialization. In Luc De Raedt (ed.), Proceedings of the Thirty-First International Joint Conference on Artificial Intelligence, IJCAI 2022, Vienna, Austria, 23-29 July 2022, pp- 5638 5645. ijcai.org, 2022. doi: 10.2466//ICAA222227786 URL https:t/p///oomm /ooii org/ 110. 24963/ijcal 2022/786.
John Wright and Yi Ma. High-dimensional data analysis with low-dimensional models: Principles, computation, and applications. Cambridge University Press, 2022.
Boqian Wu, Qiao Xiao, Shunxin Wang, Nicola Strisciuglio, Mykola Pechenizkiy, Maurice van Keulen, Decebal Constantin Mocanu, and Elena Mocanu. Dynamic sparse training ver- sus dense training: The unexpected winner in image corruption robustness. arXiv preprint arXiv:2410.03030, 2024.
13

---

Published as a conference paper at ICLR 2026
APPENDICES
A EXPERIMENTAL METHODS
In this section we describe our experimental implementation of IHT as well as a network-pruning baseline algorithm, IMP (Frankle & Carbin, 2019). Our experiments are built mixture of Py- Torch and CuPy, and our code is available at https: / www.ttuubcoom//o/o
aww.aaaaaaaaaaaaaaaaaaa IHT.
IHT updates. Our experiments use the classic IHT update rule w*+1 = H, (w= m7 V f(=*)) where fw=") is the objective function to be minimized. We do not inflate the projection sparsity level to the š required in our theoretical analysis; doing so would likely further improve performance at the cost of inflating memory usage. For most of our experiments objective with gradient V VSMSE in our experiments we use hard thresholding to enforce s-sparsity but not r-structure in the neurons. However, for our experiments on multiclass classification, we use the cross-entropy objective whose gradient is V fce(w*) = A (softmax(Aur y))
the mean squared error (MSE) (AA - y)- This yields the update rule in Equation (4), where
Memory- IHT implementation. A key strength of IHT its memory efficiency, since only the nonzero weights and their indices need to be stored during optimization. However, achiev- ing this memory efficiency requires careful implementation because each gradient V f((=* is dense vector rather than a sparse one, and because the sensing matrix A is huge. Our implementa- tion is therefore blockwise. Instead of storing the entire matrix A dx we generate each n x block on the fly it is needed. Instead of computing the entire gradient (ww")) E 4R at once, we compute each d-dimennional block and apply it to the sparse iterate ++1 before computing the next block. Mathematically this is equivalent to computing the entire gradient and performing a single hard thresholding, but can be far more memory efficient. The choice block size is a design parameter that allows our IHT implementation strategy to trade off memory and computation time, allowing large sparse models to be trained under diverse hardware constraints.
Sequential convex IHT. In our theoretical analysis, we assume that the p activation patterns (I(X,, 2 0}) are enumerated to include all possible unique activation patterns based fixed, sparse generator vectors h;, This construction produces a large-scale but convex and very sparse optimization problem. In our experiments, for computational efficiency we instead solve sequen- tially convex optimization that switches between the convex formulation in Equation (3) and the nonconvex formulation in Equation (2). We choose fixed hidden dimension m (rather than larger enumerated dimension p) for the network weights, and frequently update the construction of the sensing matrix A E Raxdm maintain consistency with the weights w E Rd as they evolve during optimization, starting from a random initialization. In between these updates to A the for- mulation is fixed and convex, hence the terminology of sequential convex optimization. We can equivalently view optimization in this sequential convex formulation system in which the sensing matrix A is really Ag, as depends on the current weight estimate w We find the best performance arises from two-stage optimization procedure in which we hold the generator vectors inside A fixed at their random initialization until completion of the first epoch (pass through the training dataset), and then allow the generator vectors to update after each subsequent IHT iteration. Intuitively, this procedure stabilizes the first phase of optimization by maintaining convexity, and then allows for refinement of the sensing matrix once IHT has had the opportunity to enter a region of attraction around the global optimum.
time-varying dynamical
Vector-output MLPs. The formulation in Equation (2) and Equation (3) assumes a scalar-output MLP in which the output layer can be fused to the hidden layer weights. In MLP with vector- valued output, we instead have
(XW)+ W. where before X Red and W Rddm for hidden dimension m. but now ò R and W mm for output dimension c. We can no longer fuse the output layer weights, so we optimize the following formulation:
Wi W
U [diag (L(Xw) 2 0})X
diag (1{ Wm 2 0})x
(6)
Wm T
14

---

Published as a conference paper at ICLR 2026
where W, E Rd is a column of W (as in the scalar-output case) and i;; E R is row of W. We use the chain rule to compute separate gradients for W and W. computed blockwise and applied on-the-fly to update a global sparse weight representation for memory efficiency.
Layerwise optimization for deeper MLPs. Although it is possible refine the formulation Equa- tion (2) for deep MLPs (Ergen & Pilanci, 2024), for simplicity of implementation we optimize deep MLPs following layerwise approach (Bengio et al., 2006; Karimi et al., 2024). For example, to optimize 3-layer MLP (2 hidden layers plus an output layer), we proceed as follows. First, we optimize a 2-layer, vector-output MLP to find sparse weights W E R and W E R We then discard W and freeze W. treating it as an input embedding while training second 2-layer MIP, this time with input dimension d = m. We note that our results for IHT on deeper MI.Ps are slightly pessimistic, our optimization procedure allocates some nonzero parameters intermediate output layers W that are not used in the final model, meaning that the final model performance is attained with strictly fewer active weights than the budgeted s. Nonetheless, IHT remains competitive despite this restriction (see Section 5).
Count sketching. For shallow scalar-output MLPs, we use the standard hard thresholding rule based weight magnitudes, retaining the gighest-magnitude entries in W at each iteration. For deeper and vector-output MLPs, we follow Aghazadeh et al. (2018) and use an intermediate count sketch data structure perform hard thresholding. Intuitively, we view the count sketch approach as a noisy but more memory efficient alternative to the deterministic sparsity inflation in our theo- retical analysis. At every gradient step. we update the count sketch maintain a noisy estimate of the full iterate, a vectorized concatenation of W and W. We use a vector of dimension 4s log(n/ /s) to represent the count sketch, balancing memory efficiency with the level of approximation error in the sketch. At each iteration of IHT, we find the s entries in the count sketch with largest estimated magnitude, and store exact values for these entries. We observe an empirical tradeoff in the use of a count sketch: for shallow scalar-output networks where we have a single weight matrix W to optimize, the approximation error introduced by the count sketch outweighs any benefit it brings by "softening" the hard thresholding operation. However, for vector-output networks or deeper net- works where IHT must implicitly decide how to allocate a fixed parameter budget among W and W, the count sketch allows IHT to make less myopic thresholding decisions that aggregate information from multiple gradient steps, the benefits of which appear to outweigh the cost of approximation error in the count sketch.
Step size selection. For our experiments with IHT, we use two different step size selection meth- ods. For shallow scalar-output networks, we fuse the output layer weights following Equation (2). We can then compute an adaptive stepsize to minimize the mean squared error (MSE) objective function
AT Auppi w (v Aw*)
1i6
|| A Auppp(u AT (y Aw
following Blumensath & Davies (2010). However, as Equation (7) does not directly apply to vector- output networks, for these we use a fixed stepsize 77s for both W and W, and manually tune n7- This manual tuning surely leaves room for improvement with adaptive strategy, which we defer to future work.
Accelerated IHT. For shallow scalar-output networks, we use an accelerated IHT following Blu- mensath (2012), which defines an accelerated IHT as any algorithm that augments the classic IHT update with refinement step to produce an iterate that both sparse and has objective value no larger than that of the iterate produced by classic IHT. Specifically, after each IHT update we take a few gradient steps restricted to the current set of nonzero weights, to lower the objective value without changing the sparse support (Blumensath, 2012). We do not find an empirical benefit to this acceleration procedure for vector-output MLPs, so we perform accelcration only for IHT on scalar-output networks.
IMP baseline. We compare our IHT approach for optimizing sparse MLPs with iterative magni- tude pruning (IMP) (Frankle & Carbin, 2019), high-performing baseline method for pruning neural
15

---

Published as a conference paper at ICLR 2026
networks that has been shown to find sparse networks that often match or exceed the quality of their dense counterparts. IMP begins by training a dense network, and then iteratively prunes (sets to zero) constant fraction of the active (nonzero) weights based magnitude, rewinds the remaining active weights to their initialization values, and retrains. IMP thus allows pruning a dense network to any desired sparsity level, but requires sufficient memory train the dense network and sufficient computation time to iteratively retrain it during pruning. Although this IMP process is computation- ally costly, it provides strong baseline of performance that can be achieved with sparse network using existing methods. Frankle & Carbin (2019) suggest pruning 20% of the active weights at each iteration to balance computation and performance; we prune only 10% of the active weights at each iteration to maximize IMP performance and provide as strong baseline as possible.
Minibatches. For both IHT and IMP, each iteration operates on a minibatch of the full dataset. For IHT, we perform minibatch update by subsampling the rows of both A and y. For fitting planted sparse MLP as well as for 10-way MNIST classification we use a minibatch size of 5000, and for binary MNIST classification we use minibatch size of 1000. These settings all correspond to 10% of each training dataset per minibatch. For fitting an implicit neural representation to MNIST and CIFAR- 10 images, we use full-batch updates as each image small.
B METR-EXPERIMENTS
Figures 6 to 10 parallel Figures to but include standard deviation over the three random trials; high-variance experiments tend to align with empirical phase transitions.
As our experiments use sequential convex updates to the matrix A whereas our theoretical analysis assumes A is static, we empirically evaluate the convergence of A under our sequential convex updating strategy. Specifically, we consider the following setting: fitting a scalar-output planted sparse 2-layer MIP with input dimension 100, hidden dimension m = 10, and sparsity level 8 500, optimized over 100 steps of IHT with random initialization and sequential convex updates to A. We update the A matrix every l steps (in our main experiments & 1), and report how the PSNR (fitting quality metric, higher is better) changes as a function of k. Note that when = 100, A is never updated because this experiment only runs for 100 steps. Results are reported in Table 1. We see that when A is updated fairly frequently (k < 5) IHT converges and matches the planted model up to numerical precision. However, as A is updated less frequently, the model no longer converges within 100 steps, and PSNR degrades with increasing k.
We also evaluate the convergence of A more directly, by checking the first step which the support of A (which columns are active) stops changing; results are reported in Table 2. We find that as k increases, A takes longer to converge. As long as 5, A still converges within 100 steps. How- ever, when i- 5. A does not converge within 100 steps, likely explaining why there is dramatic drop in PSNR in Table when ie increases from 5. we instead continue optimization longer when k- 5, A does converge at step 135, after which the PSNR increases to 161 161.44 (correspond- ing to machine precision for this task). From this set of experiments, we see that the frequency of sequential convex updates to A controls the convergence rate of IHT. This experiment also suggests a natural stopping criterion for IHT: when the support of A stops changing (though in our main experiments we use fixed number of steps, for fair comparison to IMP).
C PROOF OF LEMMA
Proof. Our goal is upper and lower bound the eigenvalues of A7 E Rexe, where the full matrix A E dx is column column-normalized version of
Idiag (T{ (11, 2 0))x
diag (T(Xh; 2 0])X] Rxddp
X R has i.ddd (O, 1) entries, and the index set S C [dp] has |S| = S n.
As singular values of As are unaffected by column permutation, without loss of generality, we assume that columns of As are ordered so that columns involving each 2y are adjacent to each other. This ordering induces block structure in AS we refer to the i, block submatrix as As)))) Using this block structure, we bound each type of entry in AI As with high probability
16

---

Published as a conference paper at ICLR 2026
PSNR
161.44
161.44
161.44
661.44
36.63
6
31.85
31.84
28.49
29.88
10
28.38
15
27.23
20
25.85
25
24.75
30
24.74
35
23.32
40
23.35
50
22.02
60
21.95
70
21,96
80
21.97
90
21.89
98
21.57
100
19.88
Table 1. Fitting a scalar-output 2-layer MLP with input dimension 100, hidden dimension mi 10. and sparsity level s 500, optimized over 100 steps IHT with random initialization and sequential convex updates to A every k steps. In our main experiments k: = 1; here we observe that fitting quality degrades with increasing k:
convergence step
43
2
64
75
84
Table 2. Fitting a scalar-output 2-layer MLP with input dimension 100, hidden dimension m 10, and sparsity level s 500, optimized over 100 steps of IHT with random initialization and sequential convex updates to A every steps. In our main experiments = 1; here we observe that A converges more slowly with incrcasing k. We report "convergence step" as the first step of optimization at which the set of active columns of A does not change, signaling support recovery.
and then use these entry-wise bounds to upper and lower bound the eigenvalues of A with high probability.
First, consider block submatrix (A].As),, on the diagonal of ATAs. Since the columns of A are normalized, the diagonal entries of As), are deterministically I for all i. The off off-diagonal en- The numerator has Eis? Dy Tr Tr D,Dy, tries of (s) take the form [D P,
xx D,,,D,r,,
while the denominator terms have expectation Tr D, and Tr Dy, respectively. Applying -rinhn (Theorem 2) to each quadratic form in this expression, we have
x DDDDD
Tr DrD +1
< TT t2
(8)
[D 11guell
Ti Dy t3
-ot? (DD,D,,47))
d2 Ty TDD
ct2
with probability at least 2exp universal constant
- 2exp
2 exp
(TD,,+1g))
for
Consider two regimes based on whether Tr Dy is less than or greater than E eV In the first regime, Tr D, Dy S e/n. We choose EVn, la 5 Tr Dj, and la 5 Tr Dy for 8 (0, 1),
17

---

Published as a conference paper at ICLR 2026
100
200
2000
100
100
-0
1000
1099
IHT
IMP
IHT
IMP
Figure Standard deviation over the three random trials (top row) for each experiment reported in Figure (bottom row).
200
200
100
200
500
500
1000
2000
IHT
IMP
IHT
IMP
Figure 7. Standard deviation over the three random trials (top row) for cach experiment reported in Figure 2 (bottom row).
which yields
of DyDe
Tr D2 Dy. teVi
<
|D,0,lll2 [D,00.112
V(1 (1-68 8) D1 0r
(c0 Tr D2) -2exp (
-00² D.).
with probability at least
2 2xxp
2 2exp
1+8
140
Since we are in the regime where Tr DrDDD S evn and, by Assumption 2, D, and Dj each have trace least en, we have
of DDDyyi:
2
D 112 Dr
[1 - 8)
with probability at least eexp
(=rrva)
exp c
18

---

Published as a conference paper at ICLR 2026
0.000
50-
20
0.006
006
0.004
630
0004
100
004
200
200
800
2000
10
03
09
06
2.6
100
200
3000
500
06
3000
200
500
06
500
1000
IHT
IMP
IHT
IMP
Figure 8. Standard deviation over the three random trials (top row) for each experiment reported in Figure (bottom row).
28
IHT
IMP
IHT
IMP
Figure 9, Standard deviation over the three random trials (top row) for each experiment reported
Figure 4 (bottom row).
In the second regime, Tr D,Dy > evn. In Equation (8) we choose f = TTrr 8T Dj. and 13 Dy,, for (0, 1). to yield
x. DDD |D 11z 11
(1 (1+ 8) Dr / Dr
-08 Tr I+3 out loss of generality, assume that Tr D3 & Tr. We are interested in upper bounding the quantity
c Tr Dr exp 1+8
2xxp - D.). With-
with probability at least 1-
=! D,D,r
[Dya,lla | D. which is maximized when all entries that take value 1 in D, also take value 1 in Dyr This choice maximizes Tr D,D for any fixed values of Tr D3 and Dyy Let Tr D) = En. By Assumption 2, Dj and D,. must take different values (one is 0 and the other is 1) in at least rm diagonal positions. Combining this with our observation of the maximizing arrangement of ones and zeros in Dj and D', for this arrangement we have that Tr Dj Tr + 7n (g + 7)n and Tr D,Dy Tr D3 ξn. Since in the regime where Tr D2 ev and, by Assumption
19

---

Published as a conference paper at ICLR 2026
300
IHT
IMP
IHT
IMP
Figure 10. Standard deviation over the three random trials (top row) for each experiment reported in
Figure 5 (bottom row).
2, Dj and Dyr each have trace at least En, we have
Df S (1+8)V
|D 112 D1 -1--) -0
with probability at least - 2 exp -0000e E
exp -066en) This upper bound is increasing in
1+0
which can take value at most 1 - 7 since Tr Dj. 2 (E 7) and by construction D, S We therefore set ç = to yield the bound:
x.. DDDDy.
(1 + OV
[| 11 11
which holds with probability least exxp regime bound is independent of n while the first-regime bound decays as n 1/2 for large n the second-regime bound dominates. For all 8 (0,1) also have that 2 exp -c55
25 E =
4 4exp I+8
(088en ). Since this second-
4 4exp
-05ten
- 2exp BEve -eepp (-ob<em) Therefore, conclude that
146
x.. DjDgI,
(1+ 1+58V1-7
0 S Diara 112 Dirrrells
1
holds with probability at least - 2 exp -c5-eym tries of diagonal block submatrix As)i, Here include deterministic lower bound = DDD,
4 exp (=ebten) for all off-diagonal en-
2 which holds because DjD diagonal matrix with entries nonneg-
[D 11 P.y2. ative.
Next, consider block submatrix (A7 Ag), that is off the diagonal of A5 The entries of x! D.D,D, D.D.
(45 take the form |1D, P where D, and Dy may be the same or different. The 0, while the denominator terms have expectation Tr and / Tr Dyy, respectively. We bound this expression with high probability by applying - each term in the denominator, and asymmetric Hanson- -riigh the numerator:
numerator has Bia Dyape
|xr D,DDar
4
(9)
[.Darll 1
Tr DJ
t2 Tr D2 to
20

---

Published as a conference paper at ICLR 2026
with probability at least 2 exp
T7 Dgr eex TT D2 ) 2 2xxp Dy for
ct
-ct
-ct
universal constant We choose l2- = 8 8' Dj and 1a 8Tr Djr for d E (0, 1), to yield
| D2Dyy
| 22 D 112 (1 -1-8% (1-8) SAIDDIVTT Dr
with probability least exp
Tr Dr,,++ ) 2 exp
-ot
D2) 2exp (
(-06- Tv D..). By
1+6
1+5
Assumption 2, Tr Dj 2 en and Tr Dj. 2 En, so we have
|o] D2 Drrsvll
4
1 2 y 112
(1- (1-8)en
with probability at least depending whether Tr Dy. is less than or greater than en.
cti 2 exp Tt Dr Dy 4 exp
Now, we consider two regimes
1+0
In the regime where Tr Dy S En, we choose t 3n3/4, yielding:
|x. ,r.D,D,
1
D,silly 1Dy (1 8'm2/4
with probability at least exx (-cen//)) ) 4exp ( (-66ee)).
In the regime where Ir Dz > en, we choose 1 8n 111-n//// 1/4 Tr DjDy. Combining this with the implication of Assumption that D, S (1 -))n, we have
| DrDer
5n - 1/4 Tr D,
8(1-1) 7)
|D ll2 |Dy. 11
(1 (1-8)En
8)en1/4
with probability at least 1 - 2 exp
2/44)
4xx
For all d E (0, 1 (-ofee and for d S _ we
n2/446
2 exp
n/
-o55n
4exp ( (=rfen)
148
S 1
exp
-o5 3/4
1+nt/4
4 exp
T-
have
(1-6)n174 2 (1 (1-8)en174- Combining these, we have that
8(1-))
| D,D,yaal
S (1- 11-21/
| 2 1D
holds with probability at least 1 - 2 exp (=<82,,/44) 4 exp ( (=obee)) ) for all 5 E (0, ( and all entries of a block submatrix
x1/446
that off the diagonal of A As.
Now that we have high probability (and in some cases deterministic) upper and lower bounds on each entry of As. we combine them into high probability bounds on the eigenvalues of AS We can decompose AT.As B + C. where B is block diagonal and C is dense except for having zeros in block-diagonal entries. First, consider a single block submatrix B., on the diagonal of B. This block submatrix has diagonal values deterministically
and off-diagonal entries bounded deterministically from below by O and bounded above by (1+6)/I- with probability at least 1
1-5
-c8" m - 4 exp ( ==d"em). We use the variational definition of the minimum and maximum
2 exp
1+6
eigenvalues, and refer to B,, as B so that subscripts may denote indices within the block submatrix:
Amii (B) = min u Bu
=
= min
Brsuf +
Σ Bijuuy it
1+ min
Buyus
||
(6)
(1+ 8)
2
1-8 8
21

---

Published as a conference paper at ICLR 2026
where in (a) use the deterministic facts that diagonal entries of B take value 1 and that ||ul|2 = 1, and in (b) we use the high-probability upper bound on the magnitude of B and the observation that the minimum is achieved by vector TL 1 (this makes the cross terms most negative). By
similar line of reasoning, we can bound
max B) =1 max Σ Euyuuuu ||
S 1+ (1 + 6) -8-1),
1- 6
where in (a) the tigh-rrobbbility upper bound the magnitude of Ba,, the observation that the maximum is achieved by vector u and the requirement that the maximum dimension of B is 8 x 8. Both of these bounds hold with probability at least 1 - 2s(s 1) exp (=e8ec/") 4s(s 1) exp (=ed" ), by union bound over all off- entries in B. The spectrum of the full block matrix B is bounded by the minimum and maximum eigenvalues of its largest block, which can have size at most s X s. Thus the bounds above on Amia (B) and Amax (B) also apply to B.
Since A BCC,, it remains bound the spectrum C' and combine the results. The structure of C' is dense, with block diagonal submatrices of value zero and each other entry bounded between (=off offf
17-8) and (7- with probability least 1 exp
-6fn2/4 21/445
4 exp
1+6 ), for all
8 E (0, - For this matrix, we use coarse bound that ||C1lop S (1-8)177
which holds with
probability least 2s(s 1) exp following union bound.
no65 n//146 4s(s 1) exp ( (=obeem) for all (0, r-1))
1+0
Combining these spectral bounds on B and C' via Weyl's inequality, we have that
Amin (AāAs) 2 21- (1 11+0vII-7
8
1-8
8/n
and
Amnax (A5As) S1 (8 1)(1
1-
(1 1-8)n1/4
E (0, with probability at least 2s(s 1) exp (-oB'e/())
for
any
1-7
1+0
2s(8
1) exp
0082, n2/446
8s(s 1) exp
_08enn 1+4
Theorem (Hanson-Wrigh) Wright (Boucheron et al., 2013)). Let x be random vector with i.i.d. zero- mean 1-sub-Guurt Garb-iaassinn entries. Let H be square matrix. Then for a universal constant e
PP P Hr E|z"HH]] Hxl| 2 t s 2exp
(-I + ||
ct2
H is a diagonal matrix with all diagonal entries equal either zero or one, |H Tr H and H|l = 1. We also use an asymmetric version of Hanson-Wright derived as follows. Let
0 H 0
H
=
;
0
v
yielding
ct
P
P|lu" Ho E Hul
2
s exxp
| +
H 11
D PROOF OF THEOREM
Proof. The proof of Theorem 1 combines Lemma 1 with Theorem (Theorem in Jain et al. (2014)), which shows that IHT with an inflated sparsity level can recover sparse signal in linear inverse problem long the sensing matrix satisfies the restricted strong convexity and restricted
22

---

Published as a conference paper at ICLR 2026
strong smoothness properties with any positive finite parameters; i.e. with an arbitrary finite stricted condition number.
Theorem 3 shows directly that, under the conditions in the theorem statement, the objective value converges as f(w*)) - f(w") I E. This implics convergence of iterates due to restricted strong convexity and the fact that v((w") O:
a | w= w= f(w*) f(w")===))))))))r())) 2 (w - w", V)(w")) + a I!uK w²
2
which proves the additional result that || - w/ 1| 2c
We note that the proof in Jain et al. (2014) also implies that, IHT projects onto the smaller sparsity level rather than the inflated sparsity level 8, each step of IHT is still guaranteed to not increase the MSE loss [w):; the requirement that > 8 allows for strict objective decrease in each step. O
Theorem (Jain al. (2014)). Assume that the objective f has restricted strong convexity param- eter a and restricted strong smoothness parameter L at sparsity level 23 + s. with > 32 (,) )²² Assume that 0* arg min 51 f(0), ie. that the true signal is s-sparse. Then IHT with projec-
2 tion (hard thresholding) to sparsity level 8 and step size = 3L- run for K = O log [100)) iterations, achieves
f(O*) f(0)) e.
E PROOF THAT ASSUMPTION 2 FOLLOWS FROM ASSUMPTION 1 WITH HIGH PROBABILITY
Let D, diag (T(Xh; 2 0)) R=x", with (D:) as the set of all such distinet activation patterns possible under Assumption 1 with data X Raxd whose entries are drawn i.id. (0 1). Assumption 2 has the following two components:
1. Tr D; 2 en for all E [p], for some E (0, 1).
2. For all 1', the diagonals of D, and D. differ in at least yyn positions, for some y
(0, 1).
E. COMPONENT 1: LOWER BOUND ON TRACE OF ACTIVATION PATTERNS D;
Component follows from Lemma 2, which does not require Assumption 1.
Lemma 2 (based Ergen al. (2019)). Let S {i x. h > 0), where 3, are i.i.d standard
Gaussian vectors distributed as La). Then with probability at least 1- (=((((t=))-(())) inf, [S] > ne. Here E E (0, 1), yp fixed numerical constant satisfying /8p > 0, satisfies n( V8p) > d, and H is the binary entropy function.
Proof. Consider the symmetric event E := :=suphae {i 27 0H 2 n(1- - =). Then
P <
Σ
Plh Sst. af 0,
VCMm (122(11-))
n
(m)(1=))
n(1 - 5)
Se (x H((e)
in which the second inequality follows from the Kinematic Formula (by flipping the sign of h).
Theorem 4 (Kinematic Formula (Amelunxen et al., 2014)). Let X be n x d i..d Gaussian matrix and = 2 with any E O. If satisfies n( - /Bp) 2 d, we have
P P[h 40 s.I. Gh > 20
3h + S.s. x 2 of pm
23

---

the indicator Published as a conference paper at ICLR 2026
Remark 3. If we further assume Assumption 1, specifically that |1hllo S s; S k. we can righten with probability at
Lemma 2 as follows. If | k S d, and set up 128 then Tr
H(E)) as long as 2 4k. Following a union bound over activation patterns, Tr D; 2 for all i E [p] with probability at least 1 pe ( ((((-=---))) -H(s)) long as n 2 4k.
least - ---
E.2 COMPONENT 2: HAMMING SEPARATION OF ACTIVATION PATTERNS
Definition (8-isometric embedding (Plan & Vershynin, 2014)). A map f X + Y is a 8- 8-isometry between metric space X with distance metric dx and metric space Y with distance metric dy for all x, x' E X, |dy (zx), (xx)) - dx (x, x') S 6.
Theorem 5 (Hamming embedding, Theorem 1.5 in Plan & Vershynin (2014)). Consider subset K S. 5 S and let > Let X be an n x d random matrix with independent ( 1) entries. Lef 2 C8- where w(K) EsupxeK Kg,=)| is the Gaussian mean width of K, with (N0,, Ta). Then with probability at least 1 - 2 (exp(_cō21), the sign map f(x) = sign(Xa) f: K 1- is 5-isometric embedding between K S with normalized geodesic dis- lance metric da(a, x.) =-s cos- x') and {-1,1]" with normalized Hamming distance metric (f(f(z) + Ei= # f(z'), Here C and denote positive absolute constants.
apply to unnormalized generator vectors from K. S. Rd and indicator-base rather than sign-based activation pattern embedding. Consider a subset K IRd and let 5 O. Let X be
Corollary 1. Theorem 5 may be restated
d random matrix with independent N entries. Let n 2 C8 (K)) where w K)) == Esupre Kig, [|=][," > is the normalized Gaussian mean width of K, with
NO0,, Ta). Then with probability at least 1 - 2exp(=c8 map f(x) = 1(Xx 2 0}, f: K + (0. is 6-isometric embedding between K S Bed with normalized geodesic distance metric do(x, ) cos-1
TT. |T and {0, 13"h with normalized 1 E f(z), # f(x'), Here and e denote
Hamming distance metric du dn(f(x), positive absolute constants.
Corollary ensures that a set of generator vectors {he]_= that are sufficiently separated in nor- malized geodesic distance will yield activation patterns D, diag (1(Xh; 2 03) whose diago- nals are separated in normalized Hamming distance, with high probability for i..d. Gaussian data X Ruxd Specifically, for the diagonals of D; and D, differ in at least yn positions for all i A i' with probability at least 2 exp(-cc52m we require a set of generator vectors (h: V that (1) include the planted first layer weights u,, and (2) are separated by at least 7 + 6 in normalized geodesic distance. In Appendix E.3 we show that both of these properties hold with high probability both of the sparse weight conditions in Assumption (recall that ot are the first layer weights and of are the second layer weights, which are fused during IHT):
(a) uf 1-1,0,113 ||u llo of I Vi [p] and kp 8, or (b) ug E Rdd. || Ilu = 8; E [smin 1c), of E f- E and E 8; S holds.
E.3 SAMPLING SPARSE ARRANGEMENTS
E.3.1 REAL-VALUED PLANTED NEURONS
We now show that random sampling of hyperplane arrangements can be guaranteed to contain the planted activation patterns, while simultaneously ensuring a packing of the Euclidean sphere Rdd
Theorem 6. Let X R have i...d. (O 1) entries. Fix m unknown vectors Wi,..., . Wm Rd each with I!lillo 8; E [Smin, kl, and an error tolerance < ē < 1. Note that the choice of affects the permissible sparsity range [Smins kt), Set
log(2 (Lox(2m)) *
12
log
where e. > is an absolute constant. Consider the set of all supports with |S] E [smins k), and draw T supports from this set uniformly random. For each randomly drawn support S, draw |S|
24

---

Published as a conference paper at ICLR 2026
values i.i.d. from (0, 1) and embed these in Rd by setting entries in S to their random Gaussian values and zero- prrddpadding outside S. Record the two collections
r r= 2 :
G (
where h denotes the zero pro-paaded generator and h, is its normalized version. There exists 8 such that, with probability at least 1 - E over the draws of X and all T generators. the following hold simultaneously
1. Coverage: [Xw, 2 0] E for every i m..
2. Minimum geodesic separation: For all distinct g. g' in G. 11g g'll2 2 5. This Euclidean
separation of unit vectors implies geodesic separation: de a 0.69 82.
1 cos cos- (g" 2
Proof. The strategy used to prove coverage to show that the cones {u sign( xu) (Xh) are not too narrow, for Gaussian i.i.d. training data X E RRad and a fixed vector E IRd. Specifi- cally, a bound on the cone sharpness developed in Kim & Pilanci (2024) implies that the probability that a uniformly sampled vector on the sphere falls into this cone is at least O O(log )))) We then apply this result to n X 8; submatrices of X translate to sparse generators, and control the error probability via the union bound. We first reintroduce the notion of cone sharpness:
Cone sharpness. For any support and non-zero u E IK;;, set Duu) := diag(I[Xsu D 0)) and define the cone (u) := { c E R : (2D(u)- Xg 2 0}. By the cone cone- -harpeess bound of Kim & Pilanci (2024) there are universal constants c, C1 > 0 such that
X (Kss((u),
v [
& C. 2 8s., C. = + 200c/ c2og(2n), 8z, := n + c
where the sharpness C(C,,) =) of a cone K. with respect to a fixed unit vector = is defined as 2'K,2) reK,,, a-v=s ||u|| + l|lll2 following Kim & Pilanci (2024). Let & be the high probability event that the cone sharpness for each the m planted neurons at most C,; & occurs with probability least 1 m(n 10 c omiin))). Now relate sharpness to the probability of sampling a specific pattern.
Spherical cap inclusion. Fix S and + and write == u/llul22- On E there exist a Ks(uu) with a - = = and ||a||2 + |lb]|2 & C.. Setting = (a + b)/2 yields (q. 2 1/(2..). Therefore the spherical cap
C₂ (u [ (v,2) 2 1/(2C-))
contained in (u)
Cap measure. For h ~ (0, Is,) the direction h/ |1 is uniform on S Standard sur- mice-measure estimates give
Ps. = P[T[X hh 2 0] = [[gu 2 o] 2 Surf Surf,,-I(C,)) (log(2n))*
Coverage probability. As consequence of the above inequality, the arrangement pattern of each planted neuron is sampled with probability least (log(2r))
(og((n) After T samples,
the probability that we have not yet sampled all m planted neuron activation patterns is at most
m 1
after T = log(2n) &
log 2m
random draws we are guaranteed to sample
log((n))
all m planted patterns with probability at least 1 2.
Packing of generators. Consider any two generator vectors h hi with supports S,S', respectively. We have
Eiesnsr hugh
WTW =
V Eies hi Eiegr h?
25

---

Published as a conference paper at ICLR 2026
where h, and i are iid. distributed (N 1). Using union bound asymmetric Hanson- Wright (Theorem 2) in the numerator and symmetric Wannonn-riigh (twice) in the denominator, have
4
ct
1_c (5)+
ctt
S exxp
+2exp
+2exp
vIsI + // 181
+
Sn +11
15
We choose l2 - 7/5), 13
715% yielding
P h h
(1 (1+
4
exx
S' +
cli
+ 2exp
(17
2exp
-c781s91
1
Next, we choose t = -11
82
(> + 1)
S||SS |
to yield
82
c(1 ) (r 1)
1 2
< 2xxp
2
S + (1- 52 ) + 1-1/SSI
(=07915) + 2 exp
(=
- 2exp
(1+>
Assuming that all planted neurons, and thus all generator vectors we need to consider, have sparsity level 8; E [smin: k:], we have
8² P hTT 2
c(1 (11-))2 (1 3 min
( (=07² smia
< exx
+ 4exp
2
k + (1
82 (7 + 1)A:
11+7
2
For sufficiently large 8min a union bound over all pairs of generators allows us to bound P [hTh'] 2 1 - 2
m(n-10 + suii) uniformly over all pairs h,h,, as required for an overall failure probability at most ē. Finally, the Fuclidean packing follows as
|/ M12 ||A112 + |M11 2 21 M > 2 (2 82) 82.
Euclidean 8-packing implies the stated separation in normalized geodesic distance follows:
g g'
do(g.g))
cos-1 cos π
Tgll2 9
(a) 1.38
(1
of g')
(b) 0.69 82_
where in (a) we use the fact that cos
-0ss-1(1
(1 x) > 1.38x for E |0, 1] and in (b) we use that
||le| ||g'll and ||g g'1l ||g1| 1|g'112
2g T g' = -2g g 82.
E.3.2 DISCRETE-VALOU PLANTED NEURONS
Theorem 7. Fix integers and < d. For subset S C [d] with entries S, and |S| = k, define the generator set
g(o)s, - Jj g(0)) 0ff s}. Gs {o(0) (-1,0011 o
i.e. all possible sign assignments on the coordinates in S and zeros elsewhere (19s) = 2* Given matrix X xdd form the associated arrangement list
Ts = {1Xg20] g c 10.
26

---

Published as a conference paper at ICLR 2026
(i) Coverage. Let wi...., tm (f1,0,,1] be k-sparse vectors (each with exactly k non- zeros). Then for every E (!....,m) 1[x 2
s.)p(w.))
iii Minimum geodesic separation. For any two distinct k-sparse vectors u, E {-1,0, 133,
1 dc(u,v) cos
u
0.69
||u/12 11112
xk
Proof. (i) Coverage. Fix (1,...,m.) m,] and set S supp(w:). Because wi has entries #1 S and zeros elsewhere, there exists E {-1, 131s1 such that w, g Gs- Hence the pattern 1[w; 2 0] belongs to s's-
(ii) Separation. Let be k-sparse vectors in {-1, 0, 13 There is an index j with u, # Uj. so Juj Vy} 2 Therefore ||u - v Σi= (u: - vz)² (u, - v,)² 2 1, implying ||u - w112 2 1.
From Euclidean separation we can infer normalized geodesic separation as follows:
1
u v || ll |1u11
da(u, v)
cos
π
(e)
1.38
uTu
||u||2 11:112
IV
1.38
||ua12 2 + ||
2 || |10|2
0.69
æk
where in (a) we use the fact that cos -1 (1 - x) 2 1.38x for x E [0, 1]. in (b) we use that ||u- 212
|ull2 + ||=12 2uTrr 2 and in (c) that ||u||2 || 7 since both and have exactly
nonzero entries each with magnitude one.
27


Reviews:
[{"review_id": "14rx7j5T0L", "reviewer": "Reviewer_RfMf", "summary": "This paper provides theoretical guarantees for recovering the sparse weights of a ReLU network using an efficient, memory-linear algorithm. The proposed method, validated on tasks on MNIST, matches or exceeds the performance of a memory-inefficient pruning baseline.", "questions": "please see weaknesses.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This work establishes the sparse recovery result for ReLU MLPs. \n\n2. Using shallow networks on Gaussian data, the experiments demonstrate that Iterative Hard Thresholding (IHT) is more effective and memory-efficient than a strong IMP baseline, recovering superior sparse networks.", "weaknesses": "Using the eq. (2), training the 2-layer MLP with MSE loss is a nonconvex optimization problem.  From the lines 162-175, the authors consider fixed generator vectors $h_i\\in R_d$ and fusing the weights via $w_i=u_iv_i$ to build the eq. (3). This transformation allows the results to be readily obtained using sparse recovery theory. However, several issues arise:\n\n1. The theory requires the matrix $A$ to satisfy the conditions of Lemma 1, which is attributed to the sparsity of $h_i$ or $u_i$. Typically, \n $X$ is assumed to be a Gaussian data matrix. How is Lemma 1 satisfied when the input data is non-Gaussian, such as the binary-valued pixels in the MNIST dataset\n\n2. The theoretical results are promising, but how does the method perform on more complex, real-world datasets like CIFAR-10 or ImageNet-200?\n\n3. This work focuses on recovering the sparsity pattern of the weights $u_i$. What can be said about the sparsity or recovery guarantees for the corresponding weights $v_i$?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ngVaQGKM2s", "reviewer": "Reviewer_4D7G", "summary": "This paper provides the first recovery guarantee for sparse ReLU neural networks by reformulating sparse MLP training as a structured compressed-sensing problem. Under Gaussian data and enumerated activation patterns, the sensing matrix satisfies restricted strong convexity and smoothness, allowing Iterative Hard Thresholding (IHT) to exactly recover sparse weights with high probability. The method requires memory linear in the number of nonzeros. Experiments on planted MLPs, MNIST, and implicit neural representations validate the theory and show that IHT performs comparably or better than Iterative Magnitude Pruning (IMP) while being significantly more memory-efficient.", "questions": "1.Can you analyze how A evolves during sequential convex updates—for example, by measuring similarity between consecutive A matrices—and test how update frequency affects convergence?\n\n2.Can you report empirical ($\\alpha,\\beta,\\beta/\\alpha$) for at least one representative setting to show whether the RSC/RSS conditions approximately hold in practice.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "A rigorous and elegant theoretical framework connecting sparse neural network recovery with compressed sensing, supported by solid experimental validation.\n1. First work proving exact recovery (Theorem 1) for sparse MLP weights with explicit high-probability bounds and convergence rates; recovery builds on Assumption 2 and the RSC/RSS properties (Lemma 1).\n2. Convex reparameterization (Equation 3) via weight fusion turns the nonconvex problem into a linear sensing framework; Assumptions 1–2 are well-motivated and verified through constructive examples in Appendix D.\n3. Memory scales linearly with sparsity (s); sample complexity grows with active weights rather than total parameters. IHT achieves competitive or superior accuracy to IMP while using less memory.\n4. Well-structured from motivation through theory to experiments, with a comprehensive appendix providing implementation and proofs.", "weaknesses": "The main limitation lies in the unquantified theory–practice gap and missing diagnostics connecting empirical behavior to theoretical parameters.\n1. Theorem 1 assumes a fixed, enumerated sensing matrix A, but experiments adopt sequential convex updates—updating A after each IHT iteration (following an initial fixed epoch). The paper acknowledges this deviation (Appendix A) but provides no convergence analysis or ablation on update frequency.\n2. Empirical values of the RSC and RSS constants and their ratio are not reported, leaving unclear how tight the theoretical conditions are in practice.\n3. Results are averages over three runs without error bars, standard deviations, or variance. Adding uncertainty plots (e.g., box plots or success-probability curves) would strengthen performance claims.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "1iXNelz30y", "reviewer": "Reviewer_XpEf", "summary": "This paper introduces and studies a method for training sparse shallow ReLU neural networks. The method reformulates the neural network optimization problem as a convex optimization problem and uses iterative hard-thresholding to obtain a sparse solution. Theoretical results guarantee the recovery under reasonable assumptions in the case of shallow networks. Empirical experiments verify that the method also works well for deep networks.", "questions": "Can the theory be extended beyond Gaussian training data?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The method for training sparse networks is innovative. The theoretical guarantees are interesting and relevant. The practical performance of the methods is good. The presentation of the paper is good, it is easy to read and clear.", "weaknesses": "The theory only applies to shallow networks. It would be interesting to try to extend it to deeper networks. Nonetheless, this is a valuable theoretical contribution.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "U9Sjc1nqU0", "reviewer": "Reviewer_A2Ki", "summary": "This paper presents the theoretical recovery guarantees for sparse ReLU neural network weights, focusing on two-layer scalar-output MLP networks. The authors prove that under certain structural conditions on sparse network weights and random Gaussian training data, an Iterative Hard Thresholding algorithm can exactly recover these weights with memory linear in the number of nonzero weights. The theoretical analysis uses convex reformulations of MLPs and establishes restricted strong convexity and restricted smoothness properties. Experimental validation on planted MLPs, MNIST classification, and implicit neural representations demonstrates competitive or superior performance compared to iterative magnitude pruning.", "questions": "- How can the proposed analysis be extended to non-Gaussian data distributions?\n- More analysis on how performance degrades when assumptions are violated would strengthen the work.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "- This work provides formal recovery guarantees for sparse neural network weights, filling an important gap between compressed sensing theory and neural network optimization.\n- The paper clearly articulates assumptions and provides detailed proofs showing how these conditions arise with high probability under Gaussian data, making the theoretical results verifiable and interpretable.\n- The experimental studies, which extend beyond the theoretical analysis, demonstrate the broader applicability of the guarantees, including vector outputs, deeper networks, and various tasks.", "weaknesses": "- The theoretical guarantees, although interesting and important, are restrictive. It is limited to two-layer, scalar-output networks with Gaussian random data. However, real-world applications may involve deeper networks, structured data distributions, and multi-dimensional outputs.\n- Assumption 1 requires concrete weight structures, i.e., binary hidden weights or binary output weights, which may not reflect realistic sparse networks.\n- The paper only compares against IMP. More recent sparse training methods, such as dynamic sparse training, pruning at initialization variants, and other memory-efficient approaches, can be compared. Another limitation is that the runtime analysis, with IHT sometimes being slower than IMP for larger problems, which may be a limitation in practice.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "wpc7VuNno6", "reviewer": "Reviewer_T4on", "summary": "The paper studies sparse training through a recovery-guarantee lens for two-layer ReLU MLPs with scalar output. Using a convex reformulation that enumerates activation patterns and fuses layer weights, the training objective becomes a structured linear sensing problem. Under Gaussian inputs and two separability conditions on activation patterns, the sensing matrix satisfies restricted strong convexity and smoothness, which makes the planted sparse weights uniquely identifiable. Experiments on planted MLP recovery, MNIST classification, and implicit neural representations show the proposed method is competitive with or better than iterative magnitude pruning while using far less memory, even beyond the exact theory regime.", "questions": "1. Extension beyond Gaussian inputs? Can the proof be adapted using sub-Gaussian or whitened real data, and what additional conditions would be needed on X?\n\n2. The experiments use randomly sampled patterns and sequential convex updates. Is there a provable guarantee when A is updated during training, perhaps via a stability or tracking argument?\n\n3. The appendix describes vector-output handling and layer-wise training. Could the convex deep-net formulations be combined with your analysis to yield recovery guarantees for deeper or multi-output networks?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1.The paper gives identifiability and efficient recovery for planted sparse weights in two-layer scalar-output MLPs under Gaussian data.\n2. The activation-pattern convexification connects ReLU training to linear sensing in a way that enables analysis with classical tools.\n3. Across planted recovery, MNIST, and image INRs, IHT is typically stronger than IMP while using memory that scales with s rather than with dense parameter counts.", "weaknesses": "1. Scope limited to shallow scalar-output MLPs in theory.\n2. The main theorem relies on i.i.d. Gaussian covariates and full enumeration of activation patterns, which can be unrealistic or computationally heavy.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)