Private Stochastic Convex Optimization and Sparse Learning with Heavy-tailed Data Revisited

1 Introduction

As one of the most fundamental problems in machine learning and statistics, Stochastic Convex Optimization (SCO) [Vapnik, 1999] with its empirical form, Empirical Risk Minimizataion (ERM), has been widely studied. Both SCO and ERM have found numerous applications in many areas such as medicine, finance, genomics and social science. However, due to the widespread concerns on privacy, how to handle sensitive data, such as biomedical datasets, has become a big hurdle for successful implementations of SCO in practice. To address the privacy issue, Differential Privacy (DP) [Dwork et al., 2006] has established itself as a canonical privacy notation for privacy-preserving data analysis.

∗Contact Author

The study of SCO and ERM under DP constraint (i.e., DPSCO and DP-ERM) has received significant attentions over the past decade. A long list of works have studied the problem from different perspectives: [Bassily et al., 2014; Bassily et al., 2019; Feldman et al., 2020] studied the problems in the low dimensional case and the central DP model, [Cai et al., 2020] considered the problems in the high dimensional sparse case and the central DP model, [Duchi et al., 2018] focused on the problems in the local DP model.

Even though there are numerous works on DP-SCO, a critical issue in most existing results is that loss function has to be assumed to satisfy the O(1)-Lipschitz property, or the underlying data distribution is assumed to be sub-Gaussian or even bounded. Despite simplifying the procedure of designing DP algorithms, such assumptions are unrealistic and may not always hold when dealing with real-world datasets, especially those from biomedicine and finance, as it has been observed that they are often heavy-tailed [Woolson and Clarke, 2011]. The heavy-tailed data could lead to unbounded gradients and thus break the Lipschitz assumption, which implies that previous algorithms may fail to provide DP guarantee. Recently, to tackle this issue, there have been several works studying DP-SCO with heavy-tailed data [Wang et al., 2020b; Kamath et al., 2021; Hu et al., 2021] or private mean estimation for heavy-tailed distributions [Barber and Duchi, 2014; Kamath et al., 2020; Liu et al., 2021]. However, all these results still need to assume that the distribution of each coordinate of the gradient of the loss function has bounded second moment, which implies the data is still well-behaved to some extent. Thus, a natural question is,

Can we further relax the bounded second moment condition to model the data distribution that are more heavytailed? And what are the theoretical behaviors of DP-SCO with more extremely heavy-tailed data?

In this paper, we revisit the problem of DP-SCO with heavy-tailed data under more relaxed assumptions. For the first time, we consider the case with more extremely heavytailed data such that the distribution of each coordinate of the loss gradient has only bounded (1 + v)-th moment for some v ∈ (0, 1]. Our contributions can be summarized as follows.

• First, we consider one-dimensional private mean estimation for heavy-tailed distributions. We propose a novel robust and (ϵ, δ)-DP estimator based on truncating

the data, which achieves an error of O(( √log 1 δ nϵ ) v 1+v ), where n is the number of data samples. We then show that our proposed estimator is optimal by providing the matching lower bound on the estimation error.

• Based on the idea in the one-dimensional case, we then extend to estimate the mean of heavy-tailed distribution in a general d-dimensional space and use it to DP-SCO. Specifically, we consider both strongly convex and general convex loss functions for DP-SCO, and propose (ϵ, δ)-DP algorithms that achieve an error of ˜O( d 1+4v 1+v (ϵn) 2v 1+v ) and ˜O( d 1+4v 2+2v (ϵn) v 1+v + d 3+12v 2+2v (ϵn) 3v 1+v ) respectively, if we omit other terms. We also provide the lower bounds under our setting in both (ϵ, δ)-DP and ϵ-DP, and show that our upper bound is optimal up to a factor of Poly(d).
• Finally, to mitigate the high dimensionality issue, we also study DP-SCO with heavy-tailed data under the sparsity constraint, i.e., DP sparse learning. Based on the previous ideas, we propose a new method under our setting and show that it is possible to achieve an error of ˜O((s∗ ) 1+2v 1+v ( log d nϵ ) 2v 1+v ), where s∗ is the underlying sparsity. We also proof that the bound is optimal up to a factor of O(s∗ ).

