On the Limits of Sparse Linear Model with Shuffled Differential Privacy

Preprint

O N THE L IMITS OF S PARSE L INEAR M ODEL WITH S HUFFLED D IFFERENTIAL P RIVACY Yan Zhang1 Jiong Tuo2 Meng Ding3 Yuheng Ma4 Di Wang5∗ Nankai University 2 Xi’an Jiaotong University 3 University of Massachusetts Boston 4 East China Normal University 5 King Abdullah University of Science and Technology (KAUST)

arXiv:2610.32822v1 [cs.LG] 23 Oct 2026

1

A BSTRACT Sparse linear regression is a fundamental problem in high-dimensional statistics, and its private counterparts leave a large gap between the central and the local model of differential privacy. Shuffled differential privacy (SDP) sits between the two, offering accuracy closer to the central model without the trust requirement of a central curator. This paper studies sparse linear regression under both non-interactive and interactive SDP protocols. In the non-interactive setting, any non-interactive (ϵ, δ)-SDP protocol is shown to incur an ℓ2 -norm estimation er√ dk ) with dimension d, sparsity k, and n users, revealing an inherent ror of Ω( nϵ dimensionality barrier. In addition, a non-interactive algorithm is proposed that q √ d k k achieves an ℓ2 -norm upper bound of Õ( n + nϵ ). To overcome this barrier, we turn to the interactive user-level setting, where each user holds m samples. A two-round interactive algorithm is developed that first recovers the support through aq heavy hitter step and then estimates the parameter, achieving an erk ror of Õ( nm + n√kmϵ ) with no polynomial dependence on d. Furthermore, by q extending the score attack framework to user-level data, a new lower bound √ Ω( k log(d/k) + k log(d/k) ) is established for any user-level (ϵ, δ)-DP, showing nm n mϵ that the proposed algorithm is near-optimal with only a gap of the factor log(d/k).

1

I NTRODUCTION

Protecting data privacy has become a paramount concern in modern information and database systems, which often store highly sensitive personal records, and (sparse) linear models, a cornerstone of both statistics and machine learning, are routinely fitted to such data. Differential privacy (DP) has emerged as the gold standard for privacy-preserving data analysis (Dwork et al., 2006; Dwork & Roth, 2014), providing mathematically verifiable protection that remains resilient to attackers with arbitrary auxiliary information. The central model of DP assumes the existence of a trusted curator who collects raw data from all individuals and then applies privacy-preserving mechanisms before releasing the results. However, in many distributed settings, users may not trust a central entity with their unencrypted sensitive information. To address this, local differential privacy (LDP) was introduced, shifting the privacypreserving step to the user’s side; each individual perturbs their own data locally before sending it to a server (Duchi et al., 2013; 2014; Kasiviswanathan et al., 2011). Despite this stronger privacy guarantee, a critical issue with LDP is its limited utility, as the significant amount of noise required to protect each individual entry often leads to much lower statistical accuracy compared to the central model. Shuffled differential privacy (SDP) has recently emerged as a promising middle ground that simultaneously addresses these concerns (Cheu et al., 2019; Erlingsson et al., 2019; Balle et al., 2019). By introducing an intermediate shuffler that de-identifies and permutes the local messages ∗

Corresponding author.

1

Preprint

before they reach the server, SDP can provide accuracy levels approaching those of the central model while maintaining a distributed architecture that does not require a trusted curator. Private linear regression has been studied both as stochastic optimization, which minimizes the empirical or population risk, and as statistical estimation, which recovers the parameter θ∗ ∈ Rd of the linear model y = ⟨x, θ∗ ⟩ + ζ under the ℓ2 -norm error (Kifer et al., 2012; Talwar et al., 2015; Cai et al., 2021). We take the estimation view and focus on the high-dimensional sparse model, where θ∗ is k-sparse with k ≪ d and d may far exceed the sample size n (Raskutti et al., 2011; Yang et al., 2014; Jain et al., 2014); here the main difficulty under privacy is to avoid an error that scales polynomially with d. In recent years, there have been many studies of the DP sparse linear model under different privacy notations. Specifically, in the central DP model, the l2 -norm estimation error for k-sparse parameters q

d d scales logarithmically with the dimension d, achieving a nearly optimal rate of O( k log + k log n nϵ ) (Cai et al., 2021). In sharp contrast, the local DP model imposes a much heavier privacy cost, where the estimation error inevitably suffers from a polynomial dependency on d; as summarized in Table 1, the improved analysis of sparse local-private regression shows that non-interactive and interactive LDP algorithms still obey strong dimension-dependent lower bounds (Zhu et al., 2024). This creates a massive gap that makes LDP algorithms essentially ill-suited for modern high-dimensional scenarios where d ≫ n. However, on the other side, for SDP, recent studies show that for many problems, such as histogram estimation Balcer & Cheu (2019), their theoretical behaviors are more similar to the central DP, while for other problems, such as parity function learning Cheu & Ullman (2021), the theoretical behaviors are more similar to local DP. Given the large gap on the limitations of sparse linear model between the two models, we have a natural question: does sparse linear regression in the shuffled model bridge the gap toward the efficiency of central DP, or does it still reflect the dimensionality barriers of local DP?

Table 1: Comparison summary of our work with related studies on (ϵ, δ)-DP linear regression under sub-Gaussian data assumptions. Here n, m, k, and d denote the number of users, samples per user, sparsity, and dimension. Note that we can regard m = 1 when the privacy unit is sample. All bounds are measured by the ℓ2 -norm estimation error. The Õ notation suppresses polylogarithmic factors. The full table can be seen in the Appendix A. Model Type Central -

Unit sample sample

Non-interactive sample Local

Setting Upper Bound q d general Õ( nd + nϵ ) q k k-sparse Õ( nk + nϵ ) √ k-sparse Õ( d√nϵk ) √

Lower Bound q d Ω( nd + nϵ ) q d d Ω( k log + k log ) n nϵ q log d Ω( dknϵ ) 2 q dk Ω( nϵ 2)

Source Cai et al. (2021) Cai et al. (2021) Zhu et al. (2024)

k-sparse Õ( k√nϵd ) Zhu et al. (2024) Interactive q k k k user-level k-sparse Õ( nm + √nmϵ ) Ω( √nmϵ ) Ma et al. (2024) q √ √ dk Non-interactive sample k-sparse Õ( nk + dnϵk ) Ω( nϵ ) Shuffle Our Work q q k log(d/k) k log(d/k) k k √ √ Interactive user-level k-sparse Õ( nm + n mϵ ) Ω( + n mϵ ) nm sample

In this work, we provide a rigorous view of sparse linear regression under SDP. Our contributions are three-fold. • We first study the non-interactive protocol. By converting robust (ϵ, δ)-SDP to pan-privacy, √ dk we establish a lower bound of Ω( nϵ ) on the estimation error measured by the ℓ2 -norm difference, revealing an inherent dimensionality barrier, which is more similarqto LDP. We √ also propose an (ϵ, δ)-SDP algorithm, which achieves an upper bound of Õ( nk + dnϵk ), characterizing the limits of non-interactive SDP. 2

Preprint

• We then study how we can overcome the theoretical barrier with relaxed settings. Specifically, we consider the interactive and user-level SDP protocol, where each user has m samples. We develop a two-round interactive algorithm, U2R (User-level Two-Round Regresq k sion), that achieves an upper bound of Õ( nm + n√kmϵ ). Moreover, we develop a novel q  k log(d/k) k log(d/k) √ user-level score attack method and show a lower bound of Ω + for nm n mϵ user-level (ϵ, δ)-DP. The two rates differ only by a factor of log(d/k), indicating that our upper bound is near-optimal and that, with interaction and user-level data, the problem is more closely related to central DP. • We evaluate U2R on synthetic and real data. Across the privacy budgets ϵ, dimension d, and user counts n, it has the lowest coefficient error and MSE among the private methods.

2

R ELATED W ORK

Private high-dimensional regression has been studied across three main privacy models. We provide a comprehensive review in Appendix A. In the central model, sparse estimation typically incurs only a logarithmic dependence on the dimension d (Talwar et al., 2015; Cai et al., 2021), and the score attack of Cai et al. (2023) gives the lower-bound technique that we extend to the user-level regime. The local model is more restrictive, with both non-interactive and interactive estimators exhibiting polynomial dependence on d (Wang & Xu, 2020; Zhu et al., 2024); user-level multi-sample regimes alleviate these costs (Ma et al., 2024). Positioned between the two, the shuffled model leverages privacy amplification (Cheu et al., 2019; Erlingsson et al., 2019; Balle et al., 2019), while multi-message protocols compute a sum with central-model accuracy (Balle et al., 2020), which is the summation step both our algorithms use. On the negative side, Cheu & Ullman (2021) relate shuffled DP to pan-privacy, the route our non-interactive lower bound follows. Shuffled protocols have been studied for primitives such as heavy hitters (Bun et al., 2019), but their application to sparse linear regression remains largely open, which is the gap we fill.

3

P RELIMINARIES

This section introduces the notations, the models, and the problem setting used throughout the paper. Notations. For an integer d, let [d] = {1, . . . , d} and let Rd be the d-dimensional Euclidean space. For a vector v ∈ Rd , ∥v∥p is the ℓp -norm, and ∥v∥0 is the number of non-zero entries. The support of θ is supp(θ) = {j ∈ [d] : θj ̸= 0}, and for a set S ⊆ [d] we write θS ∈ R|S| for the restriction of θ to the coordinates in S. Given a matrix X ∈ Rn×d , let x⊤ i be its i-th row and xij (or [X]ij ) be its (i, j)-th entry, which is also the j-th element of xi . We use ∥X∥2 for the spectral norm and ∥X∥∞,∞ = maxi,j |xij | for the max norm, and λmin (X) for the minimal eigenvalue of X. For an event A, I[A] denotes the indicator, i.e., I[A] = 1 if A occurs and I[A] = 0 otherwise. A zero-mean random variable X is sub-Gaussian if its sub-Gaussian norm ∥X∥ψ2 = inf{c > 0 : E[exp(X 2 /c2 )] ≤ 2} is bounded; a random vector x is sub-Gaussian if ⟨x, u⟩ is sub-Gaussian for e hides logarithmic factors. every unit vector u. We use O(·) and Ω(·) in the standard sense, and O(·) 3.1

D IFFERENTIAL P RIVACY

Definition 1 (Differential Privacy (Dwork et al., 2006)). A randomized algorithm M with domain Dn is (ϵ, δ)-differentially private (DP) if for all measurable sets S ⊆ Range(M) and for any two neighboring datasets D, D′ ∈ Dn that differ in exactly one record, it holds that P[M(D) ∈ S] ≤ eϵ · P[M(D′ ) ∈ S] + δ.

(1)

Definition 1 is stated with respect to a neighboring relation on Dn . In the sample-level regime a record is a single pair (xi , yi ), so D and D′ differ in one sample. In the user-level regime a record is the whole local dataset Ui of user i, so D and D′ differ in all m samples held by one user. All our user-level theorems use the latter relation. 3

Preprint

The central model assumes the curator sees the raw data. The shuffled model removes this assumption, and instead places a trusted shuffler between the users and the data collector. Definition 2 (Local Randomizer, Shuffler, and Analyzer). A protocol in the shuffled model is a triple P = (R, S, A), where the local randomizer R : D → Y t maps one record to t messages, the shuffler S : Y nt → Y nt applies a uniformly random permutation to the messages it receives, and the analyzer A : Y nt → O maps the shuffled messages to the output of the protocol. On a dataset D = (z1 , . . . , zn ), the view of the data collector is
VP (D) = S R(z1 ), . . . , R(zn ) , (2) and the output of the protocol is MP (D) = A(VP (D)). In the user-level regime, the record of user i is U
i ∈ Dm , so the local randomizer is R : Dm → Y t and the view is VP (D) = S R(U1 ), . . . , R(Un ) . Definition 3 (Shuffled Differential Privacy (Cheu et al., 2019; Balle et al., 2019)). A protocol P = (R, S, A) satisfies (ϵ, δ)-shuffled DP (SDP) if the view VP satisfies Definition 1, i.e., for any two neighboring datasets D, D′ and any measurable event E ⊆ Y nt , P[VP (D) ∈ E] ≤ eϵ · P[VP (D′ ) ∈ E] + δ.

(3)

Definition 3 implies that MP is (ϵ, δ)-DP in the sense of Definition 1. Definition 4 (Non-interactive and Interactive protocol). A protocol is non-interactive if every user applies R to their own record once and sends the resulting messages in a single round, so that R does not depend on any other user’s messages. A protocol is interactive if it proceeds in T > 1 rounds, where the randomizer used by a user in round t may depend on the aggregated output released by the analyzer in rounds 1, . . . , t − 1. 3.2

P ROBLEM S ETUP

We consider the setting of sparse linear regression. Let D = X × Y ⊆ Rd × R be the data universe and let there be n users in the population. We distinguish two data regimes. In the sample-level regime, each user i ∈ [n] holds one pair (xi , yi ) ∈ X × Y. In the user-level regime, each user i holds m i.i.d. samples Ui = {(xi,j , yi,j )}m j=1 , which we also write in matrix form as (Xi , yi ) ∈ m×d m R × R ; the total number of samples is then N = nm. Assumption 1 (Sparse Linear Model). The samples are i.i.d. realizations of the sparse linear regression model y = ⟨θ∗ , x⟩ + ζ, where the distribution of x has mean zero, the noise ζ satisfies E[ζ | x] = 0 and is sub-Gaussian with variance proxy σ 2 , and the underlying parameter lies in Θk = {θ ∈ Rd : ∥θ∥0 ≤ k} with k ≤ d and ∥θ∗ ∥2 = O(1). In the user-level regime, the samples are i.i.d. both within and across users. Assumption 2 (Covariance Structure). Let Σ = E[xx⊤ ] be the population covariance matrix. There is a constant κx > 0 such that ∥Σ−1/2 x∥ψ2 ≤ κx .

4

U NDERSTANDING THE L IMITATIONS IN N ON - INTERACTIVE SDP

In this section, we characterize sparse linear regression under non-interactive shuffled DP. We first establish a minimax lower bound on the ℓ2 -norm estimation error, and then give an algorithm whose bound is decent. e k = Θk ∩ {θ ∈ Rd : ∥θ∥1 ≤ 1}. The ℓ1 normalThroughout this section the parameter space is Θ ization strengthens the ℓ2 bound of Assumption 1 and is what keeps the sensitivity of a single record under control once each user privatizes their own statistics; it is not needed in Section 5. We also use the following conditions, which are used in this section only. Assumption 3 (Curvature and Boundedness). There is a constant κ∞ > 0 such that ∥Σw∥∞ ≥ κ∞ ∥w∥∞ for any w ̸= 0, and there are constants cx , cy = O(1) such that ∥x∥∞ ≤ cx and |y| ≤ cy almost surely. The first condition is the ℓ∞ curvature used by the soft-thresholding step, and the second lets each record be clipped at a known scale before it is privatized. 4

Preprint

4.1

L OWER B OUND

We begin with the fundamental limits of the non-interactive setting. The standard route to a minimax lower bound under local DP does not apply here. In the local model each message depends on a single record, so the KL divergence tensorizes, DKL (P n ∥U n ) = n · DKL (P ∥U ), and private Fano’s inequality or a direct mutual information bound can be used. In the shuffled model the messages of the n users are permuted before the analyzer sees them, and the mutual information between the input and the view VP can be unbounded, so these arguments do not carry over. To solve the problem, we introduce pan-privacy, which constrains the internal state of an algorithm rather than the messages of individual users. Definition 5 (Pan-Privacy (Cheu & Ullman, 2021)). An algorithm M processes a stream Z = (z1 , . . . , zn ) one record at a time: it updates an internal state from zi and the previous state, and computes an output from the final state. Write MI (Z≤t ) for the internal state after the first t records and MO (MI (Z)) for the output. Then M is (ϵ, δ)-pan-private if for every pair of streams Z, Z ′ differing in one record, every t ∈ [n], and every measurable set T of state-output pairs, 
 
 ′ P MI (Z≤t ), MO (MI (Z)) ∈ T ≤ eϵ · P MI (Z≤t ), MO (MI (Z ′ )) ∈ T + δ. (4) Pan-privacy protects the record against an adversary that observes the output together with the internal state at one point in the stream. Cheu & Ullman (2021) showed that any algorithm satisfying robust (ϵ, δ)-shuffled DP can be simulated by a (ϵ, δ)-pan-private algorithm on a constant fraction of the users; here, robustness means that the protocol remains (ϵ, δ)-shuffled DP as long as a constant fraction of the users follow it, a condition met by all protocols considered in this work. Therefore, to obtain the utility lower bound of SDP algorithms, it is enough to bound the error of pan-private algorithms. The lower bound is proved on the following family of hard instances, which specifies a distribution over the sparse parameter θ∗ and over the data (x, y). √ Definition 6 (Hard Distribution). Let ν be the target ℓ2 error and let γ = ν/ k. Define θZ = γZ with Z ∈ {+1, 0, −1}d , where each coordinate Zj is drawn independently according to P{Zj = +1} =

k , 4d

P{Zj = −1} =

k , 4d

P{Zj = 0} = 1 −

k . 2d

(5)

Given θZ , a sample (x, y) ∈ {−1, 1}d × {−1, 1} is generated by drawing x uniformly from the hypercube, so that P(x) = 2−d , and setting y = ⟨x, θZ ⟩ + ζ with P{ζ = 1 − ⟨x, θZ ⟩} =

1 + ⟨x, θZ ⟩ , 2

P{ζ = −1 − ⟨x, θZ ⟩} =

1 − ⟨x, θZ ⟩ . 2

(6)

The expected sparsity of √ θZ is k/2, and Markov’s inequality gives ∥θZ ∥0 ≤ k with probability at least 1/2; under ν ≤ 1/ k we also have ∥θZ ∥1 ≤ 1. The construction gives y ∈ {−1, 1} and E[ζ | x] = 0, and since |⟨x, θZ ⟩| ≤ 1, the noise satisfies |ζ| ≤ 2 almost surely, so the instances lie e k. in the model class of Assumptions 1–3 with parameter in Θ In Appendix D.1 we bound the sum over the d coordinates of the total variation distance between the two conditional distributions of the output given Zi = +1 and Zi = −1, which gives the following result. Theorem 1. Let MSDP be the family of all non-interactive protocols that satisfy robust (ϵ, δ)ϵ,δ shuffled DP with 0 < ϵ ≤ O(1) and δ = o(1/n). Let Z = {(xi , yi )}ni=1 consist of n i.i.d. samples from the sparse linear model of Assumption 1 with parameter θ∗ . Then we have √   dk ∗ inf sup E ∥A(Z) − θ ∥2 ≥ Ω( ), (7) SDP nϵ A∈Mϵ,δ θ ∗ ∈Θ ek where the expectation is over the data and the randomness of the protocol. p
Central DP admits the rate Ω k log d/n + k log d/(nϵ) , while non-interactive local DP obeys p
Ω dk log d/(nϵ2 ) (Cai et al., 2021; Zhu et al., 2024). Theorem 1 places non-interactive SDP 5

Preprint

√ √ √ between the two. Written as dk/(nϵ) against dk log d/( nϵ), the shuffled bound is smaller √ √ than the local one by a factor n log d, which is the gain from shuffling. The factor d, however, is unchanged. This is the difference from the central model, where the privacy noise is added once to an aggregate statistic and the dimension enters only through the log d cost of locating the support; in the local and shuffled models each user randomizes their own record, and shuffling rescales the resulting per-coordinate noise without removing it. Hence the error of any non-interactive SDP protocol carries a polynomial dependence on d, whatever the local randomizer or analyzer, which is the primary bottleneck when d ≫ n. We next give an algorithm that achieves a nontrivial rate. 4.2

U PPER B OUND

The sparse estimation is actually an ℓ0 -constrained least squares problem, and this is the fundamental challenge when proposing an efficient non-interactive SDP algorithm. To control local noise scales, we therefore relax the ℓ0 constraint to an ℓ1 -regularized problem, whose solution can be written in closed form under Assumption 3. While this relaxation will link the error to the dimension d, the shuffler narrows this gap toward central DP. Then we can obtain the solution from the aggregated statistics by solving the normal equation and applying soft-thresholding. Our non-interactive SDP algorithm is given as Algorithm 2 in Appendix B. We assume the length D = d(d + 1)/2 + d = Θ(d2 ). e k and Assumptions 1–3 hold, with 0 < ϵ ≤ O(1) and the parameters of Theorem 2. Let θ∗ ∈ Θ Algorithm 2 set as in Appendix B. If ( )! p d2 log n log(d/δ) log d d2 log d D log(D/δ) n ≥ Ω max , , , (8) ϵκ∞ κ2∞ ϵ2 then Algorithm 2 satisfies (ϵ, δ)-shuffled DP, and with probability at least 1 − O(d−2 ) its output satisfies r √ ! k d k ∗ e ∥θ̂ − θ ∥2 ≤ O + . (9) n nϵ Privacy follows in two steps. The protocol that receives the rounded statistics directly and perturbs their sum with the binomial noise Bin(M, 1/2) − M/2 is (ϵ, δ/2)-DP. By the security of the summation protocol, the aggregator’s view of the real protocol is within total variation distance e−ϵ δ/8 of this one, and Lemma F.2 of Balle et al. (2020) turns this distance into the stated (ϵ, δ)-shuffled DP guarantee. For the utility, the rounding error is of the same order as the binomial noise; we bound bound. The resulting privacy term √and then apply the soft-thresholding √ √ both in the ℓ∞ -norm d k/(nϵ) leaves a gap of d to the lower bound Ω( dk/(nϵ)) of Theorem 1. The gap comes from 2 the clipping radius: the second-moment statistic has ∥xi x⊤ i ∥F = ∥xi ∥2 = Θ(d), and to privatize it, this scale enters the noise directly. We leave closing this gap to future work; in either case, the polynomial dependence on d is not removed, which is the main barrier of the non-interactive setting.

5

OVERCOMING THE T HEORETICAL BARRIER

Section 4 shows that a polynomial dependence on d is unavoidable in the non-interactive setting. In this section we show that the barrier can be overcome when each user holds multiple samples and two rounds of interaction are allowed. Specifically, we give a two-round user-level SDP algorithm whose estimation error is independent of the dimension d once the local sample size m is large enough. In this section, we consider the following additional assumption. Assumption 4 (Restricted Eigenvalue). There is a constant κ > 0 such that v ⊤ Σv ≥ κ∥v∥22 for every S ⊆ [d] with |S| ≤ k, and every v with ∥vS c ∥1 ≤ 3∥vS ∥1 , and there is a constant C such that maxℓ∈[d] Σℓℓ ≤ 1 and λmax (ΣS,S ) ≤ C for every |S| ≤ k. This is the restricted eigenvalue condition of Bickel et al. (2009), the standard condition under which the Lasso attains its fast rate, and the cone constant 3 is the one produced by the basic inequality when the regularization parameter dominates the score (van de Geer & Bühlmann, 2009; Negahban et al., 2012). It holds with high probability for sub-Gaussian designs whose population covariance 6

Preprint

has a bounded minimal eigenvalue (Raskutti et al., 2010), and conditions of this form are standard in the DP sparse regression literature (Cai et al., 2021; Ma et al., 2024; Tian et al., 2025). Taking v supported on S gives λmin (ΣS ) ≥ κ. 5.1

U PPER B OUND

Instead of relaxing the ℓ0 constraint as in Section 4.2, we reduce the problem to a low-dimensional regression by first estimating the support of θ∗ . In the user-level regime, each user turns their m local samples into one local estimate of an active coordinate, and a heavy hitter step aggregates these estimates into the global support set (Bun et al., 2019; Ma et al., 2024). Algorithm 1 User-level Two-Round Shuffled DP Sparse Regression (U2R) 1: Input: Local datasets {(Xi , yi )}n i=1 , sparsity k, privacy parameters (ϵ, δ). 2: Stage 1: Support Recovery (users i ∈ {1, . . . , n1 }) 3: for each user i ∈ {1, . . . , n1 } do 4: Solve the Lasso, get a local candidate set Ŝi , and draw one index vi ∼ Unif(Ŝi ). 5: Send the randomized response of vi on each coordinate ℓ ∈ [d] through d parallel shufflers. 6: end for 7: Aggregator: decode the shuffled messages to get the counts {Ŵℓ }ℓ∈[d] , recover support Ŝ via heavy hitter identification, and release Ŝ, where s = |Ŝ|. 8: Stage 2: Constrained Estimation (the remaining n2 = n − n1 users) 9: for each user i ∈ {n1 + 1, . . . , n} do 10: Compute the local least squares estimate θ̃i ∈ Rs on Ŝ, clip it, and send it to the shuffler. 11: end for 12: Aggregator: decode the messages of the first half users into a coarse estimate and broadcast it; the second half messages clip around this estimate, and yield the refined estimate θ̂Ŝ ∈ Rs . 13: Output: the k-sparse estimate θ̂ ∈ Rd that equals θ̂Ŝ on Ŝ and zero elsewhere. We summarize the pipeline in Algorithm 1. In this algorithm, Stage 1: Support Recovery, Stage 2: Constrained Estimation, and their internal hyperparameters (e.g., clipping thresholds) can be found in Algorithm 3 and Algorithm 4 in the Appendix B. Algorithm 1 is two-round because we count rounds by the messages that users send. The heavy hitter step of Stage 1 (Algorithm 3) is what keeps the dimension out of the final rate, so we describe it in more detail. The heavy hitter algorithms used for the LDP model (Bun et al., 2019; Ma et al., 2024) estimate the frequency of every coordinate through a prefix tree: the users are split into ⌈log d⌉ groups, with the users of one group reporting a hashed bit of one level of the tree, and the server scans the levels in order. Two costs come with this. Every report has to be private on its own, so the randomization a user applies is calibrated to their own message, and each level of the tree is estimated from a 1/ log d fraction of the users. In the shuffled model, neither cost is needed: the shuffler removes the association between a message and its sender, so the aggregator sees nothing beyond the counts, and each user can answer at a randomization level calibrated to the view VP rather than to a single message. We therefore replace the prefix tree by one round of randomized response over [d]. User i reduces their m local samples to a single index vi , sends the indicator of vi with every coordinate flipped independently with probability γ, and the aggregator debiases the shuffled counts and keeps the coordinates that pass a threshold. To ensure Ŝ is the true support, we give Proposition 1. Proposition 1. Let Assumptions 1, 2 and 4 hold, let S ∗ = supp(θ∗ ) with |S ∗ | = k, and let the parameters of Algorithm 3 be set as in Appendix B. Suppose there is ρ > 0 with minℓ∈S ∗ |θℓ∗ | ≥ ρ, and that (i) Local sample size.  m ≥ Ω max



 σ 2 k log d k log(ed/k) , , k + log n . κ2 ρ2 κ2 7

(10)

Preprint

(ii) User number of Stage 1. ( p )! k log(dn) log(1/δ) n1 ≥ Ω max , k log(dn) . ϵ

(11)

Then with probability at least 1 − O(n−2 ), Ŝ = S ∗ ,

(12)

so that s = |Ŝ| = k. Proposition 1 gives more than the containment S ∗ ⊆ Ŝ that Stage 2 needs: although the local set Ŝi of a user may differ from S ∗ , the heavy hitter step returns S ∗ exactly. In equation 10, the first two terms make the local selection of one user in Stage 1 correct with probability at least 1 − 1/(32k), independently of n, and the third term is what Stage 2 needs for the local estimate on Ŝ. p least squares
The first term is equivalently the beta-min condition ρ ≥ Ω (σ/κ) k log d/m , which is standard for support recovery (Wainwright, 2009; Ma et al., 2024) and allows ρ to decrease as m grows. Theorem 3. Let θ∗ ∈ Θk , let the conditions of Proposition 1 hold, and let the parameters of Algorithm 1 be set as in Appendix B. If ! p m log(k/δ) nϵ ≥ Ω , (13) σ then Algorithm 1 satisfies user-level (ϵ, δ)-shuffled DP, and with probability at least 1 − O(n−2 ) its output satisfies ! r k k e ∥θ̂ − θ∗ ∥2 ≤ O + √ . (14) nm n mϵ The dimension d is absent from both terms. By Proposition 1 the second stage runs on the k coordinates of Ŝ, so the vector a user sends has k entries rather than d, and the dimension is eliminated. What remains of d is in the conditions on m and n1 in Proposition 1. In the synthetic dimension sweep, U2R keeps a low mean relative coefficient error as d increases at fixed k and m, as shown in Table 6. User-level data alone does not suffice: Appendix F.1 shows that under a non-interactive protocol the lower bound for user-level SDP still carries a polynomial dependence on d, so interaction is also essential. 5.2

L OWER B OUND

We have already shown a promising result in Theorem 3. However, a natural question is whether our upper bound can be further improved. Now we show that our upper bound is already near-optimal. Theorem 4. Let Z = {Ui }ni=1 be a dataset in which each user i ∈ [n] holds m i.i.d. samples Ui = {(xi,j , yi,j )}m j=1 from the exponential family   ⊤ yx θ − ψ(x⊤ θ) . (15) fθ (y | x) = h(y, σ) exp c(σ) Let Assumptions 2 and 4 hold and let ∥ψ ′′ ∥∞ ≤ c2 . Let Mϵ,δ be p the family of all estimators √
satisfying user-level (ϵ, δ)-DP with 0 < ϵ < 1 and δ = o log(d/k)/(n m log(nm)) . If k ≤ d,
√ k ≥ Ω(log2 (d/k)) and n m ϵ ≥ Ω k 3/4 log3/2 (d/k) , then ! r   k log(d/k) k log(d/k) ∗ √ inf sup E ∥A(Z) − θ ∥2 ≥ Ω + , (16) A∈Mϵ,δ θ ∗ ∈Θk nm n mϵ where the expectation is over the data and the randomness of the estimator. Theorem 4 defines the fundamental efficiency limits for any differential private estimator in the user-level setting. It is stated for the exponential family, which contains the sparse linear model of 8

Preprint

Assumption 1 as the special case ψ(u) = u2 /2 and c(σ) = σ 2 , so the bound applies to the model class we estimate. By Definition 3, the output of any (ϵ, δ)-shuffled DP protocol is (ϵ, δ)-DP in the sense of Definition 1, so Theorem 4 also lower bounds the error of any user-level SDP protocol. Comparing it with Theorem 3, the two rates agree in n, m, k and ϵ and only differ by the factor log(d/k): our algorithm is near-optimal, and sparse linear regression in the shuffled model is closer to central DP once interaction and enough user-level data are available. User-level Score Attack. We now outline the proof. Minimax lower bounds under privacy constraints are difficult to characterize: Fano’s inequality and mutual information bounds give no general recipe and are often loose, particularly in the user-level regime. To handle user-level data, we extend the score attack framework of Cai et al. (2023) and introduce the user-level score attack, which is our main technical tool here. Definition 7 (User-level Score Attack). For a dataset Z in which each user i ∈ [n] holds m i.i.d. samples Ui = {zi,1 , . . . , zi,m }, the user-level score attack for an estimator A(Z) is * + m X Ai := A(Z) − θ, Sθ (zi,j ) , (17) j=1

where Sθ (z) = ∇θ log fθ (z) is the score statistic. For the exponential family of Theorem 4, with mean function µ(x⊤ θ) = ψ ′ (x⊤ θ), the score statistic is Sθ (y, x) = (y − µ(x⊤ θ))x/c(σ), so the user-level score attack takes the form + * m X
1 ⊤ Ai = yi,j − µ(xi,j θ) xi,j . (18) A(Z) − θ, c(σ) j=1 Proposition 2. Let Z = {U1 , . . . , Un } be as in Theorem 4, let A(Z) be an estimator of θ, and let ′ ei . Write Z(i) be the adjacent dataset in which the i-th user’s data Ui is replaced by an i.i.d. copy U ′ ′ Ai = ⟨A(Z) − θ, Sθ (Ui )⟩ and Ai = ⟨A(Z(i) ) − θ, Sθ (Ui )⟩. Then (i) Soundness: for each i ∈ [n], E[A′i ] = 0 and r q
m E|A′i | ≤ E∥A(Z) − θ∥22 · λmax E[ψ ′′ (x⊤ θ)xx⊤ ] ; c(σ)

(19)

(ii) Completeness: under standard regularity conditions, n X

E[Ai ] =

i=1

d X ∂ E[A(Z)j ]. ∂θ j j=1

(20)

Pn With Proposition 2, we bound the expected sum i=1 E[Ai ] from above by Lemma 3 and from below by Lemma 4. Comparing the two gives Theorem 4. Numerical results. On synthetic data, U2R attains mean relative coefficient error 0.00631 at the default setting, comparable to non-private Lasso. On real-world dataset Kelmarsh, it has the lowest mean test MSE among the private methods at both ϵ = 1 and 4 for each m ∈ {10, 50, 100}; Appendix C reports the detailed setup and complete results.

6

C ONCLUSION

We studied sparse linear regression under shuffled DP. In the non-interactive setting, we showed a lower bound on the ℓ2 -norm estimation error of any (ϵ, δ)-SDP protocol and gave an algorithm with a matching form, both bounds carrying a polynomial dependence on d. In the interactive user-level setting, we gave a two-round algorithm that first applies a heavy hitter step and then estimates the parameter, whose error does not depend polynomially on d. By extending the score attack to userlevel data, we proved a lower bound for any user-level (ϵ, δ)-DP estimator, which shows that this algorithm is near-optimal. Experiments on synthetic and real-world datasets support these results. 9

Preprint

AI USE STATEMENT In this work, we used generative AI tools as coding assistants to implement the experimental algorithms. We have not used generative AI tools for drafting sections of the paper, generating synthetic datasets, or proving mathematical claims. Additionally, we used generative AI tools to revising wording and for organizing the bibliography. We have reviewed all AI-assisted work: the authors independently examined and verified all experimental scripts and mathematical derivations. We take full responsibility for the final content of this work, including text or artifacts produced with the aid of generative AI. E THICS STATEMENT Our experiments use the publicly available Kelmarsh wind-farm dataset, released by Cubico Sustainable Investments Ltd under the Creative Commons Attribution 4.0 International (CC BY 4.0) license. We use the dataset in accordance with its license terms and provide appropriate attribution. The experiments analyze wind-turbine operational measurements and involve neither human participants nor personally identifiable information. R EPRODUCIBILITY STATEMENT We have made every effort to ensure that our results can be reproduced from the paper and its appendix. For the theoretical results, the problem setting and assumptions are stated explicitly in Section 3.2 (Assumptions 1 and 2). The additional conditions used in each regime are given in Assumption 3 (Section 4) and Assumption 4 (Section 5). Complete proofs of all claims are provided in the appendix. Theorems 1 and 2 are proved in Appendix D. Proposition 1, Theorem 3, Proposition 2 and Theorem 4 are proved in Appendix E. The supplementary results on non-interactive user-level SDP and dense regression, together with their proofs, are in Appendix F. For the algorithms, Appendix B gives full pseudocode for the non-interactive protocol (Algorithm 2) and for both stages of U2R (Algorithms 3 and 4). Tables 3–5 list every parameter required by the theorems, including clipping radii, rounding precisions, binomial noise parameters, moduli, numbers of shares, and thresholds. For the experiments, Appendix C.1 describes the synthetic data-generating process, the default values and sweep ranges, the number of repetitions, and the evaluation metrics. It also gives the hyperparameters of U2R and of all baselines, and how they were selected. The real-data experiments use the publicly available Kelmarsh wind farm dataset (CC BY 4.0). The same appendix details the full preprocessing pipeline, including feature construction, the date-based public/training/test split, imputation, standardization, and the assignment of rows to simulated users. We do not release source code with this submission. We believe these descriptions are sufficient to reimplement the algorithms and reproduce the reported results.

R EFERENCES Jayadev Acharya, Clément L. Canonne, Himanshu Tyagi, and Ziteng Sun. The role of interactivity in structured estimation. In Proceedings of the 35th Conference on Learning Theory, pp. 1328– 1355. PMLR, 2022. Hilal Asi, Vitaly Feldman, Tomer Koren, and Kunal Talwar. Private stochastic convex optimization: Optimal rates in ℓ1 geometry. In Proceedings of the 38th International Conference on Machine Learning, pp. 393–403. PMLR, 2021. Victor Balcer and Albert Cheu. Separating local & shuffled differential privacy via histograms. arXiv preprint arXiv:1911.06879, 2019. Borja Balle, James Bell, Adrià Gascón, and Kobbi Nissim. The privacy blanket of the shuffle model. In Advances in Cryptology – CRYPTO 2019, pp. 638–667. Springer, 2019. Borja Balle, James Bell, Adrià Gascón, and Kobbi Nissim. Private summation in the multi-message shuffle model. In Proceedings of the 2020 ACM SIGSAC Conference on Computer and Communications Security, pp. 659–676, 2020. 10

Preprint

Raef Bassily, Vitaly Feldman, Kunal Talwar, and Abhradeep Guha Thakurta. Private stochastic convex optimization with optimal rates. In Advances in Neural Information Processing Systems 32. Curran Associates, Inc., 2019. Peter J. Bickel, Ya’acov Ritov, and Alexandre B. Tsybakov. Simultaneous analysis of lasso and dantzig selector. The Annals of Statistics, 37(4):1705–1732, 2009. Mark Bun, Jelani Nelson, and Uri Stemmer. Heavy hitters and the structure of local privacy. ACM Transactions on Algorithms, 15(4):1–40, 2019. T. Tony Cai, Yichen Wang, and Linjun Zhang. The cost of privacy: Optimal rates of convergence for parameter estimation with differential privacy. The Annals of Statistics, 49(5):2825–2850, 2021. T. Tony Cai, Yichen Wang, and Linjun Zhang. Score attack: A lower bound technique for optimal differentially private learning. arXiv preprint arXiv:2307.01801, 2023. Emmanuel Candès and Terence Tao. The dantzig selector: Statistical estimation when p is much larger than n. The Annals of Statistics, 35(6):2313–2351, 2007. Kamalika Chaudhuri, Claire Monteleoni, and Anand D. Sarwate. Differentially private empirical risk minimization. Journal of Machine Learning Research, 12:1069–1109, 2011. Albert Cheu and Jonathan Ullman. The limits of pan privacy and shuffle privacy for learning and estimation. In Proceedings of the 53rd Annual ACM Symposium on Theory of Computing (STOC), 2021. Albert Cheu, Adam Smith, Jonathan Ullman, Maxim Zeber, and Sergey Zhilyaev. Distributed differential privacy via shuffling. In Advances in Cryptology – EUROCRYPT 2019, volume 11476 of Lecture Notes in Computer Science, pp. 375–403. Springer, 2019. Amit Daniely and Vitaly Feldman. Locally private learning without interaction requires separation. In Advances in Neural Information Processing Systems 32. Curran Associates, Inc., 2019. John Duchi, Michael Jordan, and Martin Wainwright. Local privacy and statistical minimax rates. In Proceedings of the 54th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pp. 429–438, Berkeley, CA, USA, 2013. IEEE. doi: 10.1109/FOCS.2013.53. URL https: //arxiv.org/abs/1302.3203. John C. Duchi, Michael I. Jordan, and Martin J. Wainwright. Privacy Aware Learning. Journal of the ACM, 61(6):1–57, 2014. Cynthia Dwork and Aaron Roth. The Algorithmic Foundations of Differential Privacy. Now Publishers Inc., 2014. Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam Smith. Calibrating noise to sensitivity in private data analysis. In Theory of Cryptography Conference, pp. 265–284. Springer, 2006. Úlfar Erlingsson, Vitaly Feldman, Ilya Mironov, Ananth Raghunathan, Kunal Talwar, and Abhradeep Thakurta. Amplification by shuffling: From local to central differential privacy via anonymity. In Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 2468–2479. SIAM, 2019. Vitaly Feldman, Tomer Koren, and Kunal Talwar. Private stochastic convex optimization: Optimal rates in linear time. In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing, pp. 439–449, 2020. Badih Ghazi, Pasin Manurangsi, Rasmus Pagh, and Ameya Velingker. Private aggregation from fewer anonymous messages. In Advances in Cryptology – EUROCRYPT 2020, pp. 798–827. Springer, 2020. Antonious M. Girgis and Suhas Diggavi. Multi-message shuffled privacy in federated learning. IEEE Journal on Selected Areas in Information Theory, 5:12–27, 2024. 11

Preprint

Lijie Hu, Shuo Ni, Hanshen Xiao, and Di Wang. High dimensional differentially private stochastic optimization with heavy-tailed data. In Proceedings of the 41st ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems, pp. 227–236, 2022. Prateek Jain, Ambuj Tewari, and Purushottam Kar. On iterative hard thresholding methods for highdimensional m-estimation. In Advances in Neural Information Processing Systems 27. Curran Associates, Inc., 2014. Gautam Kamath, Vikrant Singhal, and Jonathan Ullman. Private mean estimation of heavy-tailed distributions. In Proceedings of the 33rd Conference on Learning Theory, pp. 2204–2235. PMLR, 2020. Gautam Kamath, Xingtu Liu, and Huanyu Zhang. Improved rates for differentially private stochastic convex optimization with heavy-tailed data. In Proceedings of the 39th International Conference on Machine Learning, pp. 10633–10660. PMLR, 2022. Shiva Prasad Kasiviswanathan, Homin K. Lee, Kobbi Nissim, Sofya Raskhodnikova, and Adam Smith. What can we learn privately? SIAM Journal on Computing, 40(3):793–826, 2011. Daniel Kifer, Adam Smith, and Abhradeep Thakurta. Private Convex Empirical Risk Minimization and High-dimensional Regression. In Conference on Learning Theory, pp. 25.1–25.40. JMLR Workshop and Conference Proceedings, 2012. Xiyang Liu, Prateek Jain, Weihao Kong, Sewoong Oh, and Arun Sai Suggala. Near optimal private and robust linear regression. arXiv preprint arXiv:2301.13273, 2023. Yuheng Ma, Ke Jia, and Hanfang Yang. Better locally private sparse estimation given multiple samples per user. In Proceedings of the 41st International Conference on Machine Learning (ICML). PMLR, 2024. Sahand N. Negahban, Pradeep Ravikumar, Martin J. Wainwright, and Bin Yu. A unified framework for high-dimensional analysis of m-estimators with decomposable regularizers. Statistical Science, 27(4):538–557, 2012. Garvesh Raskutti, Martin J. Wainwright, and Bin Yu. Restricted eigenvalue properties for correlated gaussian designs. Journal of Machine Learning Research, 11:2241–2259, 2010. Garvesh Raskutti, Martin J. Wainwright, and Bin Yu. Minimax rates of estimation for highdimensional linear regression over ℓq balls. IEEE Transactions on Information Theory, 57(10): 6976–6994, 2011. Mark Rudelson and Shuheng Zhou. Reconstruction from anisotropic random measurements. IEEE Transactions on Information Theory, 59(6):3434–3447, 2013. Adam Smith, Abhradeep Thakurta, and Jalaj Upadhyay. Is Interaction Necessary for Distributed Private Learning? In 2017 IEEE Symposium on Security and Privacy (SP), pp. 58–77. IEEE, 2017. Qiang Sun, Wen-Xin Zhou, and Jianqing Fan. Adaptive huber regression. Journal of the American Statistical Association, 115(529):254–265, 2020. Kunal Talwar, Abhradeep Guha Thakurta, and Li Zhang. Nearly Optimal Private LASSO. In Advances in Neural Information Processing Systems, volume 28. Curran Associates, Inc., 2015. Xizhi Tian, Meng Ding, Touming Tao, Zihang Xiang, and Di Wang. Differentially private sparse linear regression with heavy-tailed responses. In Machine Learning and Knowledge Discovery in Databases: Research Track (ECML PKDD), pp. 363–379. Springer, 2025. Sara A. van de Geer and Peter Bühlmann. On the conditions used to prove oracle results for the lasso. Electronic Journal of Statistics, 3:1360–1392, 2009. Prateek Varshney, Abhradeep Thakurta, and Prateek Jain. (nearly) optimal private linear regression via adaptive clipping. arXiv preprint arXiv:2207.04686, 2022. 12

Preprint

Martin J. Wainwright. Sharp thresholds for high-dimensional and noisy sparsity recovery using ℓ1 constrained quadratic programming (Lasso). IEEE Transactions on Information Theory, 55(5): 2183–2202, 2009. doi: 10.1109/TIT.2009.2016018. Di Wang and Jinhui Xu. On sparse linear regression in the local differential privacy model. IEEE Transactions on Information Theory, 67(2):1182–1200, 2020. Di Wang, Hanshen Xiao, Srinivas Devadas, and Jinhui Xu. On differentially private stochastic convex optimization with heavy-tailed data. In Proceedings of the 37th International Conference on Machine Learning, pp. 10081–10091. PMLR, 2020. Di Wang, Huangyu Zhang, Marco Gaboardi, and Jinhui Xu. Estimating smooth glm in noninteractive local differential privacy model with public unlabeled data. In Proceedings of the 32nd International Conference on Algorithmic Learning Theory, pp. 1207–1213. PMLR, 2021. Eunho Yang, Aurelie Lozano, and Pradeep Ravikumar. Elementary estimators for high-dimensional linear regression. In Proceedings of the 31st International Conference on Machine Learning, pp. 388–396. PMLR, 2014. Kai Zheng, Wenlong Mou, and Liwei Wang. Collect at Once, Use Effectively: Making Noninteractive Locally Private Learning Possible. In Proceedings of the 34th International Conference on Machine Learning, pp. 4130–4139. PMLR, 2017. Liyang Zhu, Meng Ding, Vaneet Aggarwal, Jinhui Xu, and Di Wang. Improved analysis of sparse linear regression in local differential privacy model. In The Twelfth International Conference on Learning Representations (ICLR), 2024.

A

D ETAILED C OMPARISON AND R ELATED W ORK

Table 2 provides a comprehensive comparison of our theoretical results with prior literature on (ϵ, δ)DP linear regression. This extends the summary provided in the main text by including bounds for specific edge cases (such as 1-sparse models) and detailing previous interactive local methods. Private high-dimensional regression has been studied extensively in the central model of differential privacy. Early work by Chaudhuri et al. (2011) and Kifer et al. (2012) showed that private empirical risk minimization remains feasible even in high dimension, and Talwar et al. (2015) designed a nearly optimal private LASSO algorithm under an ℓ1 constraint. The central-DP optimization literature has since developed a much sharper understanding of stochastic convex optimization and private linear regression, including optimal-rate results in different geometries and heavy-tailed regimes (Bassily et al., 2019; Feldman et al., 2020; Asi et al., 2021; Wang et al., 2020; Kamath et al., 2022; Hu et al., 2022; Liu et al., 2023; Varshney et al., 2022). On the statistical side, Cai et al. (2021) and Cai et al. (2023) established minimax rates and lower-bound techniques for private parameter estimation, showing that sample-level central DP only incurs logarithmic dependence on the ambient dimension in sparse problems. The score attack of Cai et al. (2023) is stated for datasets in which each record is one sample; the user-level score attack of Section 5.2 extends it to the regime where a record is the whole local dataset of a user, which is what Theorem 4 needs. These results provide the benchmark that one would ideally like to approach in decentralized settings. The local model is considerably more restrictive. Foundational results on local privacy and statistical minimax rates are due to Duchi et al. (2013) and Duchi et al. (2014), while later work studied the role of interaction in distributed private learning and local estimation (Smith et al., 2017; Zheng et al., 2017; Daniely & Feldman, 2019; Hu et al., 2022; Acharya et al., 2022; Wang et al., 2021). For sparse linear regression in particular, Wang & Xu (2020) initiated the study of the local model, proving lower bounds in the 1-sparse case and proposing an interactive private iterative hard-thresholding procedure under strong distributional assumptions. The later improved analysis of Zhu et al. (2024) substantially sharpened this picture in two directions: it extended the lower bounds to general ksparse regression, and it provided both the first explicit non-interactive LDP estimator for sparse regression and a corrected utility analysis for the interactive hard-thresholding line. Their results show that, even after these improvements, sample-level LDP sparse regression still exhibits a polynomial dependence on the ambient dimension d, with non-interactive and interactive upper bounds 13

Preprint

Table 2: Detailed comparison of related studies on (ϵ, δ)-DP linear regression. Data assumptions and parameter definitions are identical to Table 1. Type

Unit

Setting Upper Bound

Lower Bound

sample

general

sample

q d Õ( nd + nϵ ) q k k-sparse Õ( nk + nϵ )

q d Ω( nd + nϵ ) q d d Ω( k log + k log ) n nϵ

sample

1-sparse -

q log d Ω( dnϵ 2 )

sample

k-sparse Õ( d√nϵk )

Source

Central DP -

Cai et al. (2021) Cai et al. (2021)

Local DP Non-interactive

√

Interactive

Ω(

3 d2

Õ( √nϵ )

sample

general

sample

1-sparse -

q

dk log d ) nϵ2

-

Wang & Xu (2020) Zhu et al. (2024) Smith et al. (2017)

q Ω( nϵd2 )

q √ dk k-sparse Õ( k√nϵd ) Ω( nϵ 2) q k k k user-level k-sparse Õ( nm + √nmϵ ) Ω( √nmϵ ) sample

Wang & Xu (2020) Zhu et al. (2024) Ma et al. (2024)

Shuffle DP Non-interactive

user-level general

user-level k-sparse -

q √ √ dk k-sparse Õ( nk + dnϵk ) Ω( nϵ ) q q k log(d/k) k k √ user-level k-sparse Õ( nm + n√mϵ ) Ω( + k log(d/k) ) nm n mϵ

sample Interactive

q

d + n√dmϵ ) q nm √ Ω( k log(d/k) + n√dk ) nm mϵ

Ω(

-

14

Our Work

Preprint

√ √ √ √ e e scaling on the order of O(d k/( nϵ)) and O(k d/( nϵ)), respectively. These phenomena are closely related to the non-private high-dimensional estimation literature, where canonical sparsityaware procedures include the Dantzig selector (Candès & Tao, 2007), elementary estimators (Yang et al., 2014), and iterative hard thresholding (Jain et al., 2014); minimax rates over sparse parameter classes are summarized by Raskutti et al. (2011). Several papers have therefore looked for stronger primitives or better structural assumptions. On the robust statistics side, clipped or Huberized estimators are now standard tools for controlling heavytailed noise and stabilizing private regression (Sun et al., 2020; Kamath et al., 2020); these ideas also appear explicitly in the improved non-interactive LDP construction of Zhu et al. (2024), where clipped sufficient statistics are privatized locally before sparse recovery is performed at the server. Both of our protocols clip in the same way, and Algorithm 2 privatizes the same clipped sufficient statistics, the difference being that the noise is added to the aggregate rather than to each record. On the decentralized side, Ma et al. (2024) moved to the user-level multi-sample regime and proposed the two-stage 2-SLR framework, where each user first locally screens candidate coordinates before a second-stage low-dimensional estimator is applied. Their result is especially relevant to our work because it shows that once multiple samples are grouped per user, the severe dimension dependence of standard sample-level LDP can be substantially alleviated. Algorithm 1 follows the same twostage structure, and the local screening step of Algorithm 3 is a Lasso together with a threshold, which asks only for the restricted eigenvalue condition of Assumption 4 rather than for uncorrelated features. Shuffled differential privacy offers a middle ground between central DP and LDP. The basic model and its privacy amplification effect were formalized by Cheu et al. (2019) and Erlingsson et al. (2019), and the amplification view was later refined through the privacy-blanket perspective (Balle et al., 2019). A separate line studies multi-message protocols that compute a sum in the shuffled model with the accuracy of the central model, both for real summation and for aggregation in federated learning (Balle et al., 2020; Ghazi et al., 2020; Girgis & Diggavi, 2024). The summation step of Algorithm 2 and of Stage 2 of Algorithm 1 is of this second type: the aggregator sees a noisy sum rather than n separately randomized messages, so the noise can be calibrated to the ℓ2 sensitivity of the query. On the negative side, separations and limitations of shuffle privacy have been established for statistical tasks such as histogram estimation, together with the relation between shuffled DP and pan privacy (Balcer & Cheu, 2019; Cheu & Ullman, 2021), which is the route our non-interactive lower bound follows. Heavy-hitter primitives under local noise provide a practical mechanism for support identification in decentralized settings (Bun et al., 2019); Stage 1 of Algorithm 1 replaces the prefix tree of these protocols by one round of randomized response over [d]. However, despite this rapid progress, sparse linear regression in the shuffled model has remained largely open: the existing literature either analyzes general shuffled estimation primitives or studies local/user-level baselines rather than giving a sharp treatment of sparse regression itself. Our work fills this gap by characterizing both the non-interactive sample-level shuffled regime and the interactive user-level shuffled regime for sparse linear models.

B

D ETAILED A LGORITHMS AND PARAMETER S ETTINGS

This section gives the pseudocode of the algorithms of Sections 4 and 5, together with the parameter settings used by Theorem 2, Proposition 1 and Theorem 3. Throughout, D = d(d + 1)/2 + d = Θ(d2 ), v(·) vectorizes the upper triangle of a symmetric matrix including the diagonal, ΠB(0,R) is the Euclidean projection onto the ball of radius R, and c0 , c1 , cA , c′A , cλ , cp , Bθ , Cβ , Cγ are absolute constants. Both protocols privatize a sum of integer vectors in the same way: each user rounds its statistic to the integer grid, adds a share of a centered binomial noise, and uploads additive shares of the result through parallel shufflers, one per coordinate. The noise is calibrated to the ℓ2 sensitivity of the integer query rather than to its ℓ1 sensitivity, which is what keeps the privacy term at d in Theorem 2 and removes d altogether in Theorem 3. Algorithm 2 runs on all n users at the budget (ϵ, δ). In Algorithm 1 the three groups of users, namely the n1 users of Stage 1, the n3 users of Stage 2a and the n4 users of Stage 2b, are disjoint, so a neighboring pair of datasets changes the view of one group only and each group runs at the 15

Preprint

full budget (ϵ, δ) by parallel composition. The parameter settings below are stated accordingly. The summation step requires at least 19 participants per group. B.1

N ON - INTERACTIVE A LGORITHM

Algorithm 2 Sample-Level Non-Interactive Shuffled DP Sparse Linear Regression 1: Input: Data {(xi , yi )}n i=1 , sparsity k, privacy parameters (ϵ, δ), parameters of Table 3. 2: Local procedure (user i ∈ [n]) D 3: Compute Hi = xi x⊤ i , gi = xi yi , and set Fi = (v(Hi ); gi ) ∈ R . 4: Clip: F̄i = Fi · min{1, CF /∥Fi ∥2 }. 5: For each j ∈ [D], round: Ūi,j = ⌊pF̄i,j ⌋ + Ber(pF̄i,j − ⌊pF̄i,j ⌋). 6: For each j ∈ [D], draw ζi,j ∼ Bin(N0 , 1/2) − N0 /2 and set Wi,j = (Ūi,j + ζi,j ) mod q. 7: For each j ∈ [D], split Wi,j into r additive shares in Zq and send them, together with one unshuffled share, through the j-th shuffler. 8: Shuffler: permute the messages of each of the D streams independently. 9: Aggregator 10: For each j ∈ [D], sum the received shares modulo q to get Sj , and set Sj ← Sj − q if Sj > q/2. 11: Rescale F̂j = Sj /(np), fill Ĥ ∈ Rd×d from the first d(d + 1)/2 entries of F̂ and ĝ ∈ Rd from

the last d.

12: Solve θ̃ = Ĥ −1 ĝ, soft-threshold θ̂ = Sλn (θ̃), and project θ̂ onto {θ : ∥θ∥1 ≤ 1}. 13: Output: the k-sparse estimate θ̂.

The parameter N0 is even, so the noise of a user is integer valued, and N0 ≥ 2, so every user holds at least one Bernoulli trial per coordinate and the shares of the n users add up to a centered Bin(nN0 , 1/2) noise on each coordinate. The choice of p keeps the rounding error of the same order as the binomial noise, so that both are absorbed into λn . The modulus q is large enough that the sum √ of the uploaded values does not wrap around. Under Assumption 3 we have ∥Fi ∥2 ≤ c2x d+cx cy d, so the clipping step does not change Fi once c0 ≥ c2x + cx cy ; the step is kept so that the algorithm also applies to data that are not bounded, where the factor log n in CF covers the tail of ∥Fi ∥2 , which then exceeds CF with probability at most n−3 . Table 3: Parameters of Algorithm 2, used in Theorem 2. Symbol

Setting

Role

CF

c0 d log n √ l m nϵ p d log n log(4D/δ) √ 2pCF + D l 32(∆(p) )2 log(4D/δ) m

ℓ2 clipping radius

p (p)

∆2 M N0 q r λn

B.2

rounding precision ℓ2 sensitivity of the integer query

2

total binomial parameter

ϵ2 2⌈M/(2n)⌉

binomial parameter per user and coordinate

4n(pCF + N0 ) l 2 log (8Deϵ /δ) + log q m

modulus of the group Zq

2 2 +1 log2 n − log2 e p r CF log(4D/δ) log D  c1  log d σ + κ∞ n nϵ

number of additive shares per coordinate soft-thresholding parameter

S TAGE 1: S UPPORT R ECOVERY

The window for λloc is non-empty under the first term of equation 10: the lower end is what dominates the score ∥m−1 Xi⊤ ζi ∥∞ of one user with probability at least 1 − 1/(64k), and the upper end is what turns the ℓ2 rate of the Lasso into the coordinatewise control needed by the threshold ρ/2. 16

Preprint

Table 4: Parameters of Algorithm 3, used in Proposition 1. Symbol

Setting r

λloc

cλ σ

ρ/2

— Cγ log(1/δ) n1 ϵ2 n1 4k

γ τHH

Role log d κρ , with λloc ≤ √ m 16 k

regularization parameter of the local Lasso threshold of the local selection step flipping probability of the randomized response threshold on the debiased counts

Algorithm 3 Heavy Hitter Support Recovery (Stage 1) 1 1: Input: Local datasets {(Xi , yi )}n i=1 , sparsity k, signal strength ρ, privacy parameters (ϵ, δ), parameters of Table 4. 2: Local procedure (user i ∈ [n1 ]) 1 3: Solve the Lasso θ̂ (i) = arg minθ∈Rd 2m ∥yi − Xi θ∥22 + λloc ∥θ∥1 . (i) 4: Set Ŝi = {ℓ ∈ [d] : |θ̂ℓ | ≥ ρ/2}, and Ŝi = [d] if this set is empty. 5: Draw vi ∼ Unif(Ŝi ). 6: For each ℓ ∈ [d], send riℓ = I[vi = ℓ] ⊕ Ber(γ) through the ℓ-th shuffler. 7: Shuffler: permute the messages of each of the d streams independently. 8: Aggregator Pn1 9: For each ℓ ∈ [d], form the count Cℓ = i=1 riℓ and debias it: Ŵℓ = (Cℓ − n1 γ)/(1 − 2γ). 10: Set Ŝ = {ℓ ∈ [d] : Ŵℓ ≥ τHH } and s = |Ŝ|. 11: Output: the support estimate Ŝ.

The value of γ is the smallest one at which the counts are private. It satisfies γ ≤ 1/4 as long as 2 2 δ ≥ (dn)−Ck ϵ , so that the debiasing factor 1/(1 − 2γ) is at most 2, and under equation 11 the fluctuation of the debiased counts stays below n1 /(16k), so that τHH separates the two groups of coordinates. B.3

S TAGE 2: C ONSTRAINED E STIMATION

The n2 = n − n1 remaining users are split into two groups of sizes n3 = ⌊n2 /2⌋ and n4 = n2 − n3 . Both groups upload through the same summation step as Algorithm 2, with s parallel shufflers instead of D. The quantity  broadcast with τ is a function of the public parameters and of Ŝ, which has already been released, so broadcasting τ costs no privacy budget. The two groups of Stage 2 use the same summation step with different sizes and clipping radii. To state their parameters once, we index them by g ∈ {a, b} and set (na , Ra ) = (n3 , R0 ),

(nb , Rb ) = (n4 , τ ),

(g)

so that pg , ∆g , Mg , N0 , qg , rg in Table 5 are the parameters of Stage 2a for g = a and of Stage 2b for g = b. The constant Bθ is any known upper bound on ∥θ∗ ∥2 , which Assumption 1 provides; taking Bθ = 1 is enough on the parameter space we consider, √ so that R0 = Θ(1). The second term in the setting of pg makes ∆g ≤ 3pg Rg and ∆g ≥ s max{1, ϵ} hold for every value of s, at (g) the cost of a larger N0 , that is of more messages per user; it does not change the noise seen by thepaggregator. Indeed the precision pg cancels against the sensitivity ∆g , and that noise is p

Θ R0 s log n log(4s/δ)/(n3 ϵ) for the coarse estimate and Θ τ s log n log(4s/δ)/(n4 ϵ) for the refined one. Since τ = Θ(τβ ) under the condition equation 13, the second scale is what produces the privacy term of Theorem 3, and it involves s rather than d. 17

Preprint

Algorithm 4 Clipped Sparse Estimation (Stage 2) 1: Input: Support Ŝ with s = |Ŝ|, local datasets of the n2 users, privacy parameters (ϵ, δ), param-

eters of Table 5. 2: Local procedure (every user of Stage 2) 3: Compute the local least squares estimate θ̃i = (X ⊤ Xi,Ŝ )−1 X ⊤ yi ∈ Rs . i,Ŝ

i,Ŝ

4: Stage 2a (the first n3 users) 5: Clip θ̄i = ΠB(0,R0 ) (θ̃i ). (a)

(a)

6: For each j ∈ [s], round pa θ̄i,j to the integer grid, add Bin(N0 , 1/2) − N0 /2, reduce modulo

qa , and send ra additive shares together with one unshuffled share through the j-th shuffler.

7: Aggregator: decode the shuffled messages, rescale by n3 pa to get the coarse estimate c ∈ Rs ,

and broadcast τ = τβ + Â. 8: Stage 2b (the remaining n4 users) 9: Clip around the coarse estimate: ∆i = ΠB(0,τ ) (θ̃i − c). (b)

(b)

10: For each j ∈ [s], round pb ∆i,j to the integer grid, add Bin(N0 , 1/2) − N0 /2, reduce modulo

qb , and send rb additive shares together with one unshuffled share through the j-th shuffler. 11: Aggregator: decode the shuffled messages, rescale by n4 pb , and add back c to get θ̂Ŝ ∈ Rs . 12: Output: the estimate θ̂ ∈ Rd that equals θ̂Ŝ on Ŝ and zero elsewhere.

Table 5: Parameters of Algorithm 4, used in Theorem 3. The index g ∈ {a, b} refers to Stage 2a and Stage 2b, with (na , Ra ) = (n3 , R0 ) and (nb , Rb ) = (n4 , τ ). Symbol

Setting r

τβ

Cβ σ

R0

Bθ + τ β p r s log n log(4s/δ) s log n ′ R0 cA σ + cA n3 m n3 ϵ τβ + Â n cp png log n ϵ (√s + 2) max{1, ϵ} o p max , Rg Rg log(1/δ) √ 2pg Rg + max{1, ϵ} s + 2 l 32∆2 log(4s/δ) m g ϵ2 2⌈Mg /(2ng )⌉

 τ pg ∆g Mg (g) N0

qg rg

C

Role s log n m

(g) 4ng (pg Rg + N0 ) l 2 log (8seϵ /δ) + log q m 2 2 g

log2 ng − log2 e

uniform bound on ∥θ̃i − θŜ∗ ∥2 clipping radius of Stage 2a error bound of the coarse estimate c clipping radius of Stage 2b rounding precision ℓ2 sensitivity of the integer query; ∆g ≤ 3pg Rg total binomial parameter binomial parameter per user and coordinate modulus of the group Zqg +1

number of additive shares per coordinate

E XPERIMENTAL S ETUP AND R ESULTS

We evaluate U2R on synthetic sparse linear regression and nacelle temperature prediction using the Kelmarsh wind farm dataset. C.1

E XPERIMENTAL S ETUP

We compare U2R in Algorithm 1 with 2-SLR (Ma et al., 2024), LDP-IHT (Zhu et al., 2024), SampleSDP in Algorithm 2, and non-private Lasso fitted to all nm observations. U2R and 2-SLR protect users; LDP-IHT and Sample-SDP protect individual observations. U2R and Sample-SDP use shuffled (ϵ, δ)-DP with δ = 10−5 ; 2-SLR and LDP-IHT use pure ϵ-LDP. For the same protected unit, pure ϵ-DP implies (ϵ, 10−5 )-DP, so we compare utility at the same numerical budget while retaining the differences in privacy units and trust models. 18

Preprint

For synthetic data, each user holds m independent observations xij ∼ N (0, Id ),

∗ yij = x⊤ ij θ + ζij ,

ζij ∼ N (0, σ 2 ),

with covariates independent of residuals. We choose k coordinates uniformly, set their coefficients to 1/k, and set the others to zero. The public defaults are n = 15000, m = 50, d = 256, k = 8, σ = 0.25, and ϵ = 4; each sweep varies one quantity. The realized support is used only for evaluation. Each setting has 30 repetitions in Exp1 and Exp2 and ten in Exp3 and Exp4, with shared training data and 4000 independent test observations across methods. Matching settings and privacy budgets reuse the data and randomization, with one Lasso fit per repetition. We simulate aggregate mechanism outputs and measure relative coefficient error ∥θ̂ − θ∗ ∥2 /∥θ∗ ∥2 , test MSE, support recall |Ŝ ∩ S ∗ |/k, and exact recovery. U2R uses disjoint user groups of sizes ⌊n/2⌋ for selection p and two equal remaining groups for coarse and refined estimation. Its local Lasso penalty is 0.01σ log(dn)/m, with coefficient threshold ρ/2 and count threshold n1 /(4k), where ρ = 1/k. The response probability is the larger of 32 log(1/δ)/(n1 ϵ2 ) and the value calibrated from shuffled binary response or the bound for two bits under pure local DP. Coordinates are shuffled independently. Stage 2 follows Algorithm 4 with local OLS, clipping radius and rounding coefficients set to one, and privacy calibration for finite samples, including the summation protocol. For 2-SLR, selection uses the 64 largest positive marginal correlations, selection size 2, B = 1, and 4 bins. LDP-IHT uses T = 2 and gradient radius 4; separate development seeds select η = 1.5, τx = 1.5, τy = 0.75, k ′ = 10 for ϵ < 4, and η = 2, τx = 1, τy = 0.5, k ′ = 8 otherwise. Sample-SDP uses c0 = 0.1 and c1 = 2. Lasso uses five fold cross validation over 300 penalties, no intercept, tolerance 10−4 , and at most 3000 iterations. These rules are fixed before the sweeps. Kelmarsh contains measurements every ten minutes from six turbines over 2017 to 2022. We predict nacelle temperature from 90 predictors: 59 physical mean signals, four cyclic time features, 22 missing value indicators, and five turbine indicators. Whole UTC dates are split into public, training, and test sets in nominal proportions 20%, 70%, and 10%, leaving 1,295,987 training and 185,416 test rows after removing missing responses. Separate public dates provide 80,000 rows for preprocessing and 40,000 for validation. Public medians fill missing values; public means and standard deviations scale features and responses, with indicators centered. Each real data repetition samples nm training rows without replacement and assigns them to n = 10000 simulated users. We test m ∈ {10, 50, 100} and ϵ ∈ {1, 4}, sharing rows and groups across methods and budgets. The test set is fixed. U2R uses ρ = 0.05, a unit parameter norm bound, and σ = 0.25. Public validation over all six settings and three separate seeds selects k = 12 from {8, 10, 12}. After count thresholding, U2R keeps at most min{k, m − 2} coordinates ranked by debiased private counts and calibrates Stage 2 to the selected dimension. For real data, 2-SLR uses 16 bins. LDP-IHT uses T = 5, η = 0.1, τx = 1, τy = 0.5, and k ′ = 8. √ e + 4 d sG , where sG is Sample-SDP uses c0 = 0.03, c1 = 1, and diagonal shift max{0, −λmin (G)} the noise standard deviation per coordinate. Lasso cross validation is grouped by UTC date. Other settings remain unchanged. We report ten repetitions of MSE in the original units and its ratio to the paired Lasso fit, reused across budgets. C.2

M AIN R ESULTS

Privacy budget and number of users. Figure 1 shows relative coefficient error as ϵ varies in Exp1 and n varies in Exp2. In Exp1, U2R has the lowest mean error among the private methods at every tested budget, decreasing from 0.04260 at ϵ = 0.5 to 0.00544 at ϵ = 32. At ϵ = 4, its error is 0.00631, comparable to non-private Lasso. Most of the reduction occurs between ϵ = 0.5 and 4, with smaller gains at larger budgets. In Exp2, increasing n from 8000 to 20000 at m = 50 reduces U2R’s mean error from 0.01022 to 0.00496. Each private method has lower mean error at the larger user count. U2R recovers the exact support in every repetition, so its remaining coefficient error lies on the active coordinates. Section C.3 examines support recovery as n and m vary. 19

Preprint

U2R

2-SLR

LDP-IHT

Sample-SDP

100.5

Relative coefficient error

Relative coefficient error

100.5

Lasso

100 10−0.5 10−1 10−1.5 10−2 10−2.5

100 10−0.5 10−1 10−1.5 10−2 10−2.5

0.5

1

2

4

ε

8

16

32

8

10

12

14

n/10

(a) Privacy budget.

3

16

18

20

(b) Number of users.

Figure 1: Relative coefficient error versus ϵ in the left, and n in the right, with default n = 15000, m = 50, d = 256, k = 8 and ϵ = 4. Lines show means over 30 repetitions; shaded bands span the empirical 95% coverage of 2.5th–97.5th percentiles.

Dimension d. In Exp3, U2R’s mean error ranges from 0.00553 to 0.00659 as d increases from 64 to 512, with k = 8 and m = 50 fixed, as shown in Table 6. It remains close to non-private Lasso across these dimensions. Sample-SDP’s error increases from 0.01849 to 0.14111 over the same range.

Table 6: Exp3: relative coefficient error versus dimension at n = 15000, m = 50, and ϵ = 4. Entries are means over ten repetitions. d

U2R

2-SLR

LDP-IHT

Sample-SDP

Lasso

64

0.00643

0.49079

0.13262

0.01849

0.00450

128

0.00659

0.61461

0.12595

0.03571

0.00453

256

0.00553

0.64914

0.13354

0.07027

0.00520

512

0.00598

0.77533

0.18692

0.14111

0.00522

Real data prediction. Table 7 reports test MSE relative to non-private Lasso in Exp5. U2R has the lowest mean ratio among the private methods in all six settings. As m increases from 10 to 100, its ratio decreases from 1.48261 to 1.12978 at ϵ = 1 and from 1.43505 to 1.13213 at ϵ = 4. More observations per user bring U2R closer to Lasso. Increasing ϵ reduces its mean prediction error at m = 10 and 50; at m = 100, the ratios are similar under both budgets. C.3

A DDITIONAL R ESULTS

Support recovery. Table 8 compares the Stage 1 selectors in Exp4. At n = 15000, U2R’s mean recall increases from 3.75% at m = 25 to 100% at m = 50 and remains 100% at m = 75 and 100. It recovers the exact support in all ten repetitions at each of these three larger values of m. The 2-SLR recall increases from 26.25% to 76.25% over the same sweep. With m = 50, U2R recovers the exact support in every repetition at all tested user counts from 8000 to 20000. The 2-SLR recall ranges from 35.00% to 63.75%. 20

Preprint

Table 7: Exp5: Kelmarsh test MSE divided by the paired non-private Lasso MSE at n = 10000 and d = 90. U2R and Sample-SDP use δ = 10−5 . Entries are mean ± sample standard deviation over ten repetitions. Budget

m 10

ϵ=1

50 100 10

ϵ=4

50 100

U2R

2-SLR

LDP-IHT

Sample-SDP

1.48261 ±0.04759 1.17245 ±0.02983 1.12978 ±0.02725

11.21592 ±10.06853 3.08118 ±1.14883 3.23649 ±1.62979

3.15212 ±0.72151 2.36752 ±0.39188 2.25214 ±0.45721

3.14939 ±0.36992 1.91000 ±0.31047 1.50741 ±0.12895

1.43505 ±0.01006 1.15861 ±0.02794 1.13213 ±0.00306

3.07791 ±1.29223 2.08907 ±0.16864 2.19689 ±0.26004

2.44676 ±0.59026 2.36333 ±0.41893 2.32449 ±0.33028

2.05930 ±0.36280 1.48151 ±0.27440 1.37230 ±0.12374

Table 8: Exp4: Stage 1 support recall and exact recovery. Recall entries are mean ± sample standard deviation in percentage points; exact recovery is the number of repetitions with Ŝ = S ∗ out of ten. Value

U2R recall in %

2-SLR recall in %

U2R exact

2-SLR exact

Exp4a: vary m; n = 15000, d = 256, ϵ = 4 25

3.75 ± 8.44

26.25 ± 13.76

0/10

0/10

50

100.00 ± 0.00

63.75 ± 12.43

10/10

0/10

75

100.00 ± 0.00

72.50 ± 16.46

10/10

1/10

100

100.00 ± 0.00

76.25 ± 12.43

10/10

0/10

Exp4b: vary n; m = 50, d = 256, ϵ = 4 8000

100.00 ± 0.00

35.00 ± 12.91

10/10

0/10

10000

100.00 ± 0.00

46.25 ± 10.29

10/10

0/10

12000

100.00 ± 0.00

50.00 ± 10.21

10/10

0/10

15000

100.00 ± 0.00

63.75 ± 12.43

10/10

0/10

18000

100.00 ± 0.00

53.75 ± 15.65

10/10

0/10

20000

100.00 ± 0.00

61.25 ± 12.43

10/10

0/10

21

Preprint

Prediction error. At the default setting, U2R and Lasso have mean test MSEs of 0.0625570 and 0.0625521, respectively, close to the residual variance σ 2 = 0.0625. LDP-IHT and Sample-SDP yield 0.0648534 and 0.0631971. At ϵ = 32, U2R’s MSE is 0.0625546; its reduction in coefficient error has little effect on prediction MSE over this budget range. Prediction in original units. Table 9 gives test MSE in squared degrees Celsius for Exp5. At ϵ = 1, increasing m from 10 to 100 reduces U2R’s mean MSE from 3.55270 to 2.70655 and mean RMSE from 1.88464 to 1.64505 ◦ C; mean R2 increases from 0.87900 to 0.90782. At m = 100 and ϵ = 4, the mean MSEs of U2R, Lasso, and Sample-SDP are 2.71218, 2.39565, and 3.28753, respectively. Table 9: Exp5: Kelmarsh test MSE in squared degrees Celsius. U2R and Sample-SDP use δ = 10−5 . Entries are mean ± sample standard deviation over ten repetitions; Lasso is reused across budgets. Budget

m 10

ϵ=1

50 100 10

ϵ=4

50 100

D

U2R

2-SLR

LDP-IHT

Sample-SDP

Lasso

3.55270 ±0.11390 2.80885 ±0.07139 2.70655 ±0.06479

26.89088 ±24.15804 7.38104 ±2.75047 7.75405 ±3.90559

7.55337 ±1.72935 5.67181 ±0.93834 5.39543 ±1.09591

7.54559 ±0.87616 4.57575 ±0.74360 3.61122 ±0.30899

2.39625 ±0.00376 2.39572 ±0.00123 2.39565 ±0.00063

3.43873 ±0.02432 2.77572 ±0.06713 2.71218 ±0.00678

7.37791 ±3.10457 5.00484 ±0.40439 5.26299 ±0.62298

5.86216 ±1.40887 5.66193 ±1.00425 5.56877 ±0.79226

4.93425 ±0.86685 3.54916 ±0.65645 3.28753 ±0.29627

2.39625 ±0.00376 2.39572 ±0.00123 2.39565 ±0.00063

O MITTED P ROOFS IN S ECTION 4

This section provides proofs for the theories in Section 4 of the main text. D.1

P ROOF OF T HEOREM 1

We first reduce a robust shuffled DP protocol to a pan-private algorithm, and then prove the lower bound for pan-private algorithms by a hybrid argument over the stream and Assouad’s method over the coordinates. Let X = {±1}d+1 be the data universe of Definition 6, let U be the uniform distribution on X , and let χi (x, y) = yxi for i ∈ [d]. These functions satisfy EU [χi ] = 0 and EU [χi χl ] = I[i = l]. For
θ ∈ Rd let Qθ (x, y) = U (x, y) 1 + y⟨x, θ⟩ , which is the distribution of one sample in Definition 6 with parameter θ. For z ∈ {−1, 0, 1}d we write Pz = Qγz and Z0 (z) = {i ∈ [d] : zi = 0}. For i ∈ Z0 (z) and b ∈ {±1}, let z i,b = z + bei . Then Pzi,b − Pz = bγ χi U.

(21)

For a distribution P on X , M (P n ) denotes the distribution of the output of M on n i.i.d. records from P . Following Cheu & Ullman (2021), a protocol P is (ϵ, δ, 1/3)-robust shuffled DP if, for every γ ′ ≥ 1/3, the shuffled messages of any γ ′ n users that follow the protocol satisfy (ϵ, δ)-DP; this is the robustness in Theorem 1. Let P be such a non-interactive protocol and let n′ = n/3. By Cheu & Ullman (2021, Theorem 2.6), there is an (ϵ, δ)-pan-private algorithm M P that processes a stream of n′ records and satisfies, for every distribution P on X ,
′ n dT V M P (P n ), MP (P(2/9) ) ≤ e−c1 n , P(b) = bP + (1 − b)U, where c1 > 0 is an absolute constant. Since Qθ is linear in θ, we have 92 Qθ + 79 U = Q2θ/9 , e k whenever θ ∈ Θ e k . Let Π1 be the projection onto the unit ℓ1 -ball and define and 2θ/9 ∈ Θ 22

Preprint

M = 92 Π1 (M P ). As a post-processing of the output, M is (ϵ, δ)-pan-private. The projection does e not increase the distance to any point of the unit ℓ1 -ball, and ∥M − θ∥2 ≤ 11 2 for θ ∈ Θk , so   9   11 ′ E ∥M (Qnθ ) − θ∥2 ≤ E ∥MP (Qn2θ/9 ) − 29 θ∥2 + e−c1 n 2 2

e k. for all θ ∈ Θ

It therefore suffices to show that every (ϵ, δ)-pan-private algorithm M on n′ records satisfies √   dk ∗ n′ sup E ∥M (Qθ∗ ) − θ ∥2 ≥ c2 ′ nϵ ek θ ∗ ∈Θ

(22)

(23)

for a constant c > 0. Indeed, combining equation 22√and equation 23 gives supθ∗ E∥MP (Z) − √ 2 dk dk 2 ∗ −c1 n θ ∥2 ≥ 3 c2 nϵ − 11 . Since ϵ ≤ O(1), we have nϵ ≥ Ω( n1 ), so for n larger than an absolute 9 e √

dk constant the last term is at most 13 c2 nϵ , which gives the theorem. In the rest of the proof we fix an (ϵ, δ)-pan-private algorithm M with internal states S1 , . . . , Sn′ , and prove equation 23. We take n larger than an absolute constant, so that δ = o(1/n) gives n′ δ ≤ 1/64.

For j ∈ [n′ ], a prefix x<j = (x1 , . . . , xj−1 ) and a record x ∈ X , let qj (· | x<j , x) be the distribution of Sj when the first j records are (x<j , x). Applying Definition 5 at time t = j to two streams that differ in the j-th record shows that, for every fixed x<j , the map x 7→ qj (· | x<j , x) is (ϵ, δ)DP. By Cheu & Ullman (2021, Lemma 3.7), there is a (2ϵ, 0)-DP map x 7→ qj′ (· | x<j , x) with
dT V qj (· | x<j , x), qj′ (· | x<j , x) ≤ δ for all x. Below we drop x<j from the notation when it is fixed, and write sums over states s; for general state spaces the sums are replaced by integrals of densities with respect to a common dominating measure. Let q̄j′ (s) = EX∼U [qj′ (s | X)] and hj,s (x) =

qj′ (s | x) − 1, q̄j′ (s)

  ĥj,s (i) = EU hj,s χi .

Since qj′ is (2ϵ, 0)-DP, e−2ϵ q̄j′ (s) ≤ qj′ (s | x) ≤ e2ϵ q̄j′ (s), so |hj,s (x)| ≤ e2ϵ − 1. For a distribution P on X , write qj ◦ P for the distribution of Sj when the j-th record is drawn from P . If P ′ − P = bγχi U , then using EU [χi ] = 0,
1X X ′ qj (s | x) P ′ (x) − P (x) 2 s x γX γ = EU [qj′ (s | X)χi (X)] = Es∼q̄j′ ĥj,s (i) . 2 s 2

dT V (qj′ ◦ P ′ , qj′ ◦ P ) =

(24)

Fix z ∈ {−1, 0, 1}d , i ∈ Z0 (z) and b ∈ {±1}. For j = 0, . . . , n′ , let Hj be the distribution of the output of M when X1 , . . . , Xj ∼ Pz and Xj+1 , . . . , Xn′ ∼ Pzi,b , all independent. Then ′ ′ H0 = M (Pzni,b ) and Hn′ = M (Pzn ). In Hj−1 and Hj the prefix X<j has the same distribution ′

Pzj−1 , the suffix X>j has the same distribution Pzni,b−j and is independent of (X<j , Xj ), and only Xj differs. Since M processes the stream one record at a time, its output is a post-processing of (Sj , X>j ). Therefore hγ i
dT V (Hj−1 , Hj ) ≤ EX<j ∼Pzj−1 dT V qj ◦ Pzi,b , qj ◦ Pz ≤ EX<j ∼Pzj−1 Es∼q̄j′ ĥj,s (i) + 2δ, 2 where the second step replaces qj by qj′ at a cost of 2δ and then applies equation 21 and equation 24. By the triangle inequality, ′

n ′ ′
γX ′ j−1 Es∼q̄ ′ ĥj,s (i) + 2n δ. dT V M (Pzni,b ), M (Pzn ) ≤ E j 2 j=1 X<j ∼Pz

(25)

The distribution of X<j and the maps qj′ in equation 25 depend on z but not on i or b. P Since χ1 , . . . , χd are orthonormal in L2 (U ), Bessel’s inequality gives i∈[d] ĥj,s (i)2 ≤ EU [h2j,s ] ≤ √ P (e2ϵ − 1)2 , and by the Cauchy–Schwarz inequality i∈Z0 (z) |ĥj,s (i)| ≤ d (e2ϵ − 1). Summing 23

Preprint

equation 25 over i ∈ Z0 (z) and b ∈ {±1}, and exchanging the sum over i with the expectations, which do not depend on i, we obtain X X √ ′ ′
T (z) := Ti (z) ≤ γn′ d (e2ϵ −1)+4dn′ δ, Ti (z) := dT V M (Pzni,b ), M (Pzn ) . i∈Z0 (z)

b∈{±1}

(26) Note that equation 26 bounds the sum over coordinates only, and this is all that is used below. ′

Draw Z as in Definition 6, let θZ = γZ, draw the stream X ∼ PZn given Z, and let θ̂ = M (X). For i ∈ [d] let Gi = {∥Z−i ∥0 ≤ k − 1} and ψi = sign(θ̂i ), with sign(0) = +1. Condition on Z−i ′ and on Zi ̸= 0, and let z = (Z−i , 0), so that zi = 0. Then Zi is uniform on {±1} and X ∼ Pzni,Zi , so by Le Cam’s method and the triangle inequality,
 1

1 ′ ′ P ψi ̸= Zi | Z−i , Zi ̸= 0 ≥ 1 − dT V M (Pzni,+1 ), M (Pzni,−1 ) ≥ 1 − Ti (z) . 2 2 P Let N = i∈[d] I[Zi ̸= 0, Gi , ψi ̸= Zi ]. If N ≥ 1, some i has Zi ̸= 0 and ∥Z−i ∥0 ≤ k − 1, so ∥Z∥0 ≤ k and N ≤ ∥Z∥0 ≤ k. Each counted coordinate has |θ̂i − γZi | ≥ γ. Hence √ γN ∥θ̂ − θZ ∥2 I[∥Z∥0 ≤ k] ≥ γ N ≥ √ . k

(27)

k , Since Zi is independent of Z−i and P(Zi ̸= 0) = 2d

E[N ] ≥

d 
 k X P(Gi ) − E Ti (Z−i , 0) . 4d i=1

k ) has mean less than k2 , so Markov’s We bound the two sums in turn. First, ∥Z−i ∥0 ∼ Bin(d − 1, 2d P inequality gives P(Gci ) ≤ 21 and i P(Gi ) ≥ d2 . Second, (Z−i , 0) = Z on the event {Zi = 0}, k which is independent of Z−i and has probability 1 − 2d ≥ 21 , so by equation 26,   d d X √   
 X E Ti (Z) I[Zi = 0] ≤ 2 E T (Z) ≤ 2γn′ d (e2ϵ − 1) + 8dn′ δ. E Ti (Z−i , 0) = P(Zi = 0) i=1 i=1

Combining the two bounds, E[N ] ≥

k kγn′ (e2ϵ − 1) √ − − 2kn′ δ. 8 2 d

√

k 1 , and use n′ δ ≤ 64 , so that the We now choose γ = 8n′ (e2ϵd −1) , so that the second term equals 16 k k third term is at most 32 . This gives E[N ] ≥ 32 , and by equation 27, √ √     γ k ′ dk sup E ∥M (Qnθ∗ ) − θ∗ ∥2 ≥ E ∥θ̂ − θZ ∥2 I[∥Z∥0 ≤ k] ≥ = . ′ (e2ϵ − 1) 32 256 n ∗ ek θ ∈Θ

Finally, for 0 < ϵ ≤ c0 with an absolute constant c0 , we have e2ϵ − 1 ≤ 2e2c0 ϵ, which gives equation 23 with c2 = 512e1 2c0 and completes the proof. D.2

P ROOF OF T HEOREM 2

We first give a lemma on the binomial mechanism. It is used in the privacy part of this proof, and again in the proofs of Proposition 1 and Theorem 3. √ Lemma 1. Let s ≥ 1, ϵ > 0 and δ ∈ (0, 1) with ϵ ≤ log(2s/δ), and let ∆ ≥ s ϵ. Let B = (B1 , . . . , Bs ) have independent entries Bj ∼ Bin(M, 1/2), where 32∆2 log(2s/δ) . ϵ2 Then for any v, v ′ ∈ Zs with ∥v − v ′ ∥2 ≤ ∆ and any E ⊆ Zs , M≥

P[v + B ∈ E] ≤ eϵ P[v ′ + B ∈ E] + δ. 24

(28)

Preprint

√ Only an upper bound on ∥v −v ′ ∥2 is used, so the condition ∆ ≥ s ϵ can always be met by taking ∆ larger. Subtracting a constant from v +B, for example M/2 on each coordinate, is a post-processing and keeps the guarantee. Proof of Lemma 1. Let a = v ′ − v and A = M/2. For t ∈ {0, . . . , M }s define ℓ(t) =

s X j=1

log

P[Bj = tj ] , P[Bj = tj − aj ]

with ℓ(t) = +∞ if tj − aj ∈ / {0, . . . , M } for some j. For every output o, P[v + B = o] = ℓ(o−v) ′ e P[v + B = o], so P[v + B ∈ E] ≤ eϵ P[v ′ + B ∈ E] + P[ℓ(B) > ϵ], and it suffices to show P[ℓ(B) > ϵ] ≤ δ. p Write wj = Bj − A, let W = M log(4s/δ)/2, and let G = {maxj |wj | ≤ W }. Each wj is a sum of M independent terms in [−1/2, 1/2], so Hoeffding’s inequality gives P[|wj | > W ] ≤ √ 2 2e−2W /M = δ/(2s), and P[Gc ] ≤ δ/2. We check that W + ∆ ≤ A/2. From ∆ ≥ s ϵ and equation 28, M ≥ 32s √ log(2s/δ) ≥ 32 log(2s/δ). Since log(4s/δ) ≤ 2 log(2s/δ), we get p W ≤ M log(2s/δ) ≤ M/ 32. From equation 28 and ϵ ≤ log(2s/δ), r √ ϵ M M log(2s/δ) M ≤ ∆≤ p ≤ . 32 32 32 log(2s/δ) Hence W + ∆ ≤ M/4 = A/2. Fix a coordinate j with aj > 0 and drop the index j. On G, B = A + w and
M a Y A−w+r P[B = A + w] A+w = = .
M P[B = A + w − a] A +w−a+r A+w−a r=1 Let xr = (w −r)/A and yr = (w −a+r)/A. Since |w|+a ≤ A/2, we have xr < 1 and |yr | ≤ 1/2. Using log(1 − x) ≤ −x and − log(1 + y) ≤ −y + y 2 for |y| ≤ 1/2, and |w − a + r| ≤ |w| + a for 1 ≤ r ≤ a, log

a − 2w (|w| + a)2 A−w+r = log(1 − xr ) − log(1 + yr ) ≤ + . A+w−a+r A A2

Summing over r gives the bound (a2 − 2aw)/A + a(|w| + a)2 /A2 for this coordinate. If aj < 0, the same bound applied to M −Bj , which has the same distribution, gives (a2j −2aj wj )/A+|aj |(|wj |+ |aj |)2 /A2 ; if aj = 0 the term is zero. Summing over j and using A = M/2, on G we obtain ℓ(B) ≤

s
2 2∥a∥22 4⟨a, w⟩ 4 X − + 2 |aj | |wj | + |aj | . M M M j=1

We bound the three terms. For the first, by equation 28 and ϵ ≤ log(2s/δ), 2∥a∥22 2∆2 ϵ2 ϵ ≤ ≤ ≤ . M M 16 log(2s/δ) 16 P √ For the third, (|wj | + |aj |)2 ≤ 2W 2 + 2∆2 on G, j |aj | ≤ s ∆ and W 2 ≤ M log(2s/δ), so √ with equation 28, ∆ ≥ s ϵ and ϵ ≤ log(2s/δ), √ √ s
2 4 X 8 s ∆ log(2s/δ) 8 s ∆3 |a | |w | + |a | ≤ + j j j M 2 j=1 M M2 √ 2 √ 4 sϵ sϵ ϵ ϵ ≤ + ≤ + . 4∆ 4 128 128 ∆ log2 (2s/δ) 25

Preprint

For the second, ⟨a, w⟩ is a sum of sM independent terms, where the terms of coordinate j lie in an interval of length |aj |. By Hoeffding’s inequality and equation 28, h 4⟨a, w⟩  2(ϵM/8)2   ϵ2 M  δ ϵi δ P − ≤ exp − ≤ exp − ≤ > ≤ . M 2 M ∥a∥22 32∆2 2s 2 On G and outside this event, ℓ(B) ≤ ϵ/16 + ϵ/2 + ϵ/4 + ϵ/128 < ϵ. Hence P[ℓ(B) > ϵ] ≤ δ/2 + δ/2 = δ. Proof of Theorem 2. The parameters are set as in Table 3. We use n ≥ 20, which the summation step requires, and ϵ ≤ log(4D/δ), which holds for p d larger than an absolute constant since ϵ ≤ O(1). √ By the choice of p and CF , pCF ≥ c0 n ϵ/ log(4D/δ). Since ϵ ≤ O(1), the third term of equation 8, with a large enough constant, then gives √ √ (p) pCF ≥ D max{1, ϵ}, so D ϵ ≤ ∆2 ≤ 3pCF . (29) We first prove privacy. The rounding step of Algorithm 2 has the same distribution as Ūi,j = ⌊pF̄i,j +ui,j ⌋ with ui,j independent and uniform on (0, 1], and we use this form. Fix two neighboring datasets that differ in the sample of user i0 , and fix u = (ui,j ). The rounded vectors of all users other than i0 are the same for the two datasets. For user i0 , since |⌊a + u⌋ − ⌊b + u⌋| ≤ |a − b| + 1 and ∥F̄i0 − F̄i′0 ∥2 ≤ 2CF , √ √ (p) ∥Ūi0 − Ūi′0 ∥2 ≤ p∥F̄i0 − F̄i′0 ∥2 + D ≤ 2pCF + D = ∆2 . Pn Pn Let Tj = i=1 (Ūi,j + ζi,j ) for j ∈ [D]. Then i=1 ζi,j = Bj − nN0 /2 with Bj ∼ Bin(nN0 , 1/2) (p) independent across j, and nN0 ≥ M ≥ 32(∆2 )2 log(4D/δ)/ϵ2 . By equation 29 and ϵ ≤ (p) log(4D/δ), Lemma 1 applies with s = D, ∆ = ∆2 and δ/2 in place of δ. Hence, for fixed u, the vector T = (T1 , . . . , TD ) is (ϵ, δ/2)-DP as a function of the dataset. The variables u have the same distribution for both datasets, so T is also (ϵ, δ/2)-DP without conditioning on u. We next relate the view of the aggregator to T . Since |Ūi,j | ≤ pCF + 1, |ζi,j | ≤ N0 /2 and N0 ≥ 2, we have |Tj | ≤ n(pCF + 1 + N0 /2) < q/2. So no wrap-around occurs, and Tj is the value Sj computed by the aggregator. On coordinate j, the users run the summation protocol Ver,n of Balle et al. (2020) on the inputs W1,j , . . . , Wn,j ∈ Zq , with r shuffled shares and one more share. Since q > 4n, the setting of r gives r ≥ 3, and n ≥ 20. By Balle et al. (2020, Corollary 6.1), the protocol has worst-case statistical security with parameter σsec = log2 (8Deϵ /δ): for two input vectors with the same sum in Zq , the total variation distance between the messages of coordinate j is at most 2−σsec = e−ϵ δ/(8D). Let V sim be the view obtained when, on every coordinate j, the inputs are replaced by (Tj mod q, 0, . . . , 0). It is a post-processing of T , hence (ϵ, δ/2)-DP. The D coordinates use independent shares and separate shufflers, so coupling them one at a time gives, for every dataset,
e−ϵ δ e−ϵ δ = . dT V VP , V sim ≤ D · 8D 8 By Balle et al. (2020, Lemma F.2), VP is (ϵ, δ/2 + (1 + eϵ )e−ϵ δ/8)-DP. Since (1 + eϵ )e−ϵ /8 ≤ 1/4, Algorithm 2 satisfies (ϵ, δ)-shuffled DP. 2 We now bound the error. Under Assumption 3, ∥v(Hi )∥2 ≤ ∥xi x⊤ i ∥F = ∥xi ∥2 , so √ ∥Fi ∥2 ≤ ∥xi ∥22 + ∥xi ∥2 |yi | ≤ c2x d + cx cy d ≤ CF

for c0 ≥ c2x + cx cy and n ≥ 3. Hence F̄i = Fi for every i. The aggregator computes F̂ = T /(np), and n

F̂ =

1X Fi + E rd + E bin , n i=1

Ejrd =

Let η=

CF

n
1 X Ūi,j − pFi,j , np i=1

p log(4D/δ) log D . nϵ 26

Ejbin =

Bj − nN0 /2 . np

Preprint

The terms Ūi,j −pFi,j are independent, have mean zero and lie in (−1, 1]. By Hoeffding’s inequality and a union bound over j ∈ [D], with probability at least 1 − O(d−2 ),  √log D  rd ∥E ∥∞ ≤ O √ ≤ O(η), np √ √ d log n log(4D/δ) CF log(4D/δ) = . Each Bj − nN0 /2 is a where the second step uses √1n p ≤ nϵ c0 nϵ √ sum of p nN0 independent M + 2n. By equation 29, M ≤
terms in [−1/2, 1/2], and nN0 ≤ −2 O pCF log(4D/δ)/ϵ , so with probability at least 1 − O(d ), √   p(M + 2n) log D  log D  bin ≤O η+ √ = O(η). ∥E ∥∞ ≤ O np np Pn Pn 1 Let H̄ = n1 i=1 xi x⊤ i and ḡ = n i=1 xi yi . Since Ĥ and ĝ are read from F̂ , we obtain ∥Ĥ − H̄∥∞,∞ ≤ O(η) and ∥ĝ − ḡ∥∞ ≤ O(η). Pn 1 2 ∗ For the statistical error, the entries of xi x⊤ i are bounded by cx , and ḡ − H̄θ = n i=1 xi ζi , where, given xi , each entry xiℓ ζi is sub-Gaussian with mean zero and variance at most c2x σ 2 by Assumption 1. By Hoeffding’s inequality and a union bound, with probability at least 1 − O(d−2 ), r  r log d   log d  2 ∗ ∥H̄ − Σ∥∞,∞ ≤ O cx , ∥ḡ − H̄θ ∥∞ ≤ O cx σ . n n For any w ∈ Rd , ∥(Ĥ − Σ)w∥∞p ≤ d ∥Ĥ − Σ∥∞,∞ ∥w∥∞ . The first two terms of equation 8, with large enough constants, give d c2x log d/n ≤ κ∞ /4 and dη ≤ κ∞ /4, so by Assumption 3, κ∞ ∥Ĥw∥∞ ≥ ∥Σw∥∞ − ∥(Ĥ − Σ)w∥∞ ≥ ∥w∥∞ . 2 Hence Ĥ is invertible and ∥Ĥ −1 z∥∞ ≤ κ2∞ ∥z∥∞ for all z. Since ĝ − Ĥθ∗ = (ĝ − ḡ) − (Ĥ − H̄)θ∗ + (ḡ − H̄θ∗ ) and ∥(Ĥ − H̄)θ∗ ∥∞ ≤ ∥Ĥ − H̄∥∞,∞ ∥θ∗ ∥1 ≤ O(η), we get r   2 log d ∥θ̃ − θ ∥∞ = ∥Ĥ (ĝ − Ĥθ )∥∞ ≤ + η ≤ λn O cx σ κ∞ n √ for a large enough c1 . By Zhu et al. (2024, Lemma 29), ∥Sλn (θ̃) − θ∗ ∥2 ≤ 4 k λn . The projection onto {θ : ∥θ∥1 ≤ 1} does not increase the distance to θ∗ , which lies in this set. Hence, with κ∞ a constant,  r k log d d√k log nplog(4D/δ) log D  √ ∗ ∥θ̂ − θ ∥2 ≤ 4 k λn ≤ O σ + , n nϵ ∗

−1

∗

which gives the rate of Theorem 2. The events above hold together with probability at least 1 − O(d−2 ).

E

O MITTED P ROOFS IN S ECTION 5

This section provides the rigorous proofs for the theoretical bounds in Section 5. E.1

P ROOF OF P ROPOSITION 1

Proof of Proposition 1. Throughout the proof we take k = |S ∗ |. The parameters of Algorithm 3 are set as in Table 4, in particular γ = Cγ log(1/δ)/(n1 ϵ2 ). We use ϵ ≤ ϵ0 for an absolute constant ϵ0 , and take δ ≤ 1/8 small enough that ϵ0 /2 ≤ log(8/δ). By Assumption 4 with v = eℓ for ℓ ∈ S ∗ , we have κ ≤ Σℓℓ ≤ 1. Since k log(ed/k) ≥ log d, the second term of equation 10 gives m ≥ Ω(log d). The proof has three parts: the local Lasso of one user of Stage 1 selects S ∗ with probability at least 1−1/(32k); Stage 1 is (ϵ, δ)-shuffled DP; and the debiased counts separate S ∗ from its complement. 27

Preprint

The first part is about one user, so it needs no union bound over the users, and a user with Ŝi ̸= S ∗ is allowed. The privacy part is not needed for the statement of the proposition, and is used in the proof of Theorem 3. We first show that P[Gci ] ≤ 1/(32k) for every i ∈ [n1 ], where Gi is the event Ŝi = S ∗ . Fix i ∈ [n1 ] and let ζi ∈ Rm be the noise of user i. For ℓ ∈ [d], write xℓ = ⟨Σ1/2 eℓ , Σ−1/2 x⟩. By Assumption 2 and Σℓℓ ≤ 1, we have ∥xℓ ∥ψ2 ≤ κx . By Assumption 1, ∥ζ∥ψ2 ≤ O(σ) and E[xℓ ζ] = E[xℓ E[ζ | 1 x]] = 0. So each entry of m Xi⊤ ζi is an average of m independent sub-exponential variables with mean zero and ψ1 -norm at most O(κx σ). Since k ≤ d, we have log(128dk) ≤ 2 log(128d). Bernstein’s inequality, m ≥ Ω(log d) and a union bound over ℓ ∈ [d] give, with probability at least 1 − 1/(64k), r  1 ⊤ log d  λloc ≤ O κx σ ≤ Xi ζi , (30) m m 2 ∞ where the last step holds for cλ large enough. Let C = {v ∈ Rd : ∥v(S ∗ )c ∥1 ≤ 3∥vS ∗ ∥1 }. The rows of Xi are x = Σ1/2 (Σ−1/2 x) with Σ−1/2 x isotropic and sub-Gaussian. By Assumption 4, maxℓ Σℓℓ ≤ 1 and the result of Rudelson & Zhou (2013) for such designs, there are constants such that, when m ≥ Ω(k log(ed/k)/κ), with probability at least 1 − 2e−cm , 1 1 κ ∥Xi v∥22 ≥ v ⊤ Σv ≥ ∥v∥22 m 2 2

for all v ∈ C.

Since κ ≤ 1 and k log(ed/k) ≥ k, the second term of equation 10 gives the required m and 2e−cm ≤ 1/(64k). So the two events hold together with probability at least 1 − 1/(32k). On these two events, let ∆ = θ̂(i) −θ∗ . Since θ̂(i) minimizes the Lasso objective and yi = Xi θ∗ +ζi ,
1 1 ∥Xi ∆∥22 ≤ ⟨Xi⊤ ζi , ∆⟩ + λloc ∥θ∗ ∥1 − ∥θ∗ + ∆∥1 2m m 3λloc λloc ≤ ∥∆S ∗ ∥1 − ∥∆(S ∗ )c ∥1 , 2 2 where we used equation 30 and ∥θ∗ ∥1 − ∥θ∗ + ∆∥1 ≤ ∥∆S ∗ ∥1 − ∥∆(S ∗ )c ∥1 . The left-hand side is non-negative, so ∆ ∈ C, and then κ 1 3λloc √ ∥∆∥22 ≤ ∥Xi ∆∥22 ≤ k ∥∆∥2 . 4 2m 2 √ √ Hence ∥∆∥∞ ≤ ∥∆∥2 ≤ 6λloc k/κ ≤ 3ρ/8, where the last step uses λloc ≤ κρ/(16 k), which (i) holds under the first term of equation 10. For ℓ ∈ S ∗ , |θ̂ℓ | ≥ ρ − 3ρ/8 > ρ/2, and for ℓ ∈ / S∗, (i) ∗ c |θ̂ℓ | ≤ 3ρ/8 < ρ/2. So Ŝi = S , and therefore P[Gi ] ≤ 1/(32k). We next show that Stage 1 satisfies (ϵ, δ)-shuffled DP. For a pair of distributions, let Heϵ (P ∥ Q) = supE (P (E) − eϵ Q(E)); a mechanism is (ϵ, δ)-DP if Heϵ between its outputs on any two neighboring datasets, in both orders, is at most δ. We use that Heϵ is jointly convex, and that for product distributions Heϵ1 +ϵ2 (P1 × P2 ∥ Q1 × Q2 ) ≤ δ1 + δ2 whenever Heϵj (Pj ∥ Qj ) ≤ δj for j = 1, 2 (Dwork & Roth, 2014). Fix two datasets that differ in the local data of user i. We proceed in five steps. 1. P The ℓ-th shuffler outputs a multiset of bits, which is determined by the count Cℓ = n1 j=1 rjℓ . So the view is (C1 , . . . , Cd ). Given v = (v1 , . . . , vn1 ), the bits rjℓ are independent, and so are C1 , . . . , Cd . 2. The indices vj with j ̸= i have the same distribution under the two datasets, while vi ∼ Unif(Ŝi ) and vi′ ∼ Unif(Ŝi′ ) may differ. By joint convexity, it suffices to bound Heϵ between the views given (v−i , vi = a) and (v−i , vi = b), for any fixed v−i and a, b ∈ [d]. If a = b the two views are equal, so we take a ̸= b. 3. Given (v−i , a) and (v−i , b), the counts Cℓ with ℓ ∈ / {a, b} have the same distribution, and all counts are independent. By the composition bound above, it suffices to show that Ca , and in the same way Cb , is (ϵ/2, δ/2)-indistinguishable between the two cases. 28

Preprint

4. Consider Ca . Write ria = I[vi = a] ⊕ β with β ∼ Ber(γ), which has the same distribution in both cases; by joint convexity we fix β. For j ̸= i, since γ ≤ 1/2, the bit rja = I[vj = a] ⊕ Ber(γ) has the same distribution as a uniform bit with probability 2γ and I[vj = a] otherwise. Let J ⊆ [n1 ] \ {i} be the users that output a uniform bit, so that K = |J| ∼ Bin(n1 − 1, 2γ) does not depend on the data. Given J, X Ca = I[vj = a] + ria + Bin(K, 1/2), j ∈J, / j̸=i

and the integer part before Bin(K, 1/2) changes by one between the two cases. Lemma 1 with s = 1, ∆ = max{1, ϵ/2}, ϵ/2 in place of ϵ and δ/4 in place of δ applies, since ϵ/2 ≤ log(8/δ). It shows that, given J, Ca is (ϵ/2, δ/4)-indistinguishable as long as K≥

128 max{1, ϵ2 /4} log(8/δ) . ϵ2

(31)

By joint convexity over J, Heϵ/2 for Ca is at most P[equation 31 fails] + δ/4. 5. Since n1 ≥ 8 by equation 11, E[K] = 2γ(n1 − 1) ≥ 74 γn1 , and a Chernoff bound gives P[K < E[K]/2] ≤ e−E[K]/8 ≤ e−7γn1 /32 . We take Cγ = 256 max{1, ϵ20 /4}. With δ ≤ 1/8 we have log(8/δ) ≤ 23 log(1/δ) and log(4/δ) ≤ 34 log(1/δ), so γn1 = Cγ log(1/δ)/ϵ2 gives both E[K]/2 ≥ 78 γn1 , which is at least the right-hand side of equation 31, and e−7γn1 /32 ≤ δ/4. Hence Ca is (ϵ/2, δ/2)-indistinguishable, and the same holds for Cb . By steps 2 and 3, Stage 1 satisfies (ϵ, δ)-shuffled DP. We finally show that Ŝ = S ∗ . For ℓ ∈ [d], let qℓ = P[vi = ℓ] and wℓ = |{i ∈ [n1 ] : vi = ℓ}|. The local datasets of the users of Stage 1 are i.i.d., and each user draws vi with its own randomness, so v1 , . . . , vn1 are i.i.d. and wℓ ∼ Bin(n1 , qℓ ). On Gi , vi is uniform on S ∗ ; off Gi , vi may be any index. So for ℓ ∈ S ∗ and ℓ′ ∈ / S∗, qℓ ≥

31 P[Gi ] ≥ , k 32k

  1 qℓ′ = P vi = ℓ′ , Gci ≤ P[Gci ] ≤ . 32k

Let L = log(dn). For ℓ ∈ S ∗ , wℓ dominates Bin(n1 , 31/(32k)), and for ℓ′ ∈ / S ∗ , wℓ′ is dominated by Bin(n1 , 1/(32k)). Chernoff bounds give h h n1 i n1 i P wℓ < ≤ e−n1 /(9k) , P wℓ′ > ≤ e−n1 /(18k) . 2k 8k By a union bound over [d] and n1 ≥ Ω(kL) from equation 11, with probability at least 1 − O(n−2 ), n1 n1 wℓ ≥ for all ℓ ∈ S ∗ , wℓ′ ≤ for all ℓ′ ∈ / S∗. 2k 8k Given v, E[Cℓ | v] = n1 γ + (1 − 2γ)wℓ , so Ŵℓ − wℓ = (Cℓ − E[Cℓ | v])/(1 − 2γ). The count Cℓ is a sum of n1 independent bits, each with variance at most γ. By Bernstein’s inequality and a union bound over ℓ ∈ [d], with probability at least 1 − O(n−2 ), p for all ℓ ∈ [d]. Cℓ − E[Cℓ | v] ≤ 6n1 γL + 2L We have γ ≤ 1/4: this is n1 ≥ 4Cγ log(1/δ)/ϵ2 , which follows from equation 11 when δ ≥ 2 2 (dn)−Ck ϵ , as in Appendix B.2. So 1/(1 − 2γ) ≤ 2, and p n1 n1 n1 max Ŵℓ − wℓ ≤ 2 6n1 γL + 4L ≤ + = . 32k 32k 16k ℓ∈[d] p For the first term we used n1 γ = Cγ log(1/δ)/ϵ2 and n1 ≥ Ω(k L log(1/δ)/ϵ); for the second, n1 ≥ Ω(kL); both are in equation 11. Then for ℓ ∈ S ∗ , Ŵℓ ≥ n1 /(2k)−n1 /(16k) = 7n1 /(16k) > τHH , and for ℓ′ ∈ / S ∗ , Ŵℓ′ ≤ n1 /(8k) + n1 /(16k) = 3n1 /(16k) < τHH , where τHH = n1 /(4k). Therefore Ŝ = S ∗ and s = |Ŝ| = k. Combining the events, the statement holds with probability at least 1 − O(n−2 ). 29

Preprint

E.2

P ROOF OF T HEOREM 3

Proof of Theorem 3. The parameters are set as in Table 4 and Table 5. We use n3 , n4 ≥ 19, which the summation step requires, and ϵ ≤ log(4s/δ). The constants κ of Assumption 4 and ρ of Proposition 1 are constants of the model and are absorbed into the constants of the proof. Throughout, ΠB(0,R) does not move a point of B(0, R), and we use the following fact: if v ∈ Rs has independent, mean-zero coordinates with sub-Gaussian norm at most ν, then p
∥v∥2 ≤ O ν s + log n with probability at least 1 − n−3 , (32) P which follows from Bernstein’s inequality applied to j vj2 . We write S = Ŝ and s = |Ŝ|, and set ⊤ θ̃i = 0 when Xi,S Xi,S is singular, so that the local step is defined on every dataset. We first prove privacy, which holds for every value of Ŝ and every dataset. The three groups of users, of sizes n1 , n3 and n4 , are disjoint, and a neighboring pair of user-level datasets changes the local data of one user, hence the messages of one group only. The randomizer of Stage 2a depends on the earlier rounds through Ŝ, and the randomizer of Stage 2b through (Ŝ, τ, c), which are computed from the messages of the other groups. So it is enough to show that each group, conditionally on the released quantities of the earlier groups, satisfies (ϵ, δ)-shuffled DP; the view of the data collector is then (ϵ, δ)-DP by parallel composition. Stage 1 was treated in the proof of Proposition 1. We give the argument for Stage 2b; Stage 2a is the same with (n3 , pa , R0 , ∆a , Ma , qa ) in place of (n4 , pb , τ, ∆b , Mb , qb ), and with θ̄i in place of ∆i . Condition on (Ŝ, τ, c) and fix two datasets that differ in the local data of one user i0 of Stage 2b. As in the proof of Theorem 2, write the rounding step as Ūi,j = ⌊pb ∆i,j + ui,j ⌋ with ui,j independent and uniform on (0, 1], and fix u. Since ∥∆i0 ∥2 ≤ τ and ∥∆′i0 ∥2 ≤ τ , √ √ ∥Ūi0 − Ūi′0 ∥2 ≤ pb ∥∆i0 − ∆′i0 ∥2 + s ≤ 2pb τ + s + 2 = ∆b . √ The setting of pb gives pb τ ≥ ( s + 2) max{1, ϵ}, so √ s max{1, ϵ} ≤ ∆b ≤ 3pb τ. (33) P P (b) Let Tj = i (Ūi,j + ζi,j ) over the n4 users of Stage 2b, so that i ζi,j = Bj − n4 N0 /2 with (b) (b) Bj ∼ Bin(n4 N0 , 1/2) independent across j, and n4 N0 ≥ Mb ≥ 32∆2b log(4s/δ)/ϵ2 . By equation 33 and ϵ ≤ log(4s/δ), Lemma 1 applies with this ∆b and with δ/2 in place of δ, so T = (T1 , . . . , Ts ) is (ϵ, δ/2)-DP for fixed u, and therefore also without conditioning on u. (b)

(b)

(b)

Since |Ūi,j | ≤ pb τ +1, |ζi,j | ≤ N0 /2 and N0 ≥ 2, we have |Tj | ≤ n4 (pb τ +1+N0 /2) < qb /2, so no wrap-around occurs and Tj is the value decoded by the aggregator. On coordinate j the users run the summation protocol Verb ,n4 of Balle et al. (2020) over Zqb ; the setting of rb gives rb ≥ 3, and n4 ≥ 19. By Balle et al. (2020, Corollary 6.1) the protocol has worst-case statistical security with parameter σsec = log2 (8seϵ /δ). Let V sim be the view obtained when, on every coordinate, the inputs are replaced by (Tj mod qb , 0, . . . , 0); it is a post-processing of T , hence (ϵ, δ/2)-DP. Coupling the s coordinates one at a time gives dT V (VP , V sim ) ≤ e−ϵ δ/8 for every dataset, so by Balle et al. (2020, Lemma F.2) the view of Stage 2b is (ϵ, δ/2+(1+eϵ )e−ϵ δ/8)-DP, hence (ϵ, δ)-DP. Therefore Algorithm 1 satisfies user-level (ϵ, δ)-shuffled DP. We now bound the error. By Proposition 1, with probability at least 1 − O(n−2 ) we have Ŝ = S ∗ and s = k, and we work on this event. In particular θŜ∗ = θS∗ ∗ carries all of θ∗ , so ∥θ̂ − θ∗ ∥2 = ∥θ̂Ŝ − θŜ∗ ∥2 . Consider the local least squares estimates. The rows of Xi,S ∗ are sub-Gaussian with covariance ΣS ∗ , and λmin (ΣS ∗ ) ≥ κ by Assumption 4. By the concentration of the smallest singular value of a ⊤ −cm sub-Gaussian design, for a fixed user λmin (Xi,S ∗ Xi,S ∗ ) ≥ κm/2 with probability at least 1−2e once m ≥ Ω(k), so by the third term of equation 10 and a union bound this holds for all users of Stage 2 with probability at least 1 − O(n−2 ). On this event, for every user of Stage 2, ⊤ −1 ⊤ θ̃i − θS∗ ∗ = (Xi,S Xi,S ∗ ζi , ∗ Xi,S ∗ )

30

Preprint

which, conditionally on Xi , has mean zero and is sub-Gaussian with covariance proxy 2 ⊤ −1 σ 2 (Xi,S ⪯ 2σ ∗ Xi,S ∗ ) κm Ik by Assumption 1. By equation 32 and a union bound over the users of Stage 2, with probability at least 1 − O(n−2 ),  r k + log n  ∗ ≤ τβ for every user i of Stage 2, (34) θ̃i − θS ∗ 2 ≤ O σ m where we used k + log n ≤ 2k log n and took Cβ large enough. We work on this event as well. We check that the two clipping steps do not change the local statistics. For Stage 2a, ∥θ̃i ∥2 ≤ ∥θ∗ ∥2 + τβ = R0 by equation 34, so θ̄i = θ̃i . The aggregator of Stage 2a decodes 1 X c= θ̃i + Eard + Eabin , n3 i where the sum is over the n3 users of Stage 2a, the coordinates of Eard are averages of n3 independent (a) terms in (−1, 1] divided by pa , and the coordinates of Eabin are centered Bin(n3 N0 , 1/2) q variables (a)

divided by n3 pa . Applying equation 32 to each of the three terms, together with n3 N0 ≤ √ p √ √ √ R0 log(1/δ) 1 √ O( Ma + n3 ), Ma ≤ O(pa R0 log(4s/δ)/ϵ) from equation 33 and n3 pa ≤ c n √log n ϵ p

3

from the setting of pa , we get, with probability at least 1 − O(n−2 ),  r k + log n R p(k + log n) log(4k/δ)  0 ∗ c − θS ∗ 2 ≤ O σ + ≤ Â, n3 m n3 ϵ using k + log n ≤ 2k log n and taking cA , c′A large enough. The first term of  is at most τβ /2 once Cβ ≥ 2cA . Since n3 = Θ(n), R0 = Θ(1) and log(4k/δ) = Θ(log(k/δ)), the second term of  is at most τβ /2 under equation 13. So  ≤ τβ , hence τ = τβ +  ≤ 2τβ .

(35)

For Stage 2b, ∥θ̃i − c∥2 ≤ ∥θ̃i − θS∗ ∗ ∥2 + ∥θS∗ ∗ − c∥2 ≤ τβ + Â = τ , so ∆i = θ̃i − c and the clipping does not change it either. The output of Stage 2b is therefore θ̂Ŝ − θS∗ ∗ =

1 X θ̃i − θS∗ ∗ + Ebrd + Ebbin , n4 i

where the sum is over the n4 users of Stage 2b. Conditionally on the designs, the summands are independent, have mean zero and are sub-Gaussian with proxy 2σ 2 /(κm), so their average has proxy 2σ 2 /(κn4 m), and equation 32 gives  r k + log n 
1 X θ̃i − θS∗ ∗ ≤O σ n4 i n4 m 2 with probability at least 1 − O(n−2 ). The two error terms are bounded as in Stage 2a, with τ in place of R0 and n4 in place of n3 :  √k + log n   τ p(k + log n) log(4k/δ)  1 bin , Ebrd 2 ≤ O ≤ E bin , Eb 2 ≤ O √ n4 ϵ n4 p b cp b 2 where pthe last step uses the setting of pb and log(1/δ) ≤ log(4k/δ). By equation 35 and τβ = Cβ σ k log n/m,  σ √k log n p(k + log n) log(4k/δ)   σk log nplog(4k/δ)  bin √ √ Eb 2 ≤ O ≤O , n4 m ϵ n4 m ϵ √ √ √ √ using k k + log n ≤ k + k log n ≤ 2k log n. Since n1 ≤ n/2, we have n3 , n4 = Θ(n). Collecting the three terms,  r k log n σk log nplog(4k/δ)  r k k  ∗ e √ θ̂ − θ 2 ≤ O σ + ≤O + √ , nm nm n m ϵ n mϵ since log(4k/δ) ≤ O(log n + log(1/δ)) for k ≤ n. All the events above hold together with probability at least 1 − O(n−2 ). 31

Preprint

E.3

P ROOF OF P ROPOSITION 2

Proof of Proposition 2. Throughout, the expectation is over the data and the randomness of A, which is independent of the data. For one sample z = (x, y), the density is fθ (y | x)fx (x), and fx does not depend on θ, so
y − µ(x⊤ θ) x Sθ (z) = ∇θ log fθ (y | x) = , µ = ψ′ . c(σ) We use two standard properties of the score. First, Z Z ∇θ fθ (z) E[Sθ (z)] = fθ (z) dz = ∇θ fθ (z) dz = 0. fθ (z)

(36)

Second, Var(y | x) = c(σ)ψ ′′ (x⊤ θ) for the exponential family, so   
E Var(y | x) xx⊤  1 = E ψ ′′ (x⊤ θ) xx⊤ . Cov Sθ (z) = 2 c(σ) c(σ) Pm Write Sθ (Ui ) = j=1 Sθ (zi,j ) for the score of the local dataset of user i. Since the m samples of Ui are i.i.d., E[Sθ (Ui )] = 0 and 
m  ′′ ⊤ E ψ (x θ) xx⊤ . (37) Cov Sθ (Ui ) = c(σ) ′ ′ We first prove soundness. The dataset Z(i) does not contain Ui , so A(Z(i) ) is independent of Ui , and by equation 36, ′ E[A′i ] = E[A(Z(i) )] − θ, E[Sθ (Ui )] = 0. ′ For the first absolute moment, we condition on A(Z(i) ) and use equation 37: r h q


i ′ ) − θ ⊤ Cov S (U ) A(Z ′ ) − θ E|A′i | ≤ E[(A′i )2 ] = E A(Z(i) θ i (i) r q
m ′ ) − θ∥2 · λmax E[ψ ′′ (x⊤ θ)xx⊤ ] . ≤ E∥A(Z(i) 2 c(σ)

ei are i.i.d., Z ′ and Z have the same distribution, so E∥A(Z ′ ) − θ∥2 = E∥A(Z) − Since Ui and U 2 (i) (i) 2 θ∥2 , which gives the bound in (i). Pn We next prove completeness. The samples of all users are i.i.d., so i=1 Sθ (Ui ) = ∇θ log fθ (Z), where fθ (Z) is the joint density of the nm samples. By equation 36 applied to Z, the term θ does not contribute: n X E[Ai ] = E A(Z) − θ, ∇θ log fθ (Z) = E A(Z), ∇θ log fθ (Z) . i=1

Writing the inner product over coordinates and exchanging differentiation and integration under the regularity conditions, d Z d i X h X   ∂ ∂ log fθ (Z) = E A(z)j fθ (z) dz E A(Z)j ∂θj ∂θj j=1 j=1 =

d X  ∂  E A(Z)j , ∂θj j=1

where the inner expectation is over the randomness of A. This proves (ii). The proof of Theorem 4 uses a sparse version of the user-level score attack, in which the score is kept only on the support of θ. For T = supp(θ), define

′ ATi = A(Z) − θ T , Sθ (Ui )T , ATi ′ = A(Z(i) ) − θ T , Sθ (Ui )T . (38) 32

Preprint

The proof above applies to ATi with every vector restricted to the coordinates in T . For soundness, E[ATi ′ ] = 0 and r q
m E|ATi ′ | ≤ E∥(A(Z) − θ)T ∥22 · λmax E[ψ ′′ (x⊤ θ)xT x⊤ T] . c(σ) For completeness, viewing fθ as a function of θT with θT c = 0 held fixed, n X

E[ATi ] =

i=1

E.4

X ∂   E A(Z)j . ∂θj

j∈T

P ROOF OF T HEOREM 4

Throughout this subsection we work on the parameter space  Θ = θ ∈ Rd : ∥θ∥0 ≤ k, ∥θ∥∞ ≤ 1 ⊆ Θk , we write T = supp(θ), and ATi , ATi ′ are the sparse user-level score attack of equation 38. The first lemma turns the differential privacy of the estimator into an upper bound on the score attack of one user; it is the user-level form of Proposition 2.1 of Cai et al. (2023). Lemma 2. Let A satisfy user-level (ϵ, δ)-DP with 0 < ϵ < 1, and let θ ∈ Θ. Then for every t0 > 0 and every i ∈ [n], Z ∞   T′ T P |ATi | > t dt. E[Ai ] ≤ 2ϵ E Ai + 2δt0 + t0

Proof of Lemma 2. Write A = ATi and A′ = ATi ′ , and let Z + = max(Z, 0) and Z − = ′ − min(Z, 0). The datasets Z and Z(i) are user-level neighbors, and A and A′ are the same ei ) the privacy of A gives function of (A(Z), Ui ) and of (A(Z ′ ), Ui ), so conditionally on (Ui , U (i)

P[A+ >R t] ≤ eϵ P[A′+ > t]R + δ for every t, and the same with the two sides exchanged. Using ∞ ∞ E[A] = 0 P[A+ > t] dt − 0 P[A− > t] dt and eϵ ≤ 1 + 2ϵ for 0 < ϵ < 1, Z ∞ Z t0 Z ∞
P[A+ > t] dt ≤ eϵ P[A′+ > t] + δ dt + P[A+ > t] dt 0 0 t Z ∞ Z ∞ 0 Z ∞ ′+ ≤ P[A > t] dt + 2ϵ P[A′+ > t] dt + δt0 + P[|A| > t] dt. 0

0

−

−ϵ

t0

′−

−ϵ

In the same way, from P[A > t] ≥ e (P[A > t] − δ) and e ≥ 1 − 2ϵ, Z ∞ Z ∞ Z ∞ − ′− P[A > t] dt ≥ P[A > t] dt − 2ϵ P[A′− > t] dt − δt0 . 0

0

0

Subtracting the second display from the first gives E[A] ≤ E[A′ ] + 2ϵ E|A′ | + 2δt0 + t] dt, and E[A′ ] = 0 by the soundness part of Proposition 2.