Due to space limit, all the proofs are included in Appendix.

2 RelatedWork

As mentioned in Section 1, there are a vast number of works on DP-SCO and DP-ERM. Due to space limit, here we just discuss the results that are the most related to ours. For DPSCO with heavy-tailed data, [Wang et al., 2020b] first studies the problem by proposing three methods based on different assumptions. However, all the three methods need to assume the distribution of gradient of the loss is sub-exponential or has at least bounded second moment. [Kamath et al., 2021] recently revisits the problem under the same assumption as in [Wang et al., 2020b] and improves the (expected) excess population risk for both convex and strongly convex loss functions. It also provides the lower bounds in both (ϵ, δ)-DP and ϵ-DP models. Note that although some ideas of our algorithms and proofs are the same as theirs, their methods cannot be used in our relaxed setting and there are several differences. See Remark 1 and 2 for details.

DP sparse learning has been studied previously [Wang and Gu, 2019; Wang and Xu, 2019; Wang and Xu, 2021]. However, all of the previous methods need either the loss function be Lipschitz, or the data distribution be sub-Gaussian or even bounded [Cai et al., 2019]. [Hu et al., 2021] recently extends to the heavy-tailed case where each coordinate of the gradient has bounded second moment. However, due to their private estimator, it is impossible to use their methods to the bounded (1 + v)-th moment case. See Remark 3 for details.

3 Preliminaries

Definition 1 (Differential Privacy (DP)[Dwork et al., 2006]). A randomized algorithm M : X n ↦→ Y satisfies (ϵ, δ)- differential privacy if for every pair of neighbouring datasets X,X′ ∈ X n (i.e., datasets that differ in exactly one entry), it holds for ∀Y ⊆ Y that

When δ = 0, we call M as ϵ-DP.

Lemma 1 (Adaptive Composition Theorem). Given target privacy parameters 0 < ϵ < 1 and 0 < δ < 1, to ensure (ϵ, δ)-DP over m mechanisms, it suffices that each mechanism is (ϵ′ , δ′ )-DP, where ϵ′ = ϵ 2√2mln(2/δ) and δ′ = δ 2m .

Lemma 2 (Laplacian Mechanism). Given a dataset D ∈ X n and a function q : X n → Rd , the Laplacian Mechanism is defined as q(D) + (Y1, Y2, · · · , Yd), where each Yi is i.i.d. sampled from the Laplacian Distribution Lap( ∆1(q) ϵ ), where ∆1(q) is the ℓ1-sensitivity of the function q, i.e., ∆1(q) = supD∼D′ ||q(D) - q(D′ )||1. The density of the Laplacian distribution with parameter λ is Lap(λ)(x) = 1 2λ exp(- xλ ). Laplacian mechanism preserves ϵ-DP.

Lemma 3 (Gaussian Mechanism). Given a dataset D ∈ X n and a function q : X n → Rd , the Gaussian mechanism is defined as q(D) + ξ where ξ ∼ N(0, 2∆22 (q) log(1.25/δ) ϵ2 Id), where ∆2(q) is the ℓ2-sensitivity of the function q, i.e., ∆2(q) = supD∼D′ ‖q(D) - q(D′ )‖2. Gaussian mechanism preserves (ϵ, δ)-DP.

Definition 2 (DP-SCO [Bassily et al., 2014]). Let D be some unknown distribution over the data universe X and X = {x1, . . . , xn} ⊆ X n be i.i.d samples from the distribution D. Given a convex constraint set W ⊆ Rd and a convex loss function ℓ : W×X ↦→ R. Differentially Private Stochastic Convex Optimization (DP-SCO) is to find wpriv so as to minimize the population risk, i.e., LD(w) = Ex∼D[ℓ(w, x)] with the guarantee of being differentially private. The utility of an algorithm A for DP-SCO is measured by the expected excess population risk, which is defined as follows:

4 One-dimensional Private Mean Estimation for Heavy-tailed Distributions

Before studying DP-SCO with heavy-tailed data, we start from the one-dimensional private mean estimation for heavytailed distributions to illustrate the idea of our following methods more clearly. Here we assume the data distribution has bounded (1 + v)-th moment with some v ∈ (0, 1]. Formally, we are given a dataset X = {x1, . . . , xn} with each xi ∈ R sampled from the one dimension distribution D such that Ex∼D[|x|1+v ] ≤ u = O(1). 1 We aim to privately estimate the mean of the distribution, i.e., µ = Ex∼D[x]. Note that here we use the raw moment, which also implies that its central moment is bounded, i.e., Ex∼D[|x-Ex∼D[x]|1+v ] ≤ O(1), and vice versa. See Lemma 13 in Appendix for details.

We propose a private and robust estimator based on truncation. Specifically, for each data sample xi, if its magnitude is within the designed threshold B then we will keep it, otherwise we set it be 0. After the preprocessing, each sample now Algorithm 1 Truncation Based DP Mean Estimator: DPODME Tϵ,δ,ξ(X)

1 Throughout the whole paper we assume that v and u are known.

Input: Data samples X = {xi}ni=1 , xi ∈ R. Parameters ϵ, δ, ξ

Output: A private mean estimator ˜ µ ∈ R

1: Truncate each data sample xi by x′i ← xi ·1|xi|≤B, where

• 2: Get the robust mean estimator ̂ µ ← 1n ∑ni=1 x′i .
• 3: return ˜ µ ← ̂ µ+ ν, where ν ∼N(0, 8B2 n2 ϵ2 log 1.25 δ ).

is bounded in [-B,B]. Thus we can add Gaussian noise to the mean of the truncated data. See Algorithm 1 for details.

Theorem 1. For any 0 < ϵ, δ ≤ 1, Algorithm 1 DPODME Tϵ,δ,ξ(X) satisfies (ϵ, δ)-DP. Moreover, given any failure probability ξ, with probability at least 1 - ξ, the output ˜ µ satisfies

Remark 1. When v = 1, the error becomes O( 1 √nϵ ) (if we omit other terms), which matches the bound in [Wang et al., 2020a] and is optimal [Kamath et al., 2020]. However, previous results are for the case where the distribution has bounded second moment. Here we relax it to the case where the distribution only has its (1 + v)-th moment bounded. Thus, our method can be thought as a generalization of the previous results. Recently [Tao et al., 2021] also considers a similar problem. However, there are several differences: First, [Tao et al., 2021] focuses on the online setting in the ϵ-DP model while we consider the offline setting and the (ϵ, δ)-DP model. Secondly, the methods in [Tao et al., 2021] are based on the tree mechanism and Laplacian mechanism, while we mainly use the Gaussian mechanism. Thus, our error bound is lower. Thirdly, besides the upper bound, we also show that the bound of O((u 1 v+1 ( 1 nϵ ) v 1+v )) in Theorem 1 is optimal by showing its lower bound in the Theorem 2 below, which has not been studied in [Tao et al., 2021].

Theorem 2. There exists a distribution D with mean µ and its (1+v)-th raw moment is bounded by u. For any (ϵ, δ)-DP algorithm, its output ˜ µ satisfies the following with at least a constant probability

When v = 1 and u = 1, our result will be equivalent to the lower bound in [Kamath et al., 2020]. Thus, Theorem 2 is a generalization of the previous result on the lower bound of private mean estimation for heavy-tailed distributions. Moreover, besides the (ϵ, δ)-DP model, the lower bound also holds for any ϵ-DP algorithm.

5 DP-SCOwith Heavy-Tailed Data

In this section, based on the previous one dimensional private mean estimator, we provide our methods for DP-SCO with heavy-tailed data. Before that, we provide the assumptions that will be used throughout the section.

Definition 3 (Lipschitz). A function f : W → R is LLipschitz if for all w1, w2 ∈ W we have |f(w1) - f(w2)| ≤ L‖w1 - w2‖2.

Definition 4 (Strong convexity). A function f is α-strongly convex on W if for all w1,w2 ∈ W we have f (w1) ≥ f (w2) + 〈∇f (w2) , w1 - w2〉 + α2 ‖w1 - w2‖22 .

Definition 5 (Smoothness). A function f is β-smooth on W if for all w1,w2 ∈ W we have f(w1) ≤ f(w2) +

Assumption 1. We make the following assumptions:

1. The parameter space W is convex and bounded with diameter ∆, i.e., for ∀w1,w2 ∈W, ‖w1 -w2‖2 ≤∆.
2. The loss function ℓ(w, x) is non-negative and differentiable for all w ∈ W and x ∈ D.
3. The population risk LD(·) is β-smooth over W. For any w ∈ W, the gradient of the population risk function satisfies ‖∇LD(w)‖2 ≤ R = O(1). Moreover, the optimal solution w∗ = argminw∈W LD(w) satisfies that ∇LD(w∗ ) = 0.
4. For any w ∈ W, the distribution of each coordinate of the gradient of the loss function has bounded (1 + v)- th (raw) moment with some v ∈ (0, 1], i.e., there is a constant u > 0 such that Ex∼D[|∇jℓ(w, x)|1+v ] ≤ u for all j ∈ [d].

There are several notes on the terms in Assumption 1. Firstly, the first three terms in Assumption 1 are commonly used in the previous works on DP-SCO with heavy-tailed data [Kamath et al., 2021; Wang et al., 2020b]. The fourth condition assumes that the gradient of the loss is heavy-tailed, which is a commonly used assumption in the study of robust learning with heavy-tailed data. However, as mentioned previously, most of those works only assume that the gradient has at least bounded second moment while here we relax to the (1+v)-th moment. Moreover, we can see it is a relaxation of the Lipschitz condition that ‖∇LD(w)‖2 ≤ L for all w.

In Algorithm 2, we propose our framework. The idea of the algorithm is quite straightforward: we use the private version of Projected Gradient Descent (PGD). Specifically, in each iteration t, we first privately estimate the vector ∇LD(wt-1) by using gradients {∇ℓ(wt-1, xi)}ni=1 , where wt-1 is the current parameter. Then we update the parameter via the PGD. In the classical setting where the data is regular or the loss function is Lipschitz, we can use the Gaussian mechanism to the average of the gradients ∑ni=1 ∇ℓ(wt-1,xi) n to get a private estimation of ∇LD(wt-1). However, due to the heavy-tailed assumption on gradients in Assumption 1, in our problem we cannot use the same approach directly as now the ℓ2-norm sensitivity maybe infinite. Thus the main difficulty is designing private estimator for ∇LD(w) = Ex∼D[∇ℓ(w, x)],

Algorithm 2 DP-SCOϵ,δ,η,T,τ

Input: Data samples X = {xi}ni=1 ⊆ Rd , parameters ϵ, δ, η, T, τ . Output: Private minimizer wpriv . 1: for t ← 1, · · · , T do 2: if ℓ is convex then 3: Set Xt =X, ϵ′ ← ϵ 2√2T log(2/δ) , δ′ ← δ 2T . 4: else if ℓ is strongly convex then 5: Set Xt = {x(t-1)n/T+1, · · · , xtn/T }. 6: Set ϵ′ ← ϵ, δ′ ← δ. 7: end if 8: ∇˜LD(wt-1) ←DPHDMEϵ′ ,δ′ ,τ({∇ℓ(wt-1, x)}x∈X) 9: wt ←ProjW(wt-1 - η∇˜LD(wt-1)) 10: end for 11: if LD(·) is strongly convex then 12: Set wpriv ←wT 13: else if LD(·) is convex then 14: Set wpriv ← 1 T ∑t∈[T ] wt 15: end if 16: return wpriv