R∞ t0

P[|A| >

The next lemma sums Lemma 2 over the n users. The tail term is where the user-level regime differs from the sample-level one: the score of a user is a sum of m terms, so its size is governed by the empirical second moments of that user’s m design vectors. Lemma 3. Let Assumptions 2 and 4 hold, let ∥ψ ′′ ∥∞ ≤ c2 , and let θ ∈ Θ. Let A satisfy user-level (ϵ, δ)-DP with 0 < ϵ < 1 and take values in Θ, and let L = log(2k/δ) ≤ m. Then s r n q X  T mCc2 c2 √ 2 E Ai ≤ 2nϵ E∥A(Z) − θ∥2 + C1 nδkκx mL, c(σ) c(σ) i=1 where C1 is an absolute constant. 33

Preprint

Proof of Lemma 3. By the soundness part of Proposition 2 applied to ATi ′ , together with ∥ψ ′′ ∥∞ ≤ c2 and λmax (ΣT,T ) ≤ C from Assumption 4, s q mCc2 T′ 2 E Ai ≤ E∥A(Z) − θ∥2 . (39) c(σ) It remains to choose t0 in Lemma 2 and to bound the tail. Write wi,j = yi,j − ψ ′ (x⊤ i,j θ) and Pm Ui,ℓ = j=1 wi,j xi,j,ℓ for ℓ ∈ T , so that ATi =

1 c(σ)

ATi ≤

A(Z) − θ T , Ui,T ,

(A(Z) − θ)T c(σ)

1

max Ui,ℓ ≤ ℓ∈T

2k max Ui,ℓ , c(σ) ℓ∈T

where the last step uses |T | ≤ k and ∥A(Z)∥∞ , ∥θ∥∞ ≤ 1. For the exponential family, E[eλw/c(σ) | 2 x] ≤ ec2 λ /(2c(σ)) , so conditionally on x the variable w has mean zero and is sub-Gaussian with variance proxy c(σ)c2 . By Assumption 2 and Σℓℓ ≤ 1 we have ∥xℓ ∥ψ2 p ≤ κx , so each product wi,j xi,j,ℓ has mean zero and sub-exponential norm at most K = O(κx c(σ)c2 ). Bernstein’s inequality for the sum of the m independent products, followed by a union bound over ℓ ∈ T , gives h i  n u2 u o P max |Ui,ℓ | > u ≤ 2k exp − c min , . ℓ∈T mK 2 K √ √ Set a = 2kK/c(σ) and t0 = a( mL + L)/ c, so that P[|ATi | > t0 ] ≤ 2ke−L = δ and, integrating the same bound, Z ∞ p √

  P |ATi | > t dt ≤ O a m/L + a δ ≤ O a mL δ. t0

√ √ √ R∞ Since L ≤ m, we have mL + L ≤ 2 mL, so 2δt0 + t0 P[|ATi | > t]dt ≤ O(aδ mL). p p √
Substituting a = 2kK/c(σ) and K = O(κx c(σ)c2 ) bounds this by O δkκx c2 /c(σ) mL . Combining with equation 39 in Lemma 2 and summing over the n users gives the claim. The last lemma bounds the same sum from below, on average over a prior supported on sparse vectors. The prior is the one of Cai et al. (2023): let θe1 , . . . , θed be i.i.d. from the N (0, γ 2 ) distribution truncated to [−1, 1] with 1 1 γ2 = ≍ , (40) 4 log(d/4k) log(d/k) let S be the set of indices of the k largest |θej |, and set θj = θej I[j ∈ S]. We write π for the law of θ, and pSl for its density conditioned on S = Sl , where S1 , . . . , S(d) enumerate the subsets of [d] of k size k. √ Lemma 4. Let k ≤ d and k ≥ C0 log2 (d/k) for a large enough constant C0 . Let A take values in Θ and satisfy EZ|θ ∥A(Z) − θ∥22 ≤ R for every θ ∈ Θ, where √ c1 k R≤ (41) log(d/k) for a small enough constant c1 . Then n X

  Eπ EZ|θ ATi ≥ c k log(d/k)

i=1

for an absolute constant c > 0. Proof of Lemma 4. Let g(θ) = EZ|θ [A(Z)]. By the completeness part of equation 38, for every θ Pn P with supp(θ) = T , i=1 EZ|θ [ATi ] = j∈T ∂gj (θ)/∂θj . The prior π is not absolutely continuous
−1 P on Rd , so we condition on the support: π = kd l pSl , and n X

X ∂gj (θ)   1 X Eπ EZ|θ ATi = d
EpSl . ∂θj k i=1 l

34

j∈Sl

Preprint

Within a fixed Sl the law of θ has a density, so Stein’s lemma applies coordinatewise and gives EpSl

X ∂gj (θ) j∈Sl

∂θj

≥ EpSl

hX j∈Sl

−θj

hX p′Sl ,j (θj ) i p′Sl ,j (θj ) i + B, − E pS l gj (θ) − θj pSl ,j (θj ) pSl ,j (θj ) j∈Sl

where B collects the boundary terms at ±1. Each θj is an order statistic of dptruncated nor2 mals, so its density vanishes at −1, and p √ at +1 it is bounded by dϕ(1) = O(k log(d/k)/d), which is O( d. As A takes values in Θ, |gj (θ) − θj | ≤ 2, so p log(d/k)) since k ≤ |B| = O(k log(d/k)) ≪ k log(d/k). The two remaining terms depend on the prior only. By Lemma C.1 of Cai et al. (2023), the first term √ is at least c′ k log(d/k) for an absolute constant c′ , using equation 40 and k ≤ d. For the second term, the Cauchy–Schwarz inequality and the bound on the score of the order statistics in the proof of Lemma C.2 of Cai et al. (2023) give v u hX X  p′S ,j (θj ) 2 √
p′Sl ,j (θj ) i √ u l ≤ R tEpSl EpSl gj (θ)−θj ≤ R O k 3/4 log3/2 (d/k) , pSl ,j (θj ) pSl ,j (θj ) j∈Sl

j∈Sl

√ where we used EZ|θ ∥A(Z) − θ∥22 ≤ R and j≤k j log3 (d/k) = O(k 3/2 log3 (d/k)). By equa√ tion 41 this is O( c1 k log(d/k)), which is at most c′ k log(d/k)/4 for c1 small enough. The same bounds hold for every Sl , so averaging over l gives the claim with c = c′ /2. P

Proof of Theorem 4. It is enough to bound inf A supθ∈Θ E∥A(Z) − θ∥22 from below, since Θ ⊆ Θk . We may also assume that A takes values in Θ: replacing A(Z) by a closest point of Θ is a post-processing, hence preserves privacy, and at most doubles the error at every θ ∈ Θ. Let R = supθ∈Θ E∥A(Z) − θ∥22 and let L = log(2k/δ). We first p 3. The map
pcheck that the conditions of the theorem control the tail term of Lemma x 7→ x log(2k/x) is increasing on (0, 2k/e), so the condition δ = o log(d/k)/(n m log(nm)) together with k ≤ nm gives s √ c(σ) 1 ′ nδ mL ≤ c0 log(d/k) (42) c2 κx for a constant c′0 as small as we wish, and m ≥ log(nm) gives L ≤ m once δ ≥ 2k/(nm)C ; for smaller δ the left-hand side of equation 42 only decreases, and L ≤ m can be assumed by replacing δ by the larger value. Hence the tail term of Lemma 3 is at most 2c k log(d/k) with c the constant of Lemma 4. √
√ Suppose first that R > c1 k/ log(d/k). The condition n m ϵ ≥ Ω k 3/4 log3/2 (d/k) gives  √k  k 2 log2 (d/k) ≤O , n2 mϵ2 log(d/k) so R is already of√the order claimed by the theorem and there is nothing to prove. We may therefore assume R ≤ c1 k/q log(d/k), so that Lemma 4 applies. Combining it with Lemma 3, averaged √ over π and using Eπ EZ|θ ∥A(Z) − θ∥22 ≤ R, √ c k log(d/k) ≤ 2nϵ R and therefore R≥Ω

s

mCc2 c + k log(d/k), c(σ) 2

 c(σ) k 2 log2 (d/k)  · . Cc2 n2 mϵ2

The privacy constraint is not the only one: any estimator, private or
not, is bounded below by the non-private minimax risk over Θ, which is Ω c(σ)k log(d/k)/(nm) for the N = nm samples; for 35

Preprint

the linear model ψ(u) = u2 /2 this is the rate of Raskutti et al. (2011), whose design conditions are supplied by Assumption 4. Adding the two bounds and absorbing C, c2 and c(σ) into the constants,  k log(d/k) k 2 log2 (d/k)  2 inf sup E A(Z) − θ∗ 2 ≥ Ω , + A∈Mϵ,δ θ ∗ ∈Θk nm n2 mϵ2 and taking square roots gives the statement of Theorem 4.

F

A DDITIONAL R ESULTS

F.1

N ON - INTERACTIVE U SER - LEVEL S HUFFLED DP

Theorem 3 removes the dimension from the rate by using two rounds. We show here that the second ingredient, interaction, is also needed: under a non-interactive protocol, user-level data alone leaves a polynomial dependence on d in the privacy term, whatever the local sample size m is. Theorem 5. Let MSDP ϵ,δ be the family of all non-interactive protocols that satisfy robust user-level (ϵ, δ)-shuffled DP with 0 < ϵ ≤ O(1) and δ = o(1/n). Let Z = {Ui }ni=1 be a dataset in which each user holds m i.i.d. samples from the sparse linear model of Assumption 1 with parameter θ∗ . If n3/2 ϵ2 ≥ Ω(dk), then r inf

  sup E ∥A(Z) − θ∗ ∥2 ≥ Ω

∗ A∈MSDP ϵ,δ θ ∈Θk

(43) ! √ k log(d/k) dk + √ . nm n mϵ

√ Comparing with Theorem √ 3, the privacy term of the non-interactive protocol carries d, while the interactive one carriesp k in its place. The first term is the non-private minimax rate for N = nm samples, which is Ω( k log(d/k)/(nm)) (Raskutti et al., 2011); we only prove the privacy term below. Proof of Theorem 5. We use the notations of Appendix D.1: X = {±1}d+1 , U is uniform on X , χi (x, y) = yxi , and Qθ = U (1 + y⟨x, θ⟩). For z ∈ {−1, 0, 1}d and a scale β > 0 we write Dz(β) = Q⊗m βz , the law of the m samples of one user at parameter βz, and Z0 (z), z i,b are as in Appendix D.1. All densities below are taken with respect to U ⊗m , and ⟨·, ·⟩ and ∥ · ∥ are the inner product and the norm of L2 (U ⊗m ). Fix a scale γ with mkγ 2 ≤ 1, (44) ′ and let γ = 2γ/9. Let P be a protocol in MSDP , which we take to be (ϵ, δ, 1/3)-robust, and let ϵ,δ n′ = n/3. By Cheu & Ullman (2021, Theorem 2.6) applied with the record universe X m and the fixed distribu(0) tion U ⊗m = Dz , there is an (ϵ, δ)-pan-private algorithm M that processes n′ records and satisfies, for every z, 
 ′
7 2 (45) dT V M (Dz(γ) )⊗n , MP Q⊗n ≤ e−c1 n , Qz = Dz(γ) + U ⊗m . z 9 9 The mixture Qz is not a member of the model, and this is the one point where the user-level regime differs from Appendix D.1: the law of a record is Q⊗m βz , which is not linear in the parameter. We (γ ′ )

therefore compare Qz with the member of the model that agrees with it to first order, namely Dz Writing d X X Y zl χl , ej = χz (wt ), χz = T ⊆[m], |T |=j t∈T

l=1

Pm P (β) we have Dz = j=0 β j ej and Qz = 1 + 29 j≥1 γ j ej , so the two first-order terms agree and ′

Rz := Qz − Dz(γ ) =

X 2 j≥2

36

9

−

 2 j  9

γ j ej .

.

Preprint

The characters appearing in ej for different j are distinct and ∥ej ∥2 = √ with u = γ mk ≤ 1 from equation 44, X uj X 1 √ ≤ u2 √ ≤ 3mkγ 2 . ∥Rz ∥ ≤ j! j! j≥2 j≥0

j m j j ∥z∥0 ≤ (mk) /j!, so

Since Qz ≥ 97 pointwise, the squared Hellinger distance satisfies Z
2

9 Rz2 2 (γ ′ ) 2 (γ ′ ) Qz = dH Dz , Qz ≤ χ Dz dU ⊗m ≤ ∥Rz ∥2 ≤ 12 mkγ 2 , Qz 7 √ and by the subadditivity of d2H over product measures and dT V ≤ 2 dH , 
 √ (γ ′ ) ⊗n (46) dT V Q⊗n ≤ 6 n mkγ 2 . z , Dz √ Combining equation 45 and equation 46, and writing σn = e−c1 n + 6 n mkγ 2 , we get for every z 
 ′
′ dT V M (Dz(γ) )⊗n , MP (Dz(γ ) )⊗n ≤ σn . (47) So it is enough to bound the total variation distances between the outputs of the pan-private algorithm M , and to pay σn per comparison. We bound these distances as in Appendix D.1. Fix z with ∥z∥0 ≤ k−1 and let i ∈ Z0 (z), b ∈ {±1}. For j ∈ [n′ ], the state Sj of M is, for every fixed prefix, an (ϵ, δ)-DP function of the j-th record, so by Cheu & Ullman (2021, Lemma 3.7) there is a (2ϵ, 0)-DP kernel qj′ within total variation distance (γ)

(γ)

δ of it. Running the hybrid in which the first j records are drawn from Dz and the rest from Dzi,b , the output is a post-processing of (Sj , X>j ), the prefix has the same law in the two hybrids, and the suffix is independent of the j-th record, so
1 (γ) (48) Dib = Dzi,b − Dz(γ) , dT V Hj−1 , Hj ≤ EX<j Es∼q̄j′ hj,s , Dib + 2δ, 2 where q̄j′ (s) = EU ⊗m [qj′ (s | W )] and hj,s = qj′ (s | ·)/q̄j′ (s) − 1 satisfies ∥hj,s ∥∞ ≤ e2ϵ − 1. Neither the law of the prefix nor the kernel depends on i. Two factsQ about Dib areQ used. First, expanding the two products, every character appearing in Dib is of the form t∈T χi (wt ) t∈T ′ χlt (wt ) with T ̸= ∅, T ∩T ′ = ∅ and lt ∈ supp(z). Since i ∈ / supp(z), the characters of Dib and Dib′ are distinct for i ̸= i′ in Z0 (z). Second, writing k ′ = ∥z∥0 ,
m
m 2 ′ 2 Dib = 1 + γ 2 (k ′ + 1) − 1 + γ 2 k ′ ≤ eγ m(k +1) mγ 2 ≤ e2 mγ 2 , where the last step uses equation 44. By the first fact and the Cauchy–Schwarz inequality, for every h with ∥h∥∞ ≤ e2ϵ − 1, X √ √ √ s X
⟨h, Dib ⟩2 ≤ d ∥h∥ max Dib ≤ e dm γ e2ϵ − 1 . ⟨h, Dib ⟩ ≤ d i

i∈Z0 (z)

i∈Z0 (z)

Summing equation 48 over j ∈ [n′ ], over b ∈ {±1} and over i ∈ Z0 (z), and then adding 4σn per coordinate from equation 47, X √
Ti (z) ≤ e γ n′ dm e2ϵ − 1 + 4dn′ δ + 4dσn , (49) i∈Z0 (z)

where Ti (z) =

P

b dT V

(γ ′ ) ⊗n

(γ ′ )

MP ((Dzi,b )⊗n ), MP ((Dz

)

) is stated for the protocol.

The rest is Assouad’s method, exactly as in Appendix D.1. Draw Z as in Definition 6, let θZ = γ ′ Z, (γ ′ ) draw the n local datasets from DZ given Z, and let θ̂ = MP (Z). Let Gi = {∥Z−i ∥0 ≤ k − 1}, P let ψi = sign(θ̂i ) and let N = i I[Zi ̸= 0, Gi , ψi ̸= Zi ]. Then N ≤ k whenever N ≥ 1, and √ √ P ∥θ̂ − θZ ∥2 ≥ γ ′ N ≥ γ ′ N/ k. Le Cam’s method, P(Zi ̸= 0) = k/(2d), i P(Gi ) ≥ d/2 by Markov’s inequality, and equation 49 together give  √
k d E[N ] ≥ − 2eγn′ dm e2ϵ − 1 − 8dn′ δ − 8dσn . 4d 2 37

Preprint

√ 
√ We now fix γ. Take γ = d C2 n′ m(e2ϵ − 1) with C2 a large enough constant, so that the second term is at most d/16; the third term is at most d/16 because δ = o(1/n); and the fourth is at most d/16 once σn ≤ 1/128, which holds for n larger than an absolute constant together with √

n mkγ 2 =

C22 n′2

√

1 dk ≤ , 768 n (e2ϵ − 1)2

which is equation the small signal condition equation 44 √ 43 for ϵ ≤ O(1). For this choice of2 γ, 2 reads nϵ ≥ Ω( dk), which equation 43 implies since n ϵ ≥ n3/2 ϵ2 . Hence E[N ] ≥ k/32, and ! √ √ ′     k dk γ √ sup E ∥θ̂ − θ∗ ∥2 ≥ E ∥θ̂ − θZ ∥2 ≥ =Ω , 32 n mϵ θ ∗ ∈Θk using γ ′ = 2γ/9, n′ = n/3 and e2ϵ − 1 ≤ O(ϵ). F.2

U SER - LEVEL DP IN D ENSE R EGRESSION

The user-level score attack is not tied to sparsity. We record here what it gives for the dense problem, where the parameter is only assumed to lie in the box Θbox = {θ ∈ Rd : ∥θ∥∞ ≤ 1}. The rate is the one of Theorem 4 with k log(d/k) replaced by d. Theorem 6. Let Z = {Ui }ni=1 be a dataset in which each user holds m i.i.d. samples from the exponential family of Theorem 4. Let Assumption 2 hold with λmax (Σ) ≤ C, and let ∥ψ ′′ ∥∞ ≤ c2 . Let Mϵ,δ of all estimators satisfying user-level p be the family √ (ϵ, δ)-DP with 0 < ϵ < 1 and δ =
√ o 1/(n m log(nm)) . If m ≥ log(nm) and n m ϵ ≥ Ω( d), then ! r   d d ∗ . + √ inf sup E ∥A(Z) − θ ∥2 ≥ Ω A∈Mϵ,δ θ ∗ ∈Θbox nm n m ϵ Throughout this subsection Ai is the user-level score attack of Definition 7, which is equation 38 P with T = [d]. The first lemma bounds i E[Ai ] from below through Stein’s lemma, under a prior that is supported on Θbox and vanishes at its boundary. 2 2 Lemma 5. Let π be the product prior on Θbox with marginal densities πj (t) = 15 16 (1 − t ) on 2 [−1, 1]. Let A take values in Θbox and satisfy EZ|θ ∥A(Z) − θ∥2 ≤ R for every θ ∈ Θbox , with R ≤ d/40. Then n X   d Eπ EZ|θ Ai ≥ . 2 i=1

Proof of Lemma 5. Let g(θ) = EZ|θ [A(Z)]. By the completeness part of Proposition 2, Pd Pn i=1 EZ|θ [Ai ] = j=1 ∂gj (θ)/∂θj . Each πj is continuously differentiable on (−1, 1) and vanishes at ±1, so Stein’s lemma applies coordinatewise with no boundary term, and, as in Proposition 2.2 of Cai et al. (2023), v u d d d  ′ hX hX X πj (θj ) 2 i πj′ (θj ) i u ∂gj (θ) Eπ ≥ Eπ −θj − tR Eπ . ∂θj πj (θj ) πj (θj ) j=1 j=1 j=1 For this prior, πj′ (t)/πj (t) = −4t/(1 − t2 ), so Z Z 1 h h π ′ (t) 2 i πj′ (t) i 15 1 2 j 2 Eπj − t = t (1 − t ) dt = 1, Eπj = 15 t2 dt = 10. πj (t) 4 −1 πj (t) −1 √ The first term is therefore d and the second is 10Rd ≤ d/2 for R ≤ d/40. The second lemma bounds the same sum from above through privacy. It is the dense form of Lemma 3, with ∥(A(Z) − θ)∥1 ≤ 2d in place of 2k and the union bound taken over all d coordinates. 38

Preprint

Lemma 6. Under the assumptions of Theorem 6, let A take values in Θbox and satisfy user-level (ϵ, δ)-DP with 0 < ϵ < 1, and let L = log(2d/δ) ≤ m. Then s r n q X   mCc2 c2 √ E Ai ≤ 2nϵ E∥A(Z) − θ∥22 + C1 nδdκx mL, c(σ) c(σ) i=1 where C1 is an absolute constant. Proof of Lemma 6. The proof of Lemma 3 applies p verbatim with T p = [d]. The soundness part ′ 2 of Proposition 2 and λmax (Σ) ≤ C give E|Ai | ≤ E∥A(Z) − θ∥2 mCc2 /c(σ). For the tail, ∥(A(Z) − θ)∥1 ≤ 2d because both vectors lie in Θbox , so |Ai | ≤ 2d maxℓ∈[d] |Ui,ℓ |/c(σ) with Pm Ui,ℓ = j=1 (yi,j − ψ ′ (x⊤ i,j θ))xi,j,ℓ . Bernstein’s inequality and a union bound over ℓ ∈ [d] give the p √ same tail bound as there with k replaced by d, and the choice t0 ≍ dκx c2 /c(σ) mL in Lemma 2 gives the claim. Proof of Theorem 6. We may assume that A takes values in Θbox , since projecting onto this set is a post-processing and does not increase the error at any θ ∈ Θbox . Let R = supθ∈Θbox E∥A(Z) − θ∥22 and L = log(2d/δ). √ √ Suppose first that R > d/40. The condition n m ϵ ≥ Ω( d) gives d2 /(n2 mϵ2 ) ≤ O(d), so R already exceeds the bound claimed by the theorem and there is nothing to prove. We may therefore assume R ≤ d/40, so that Lemma 5 applies. As in the proof of Theorem 4, the condition on δ together with m ≥ log(nm) gives L ≤ m and r d c2 √ C1 nδdκx mL ≤ . c(σ) 4 Combining Lemma 5 with Lemma 6, averaged over π, s   √ c(σ) mCc2 d d2 d ≤ 2nϵ R + , hence R≥Ω · 2 2 . 2 c(σ) 4 Cc2 n mϵ
Adding the non-private minimax rate Ω c(σ)d/(nm) for N = nm samples, absorbing C, c2 and c(σ) into the constants and taking square roots gives the statement.

39

Related documents

• On the Limits of Sparse Linear Model with Shuffled Differential Privacy#1448338
• Privacy-Preserving Fully Distributed Gaussian Process Regression#1448018
• High-Dimensional Asymptotics of Differentially Private PCA#1187270
• Exact Privacy--Distortion Design of Support-Constrained Randomized Response#1327562
• Minimax and Adaptive Covariance Matrix Estimation under Differential Privacy#1104987