which could be seen as an instance of the private mean estimation in the d-dimensional space.

In Section 4 we considered the case where d = 1. Now we will use its idea in the general high dimensional case. Our estimator is presented in Algorithm 3. For each coordinate, we first partite the whole dataset into m subgroups. Then in each sub-dataset, for each coordinate, we truncate the data and calculate the mean of the truncated data. Finally, we use the traditional Median of Means (MoM) method in each coordinate, i.e., for each coordinate, we calculate the median among the means of these m subgroups, and add Gaussian noise to ensure (ϵ, δ)-DP.

Theorem 3. For any 0 < ϵ, δ < 1, Algorithm 3 is (ϵ, δ)-DP. Moreover, assume each data xi ∼ D where the distribution D satisfies: (1) D has the mean µ ∈ Rd and ‖µ‖2 ≤ R = O(1); (2) Exi∼D|[xi]j|1+v ≤ u for each j ∈ [d]. Then for any given truncation parameter τ ∈ R and failure probability ξ, with probability at least 1 - ξ, Algorithm 3 outputs a private mean estimator ˜ µ ∈ Rd such that,

where the Big-O notation omits the term of R.

Remark 2. To privately estimate the mean of heavy-tailed distributions with bounded second moments, in general there are three approaches: The first one is directly using one dimensional private estimator to each coordinate [Wang et al., 2020b; Wang et al., 2020a]. However, the bound of this approach is only sub-optimal. [Kamath et al., 2020] proposes a method which aggressively truncate the distribution around a

Algorithm 3 DP High-Dimension Mean Estimator DPHDMEϵ,δ,τ (X) Input: Data samples X = {xi}ni=1 , xi ∈ Rd . Parameters ϵ, δ, τ Output: ADP mean estimator ˜ µ ∈ Rd 1: m← 4 log ( 2d ξ ). 2: for j ← 1, · · · , d do 3: for k ← 1, · · · ,mdo 4: for i ← (k -1) · nm +1, · · · , k · nm do 5: [x′i ]j ← [xi]j · 1|[xi]j|≤τ 6: end for 7: ̂ µkj ← mn ∑ni=1 [x′i ]j 8: end for 9: ̂ µj ← median(̂ µ1j , ̂ µ2j , · · · , ̂ µmj ) 10: end for 11: ̂ µ ← (̂ µ1, ̂ µ2, · · · , ̂ µd) 12: return ˜ µ ← ̂ µ + ν, where ν ∼ N (0, 8τ 2m2 d ln 1.25 δ ϵ2 n2 Id)

point, and compute the noisy empirical mean. However, their theoretical guarantee only holds with constant probability. Instead of these two approaches, here we adopt the idea of the third one, which is a private version of the MoM method and was recently proposed by [Kamath et al., 2021]. However, there are two crucial differences: First, the truncation step is quite different, [Kamath et al., 2021] truncates each data into an interval [a, b], i.e., for the sample x, if x > b then we will let x′ = b and if x < a then we will let x′ = a. However, our approach could be seen as thresholding, i.e., when |x| > τ then x′ = 0. Second, to get the theoretical guarantee, [Kamath et al., 2021] needs Lemma 4.4 in [Kamath et al., 2020] (or Lemma A.2 in [Kamath et al., 2021]), which only holds for data distributions that have at least bounded second moment. Thus, we cannot adopt their approach in our setting. To overcome the challenge we provide the following lemma on the concentration of heavy-tailed distributions, which is the key to prove Theorem 3.

Theorem 4. Let x1, x2, · · · , xn ∈ R be i.i.d. random variables with bounded (1 + v)-th moment, i.e., for ∀i ∈ [n], we have E[|xi|1+v ] ≤ M for some constant v ∈ (0, 1]. Let µ be E[xi], ∀i ∈ [n]. Then

Based on Theorem 3 and the convergent rate of PGD, we give the accuracy guarantees of Algorithm 2 by considering both general convex and strongly convex loss functions.

Theorem 5. For any 0 < ϵ, δ < 1, Algorithm 2 is (ϵ, δ)-DP. Theorem 6 (General Convex Case). Suppose we have a DPSCO problem satisfying Assumption 1. Taking T = R2 ϵ2 n2 τ 2 d4 ,

where the Big- ˜O notation omits all the logarithmic terms and R.

Theorem7 (Strongly Convex Case). Suppose we have a DPSCO problem satisfying Assumption 1, and additionally the loss function ℓ(·, x) is α-strongly convex for every x ∈ X. Taking parameters T = log ( (α+β)G αβ )/ log ( α2+β2+αβ (α+β)2 ), ) 1

η = 1 α+β and τ = ( uϵn d 3 2 √T 1+v in Algorithm 2, then the output wpriv = wT satisfies

where the Big- ˜O notation omits all the logarithmic terms and R.

˜O( d3/2 nϵ ) for convex and strongly convex case respectively. These bounds are consistent with the best known result in [Kamath et al., 2021]. However, the methods in [Kamath et al., 2021] cannot be extended to the case where v ∈ (0, 1). In the following we show that the term of O( 1 (ϵn) v 1+v ) and

O( 1 (ϵn) 2v 1+v ) is optimal for convex and strongly convex loss in (ϵ, δ)-DP respectively. Moreover, we provide lower bounds in the ϵ-DP model.

Theorem 8 (Lower Bound of Strongly Convex Loss). Assume W is the unit norm ball. For any v ∈ (0, 1], there exists a strongly convex and smooth loss function ℓ : W × Rd ↦→ R such that, for any ϵ-DP algorithm A, there is a distribution D over Rd such that for any w, supj∈[d] Ex∼D[|∇jℓ(w, x)|1+v ] ≤ u, the output wpriv of A satisfies the following if n ≥ Ω(u 1v d 1+3v 2v /ϵ)

For any (ϵ, δ)-DP algorithm A with ϵ ≪ log 1δ , there is a distribution D over Rd such that for any w, supj∈[d] Ex∼D[|∇jℓ(w, x)|1+v ] ≤ u, its output wpriv satis-

fies the following if n ≥ Ω(u 1v √log 1δ d 1+2v 2v /ϵ)

Theorem 9 (Lower Bound of Convex Loss). Assume W is the unit norm ball. For any v ∈ (0, 1], there exists a convex and smooth loss function ℓ : W×Rd ↦→R such that, for any ϵ-DP algorithm A, there is a distribution D over Rd such that for any w, supj∈[d] Ex∼D[|∇jℓ(w, x)|1+v ] ≤ u, its output wpriv satisfies the following when n ≥ Ω(d/ϵ)

For any (ϵ, δ)-DP algorithm A with ϵ ≪ log 1δ , there is a distribution D over Rd such that for any w, supj∈[d] Ex∼D[|∇jℓ(w, x)|1+v ] ≤ u, its output wpriv satisfies the following when n ≥ Ω(√d log 1δ /ϵ)

When v = 1, all the rates in Theorem 8 and 9 match the lower bounds in [Kamath et al., 2021]. Thus, our results can be seen as extensions of the previous results. In Theorem 8, the gap between the rates in ϵ-DP and (ϵ, δ)-DP is O(d v 1+v ), while it is O(d v 2(1+v) ) in Theorem 9. This is quite different with the case when the loss is Lipschitz [Bassily et al., 2014]. To prove the lower bounds, we first reduce the problem to mean estimation, and then we use the private version of Fano's lemma in [Acharya et al., 2021; Kamath et al., 2021], based on the packing of distributions in [Barber and Duchi, 2014].

6 Differentially Private Sparse Learning with Heavy-tailed Data

In the previous section, we studied the general case of DPSCO under the assumption that the distribution of each coordinate of the loss gradient has (1 + v)-th moment. However, one weakness of our previous results is that, all the error bounds are in the form of O(Poly(d, 1n , 1ϵ )), which indicates that the error will be large in the high dimensional case where d ≫ n. Moreover, we also showed that in general these polynomial dependencies are unavoidable. Thus, to address the high dimensionality issue, in this section, we focus on some special cases. Specifically, we will study the problem of DP-SCO under sparsity constraints, which is also called DP sparse learning, i.e., W is defined as W = {w : ‖w‖0 ≤ s∗ }. We note that such a formulation encapsulates several important problems such as the ℓ0-constrained linear/logistic regression [Bahmani et al., 2013]. In this section, unlike the previous results on DP sparse learning which need strong assumptions on data distribution, we study the problem under the assumption that the gradient has only (1 + v)-th moments. We first introduce some definitions to the loss functions, which are commonly used in previous research on sparse learning.

Definition 6 (Restricted Strong Convexity, RSC). A differentiable function f(x) is restricted µr-strongly convex with parameter r if for any x, x′ with ‖x - x′ ‖0 ≤ r, we have f (x) - f (x′ ) - 〈∇f (x′ ), x - x′ 〉 ≥ µr 2 ‖x-x′ ‖22 .

Definition 7 (Restricted Strong Smoothness, RSS). A differentiable function f(x) is restricted γs-strongly smooth with parameter r if for any x, x′ with ‖x - x′ ‖0 ≤ r, we have f (x) - f (x′ ) - 〈∇f (x′ ), x - x′ 〉 ≤ γr 2 ‖x-x′ ‖22 .

Note that RSC and RSS are weaker than the strong convexity and smoothness. Next we propose the assumptions that will be used in this section.

Assumption 2. We assume that the objective function LD(·) is µr-RSC and ℓ(w, x) is γr-RSS with parameter r = 2s+s∗ ,

Algorithm 4 Peeling(v, s, ϵ, δ, λ)[Cai et al., 2019]

Input: A vector v ∈ Rd of a dataset X, sparsity s, privacy parameter ϵ, δ, and noise scale λ.

• 1: Initialize S = ∅.
• 2: for i ← 1, · · · , s do
• 3: Generate wi ∈ Rd with wi,1, · · · ,wi,d ∼ Lap( 4λ√2s log 1 δ ϵ ).
• 4: Append j ∗ = argmaxj∈[d]\S |vj| +wi,j to S.
• 5: end for
• 6: Generate ˜ w ∈ Rd with ˜ w1, · · · , ˜ wd ∼ Lap( 4λ√2s log 1 δ ϵ ).
• 7: return vS + ˜ wS .

where s = O(( γr µr )2 s∗ ). We also assume for any w ∈ W′ , the distribution of each coordinate of the gradient of the loss function has bounded (1 + v)-th (raw) moment with some v ∈ (0, 1], i.e., for each j ∈ [d], E[|∇jℓ(w, x)|1+v ] ≤ u, where W′ = {w|‖w‖0 ≤ s}.

There are many problems satisfying Assumption 2, e.g., mean estimation and ℓ2-norm regularized generalized linear loss where LD(w) = E[ℓ(〈w, x〉)] + λ2 ‖w‖22 . If |ℓ′ (·)| ≤ O(1), |ℓ′′ (·)| ≤ O(1) (such as the logistic loss) and xj has bounded (1 + v)-th moment, then we can see that it satisfies Assumption 2.

Our method can be found in Algorithm 5, which is built upon the ideas of our previous private one dimensional mean estimator for heavy-tailed distributions and the Iterative Hard Thresholding method [Blumensath and Davies, 2009]. In detail, in each iteration, we first perform the truncation step to each coordinate of the gradient of the loss to get onedimensional mean estimator. Next, unlike the one dimensional private mean estimator in Algorithm 1 where we add Gaussian noise to the mean of the truncated gradients, here we privately select top s indices via the Peeling mechanism (shown in Algorithm 4). This is due to that if we use Algorithm 1 to each coordinate, then the magnitude of the noise we add will depend on polynomial of d, which is large. However, using the Peeling mechanism will introduce an error that only depends on polynomial of s and log d. In the following we show the theoretical guarantee of our algorithm.

Theorem 10. For any 0 < ϵ, δ < 1, Algorithm 5 is (ϵ, δ)-DP. Moreover, under Assumption 2, given any failure probability ξ, if we set T = ˜O( γr µr logn), s = O(( γr µr )2 s∗ ), η0 = 2 3γr

where the Big-O notation omits γr and µr.

Compared with the results in Section 5, we can see that in Theorem 10 the bound is only logarithmic in d and polynomial in s∗ , 1ϵ and 1n , which means it is more suitable to the high dimensional case.

Algorithm 5 Heavy-Tailed Private Sparse Optimization

Input: Data samples X = {xi}ni=1 ⊆ Rd , parameters s, T, η, initial s-sparse parameter w1 , privacy parameter ϵ, δ.

Output: Private minimizere wpriv

• 1: Split X into T parts {Xt}Tt=1 with |Xt| =m= n T .
• 2: for t ← 1, · · · , T do.
• 3: for each dimension j ∈ [d] do

4:

5:

∇j ℓ′ (wt , x) ← ∇j ℓ(wt , x)1|∇j ℓ(wt ,x)|≤B

• 6: end for 7: end for
• 8: Get the robust gradient estimator ˜ gt (wt , Xt):

• 9: Denote wt+0.5 ←wt - η0˜ gt (wt , Xt)
• 10: Let wt+1 ←Peeling(wt+0.5 , s, ϵ, δ, 2Bη0 m ).
• 11: end for
• 12: return wpriv ←wT+1 .

Remark 3. For DP sparse learning with Lipschitz loss or regular data, [Wang and Xu, 2019] provided an upper bound of ˜O( s∗ n2 ϵ2 ). Moreover, for high dimensional sparse mean estimation and Generalized Linear Model (GLM) with the Lipschitz loss and sub-Gaussian data, [Cai et al., 2020; Cai et al., 2019] provided optimal rates of ˜O( s∗ log d n + (s∗ log d)2 (nϵ)2 ). We can see that compared with these results, the error bound now becomes ˜O( (s∗ ) 1+2v 1+v u 2 1+v (nϵ) 2v 1+v ) due to data irregularity. When

v = 1, the error bound now becomes to ˜O( (s∗ ) 3 2 u (nϵ) ), which matches the result in [Hu et al., 2021]. Thus, our result can be seen as a generalization of the previous ones.

One open question is whether we can further improve the rate of error in Theorem 10. In the following we show that the bound is optimal up to a factor of ˜O(s∗ ).

Theorem 11. For W = {w ∈ Rd : ‖w‖0 ≤ s∗ } and any v ∈ (0, 1], there exists a strongly convex and smooth loss function ℓ : W × Rd ↦→ R such that, for any ϵ-DP algorithm A, there is a distribution D over Rd such that for any w, supj∈[d] Ex∼D[|∇jℓ(w, x)|1+v ] ≤ u, its output wpriv satisfies the following when n ≥ Ω(s∗ log d/ϵ)

For any (ϵ, δ)-DP algorithm A with ϵ ≪ log 1δ , there is a distribution D over Rd such that for any w, supj∈[d] Ex∼D[|∇jℓ(w, x)|1+v ] ≤ u, its output wpriv satisfies the following when n ≥ Ω(√s∗ log d log 1δ /ϵ)

for each data sample x ∈ Xt do

Acknowledgments

• Di Wang and Yulian Wu were support in part by the baseline funding BAS/1/1689-01-01, funding from the CRG grand URF/1/4663-01-01 and funding from the AI Initiative REI/1/4811-10-01 of King Abdullah University of Science and Technology (KAUST). Xiuzhen Cheng and Youming Tao were partially supported by National Key R&D Program of China (No. 2019YFB2102600).