On the Computational Complexity of Private High-dimensional Model Selection

Authors: saptarshi  roy
zehua  wang
ambuj  tewari
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 19:01:53

Abstract: We consider the problem of model selection in a high-dimensional sparse linear regression model under privacy constraints. We propose a differentially private (DP) best subset selection method with strong statistical utility properties by adopting the well-known exponential mechanism for selecting the best model. To achieve computational expediency, we propose an efficient Metropolis-Hastings algorithm and under certain regularity conditions, we establish that it enjoys polynomial mixing time to its stationary distribution. As a result, we also establish both approximate differential privacy and statistical utility for the estimates of the mixed Metropolis-Hastings chain. Finally, we perform some illustrative experiments on simulated data showing that our algorithm can quickly identify active features under reasonable privacy budget constraints.

Main paper: 
 introduction: 
 In this paper, we consider the problem of private model selection in high-dimensional sparse regression which has been one of the central topics in statistical research over the past decade. Consider n observations {(x i , y i )} n i=1 ⊆ X × Y following the linear model:EQUATIONwhere {x i } i∈ [n] are fixed p-dimensional feature vectors, {w i } i∈[n] are i.i.d. mean-zero σ-sub-Gaussian noise, i.e., E exp(λw i ) ≤ exp(λ 2 σ 2 /2) for all λ ∈ R and i ∈ [n], and the signal vector β ∈ R p is unknown but is assumed to have a sparse support. In matrix notation, the observations can be represented as y = Xβ + w, where y = (y 1 , . . . , y n ) ⊤ , X = (x 1 , . . . , x n ) ⊤ , and w = (w 1 , . . . , w n ) ⊤ . We consider the standard high-dimensional sparse setup where n < p, and possibly n ≪ p, and the vector β is sparse in the sense that ∥β∥ 0 := p j=1 1(β j ̸ = 0) = s, which is much smaller than p. The main goal of variable selection is to identify the active set γ * := {j : β j ̸ = 0}.For the past two decades, there has been ample work on model selection problem in the non-private setting for ℓ 1 -penalized methods [57, 54, 46, 23] , concave regularized methods [55, 53, 12, 16] , ℓ 0penalized/constrained methods [13, 1, 37, 39] in high-dimensional setting. On the computational side, recent advancements related to mixed integer optimization (MIO) in [4, 5] and [20] have pushed the computational barrier of best subset selection (BSS) in terms of solving problems of large dimensions (large p), and consequently, simulation studies in [19] have revealed the improved performance of BSS over its computational surrogates like LASSO, SCAD, and MCP. Despite these theoretical and computational advancements related to BSS, to the best of our knowledge, there is no computationally efficient private algorithmic framework for BSS for highdimensional sparse regression setup (1) . This is especially surprising as private model selection is important in many contemporary applications involving sensitive data including genetics [21] , neuroimaging [33] , and computer vision [56] . One major reason for this could be the lack of DP mechanisms for MIO problems which restricts us from exploiting the MIO formulation of BSS introduced in [4] . Secondly, the apparent computational burden stemming from the requirement of exponentially large numbers of search queries in private BSS has eluded the majority of the machine learning and statistics community. In this paper, we address the latter issue by mainly focusing on the utility and computational complexity of BSS under privacy constraints. To be specific, we make the following contributions listed below:1. We adopt the exponential mechanism [31] to design a DP BSS algorithm, and we establish its good statistical or utility guarantee under high-privacy regime whenever β min := min j∈γ * |β j | ≳ σ{(s log p)/n} 1/2 . 2. Under the low-privacy regime, we show that accurate model recovery is possible whenever β min ≳ σ{(log p)/n} 1/2 , which is the minimax optimal β min requirement for model recovery under non-private setting. Therefore, this paper points out an inflection phenomenon in the signal strength requirement for the model consistency across different privacy regimes. 3. In addition, we design an MCMC chain that converges to its stationary distribution that matches the sampling distribution in the exponential mechanism. As a consequence, the model estimator generated by the MCMC also enjoys (approximate) DP. Furthermore, under certain regularity conditions on the design, we show that the MCMC chain enjoys a polynomial mixing time in (n, p, s) to the stationary distribution with good utility guarantee.In summary, this paper proposes a DP version of BSS that generates a private model estimator of γ * with strong model recovery property within polynomial time in the problem parameters n, p, s. In the next section, we will discuss some prior related works on DP model selection and discuss some of their limitations.In the past decades, there has been a considerable amount of work studying DP sparse regression problems. However, most of these works focus either on empirical risk minimization [25, 44, 26, 48] or establishing ℓ 2 -consistency rate [47, 6] which are not directly related to the task of model selection.To the best of our knowledge, there are only three works considering the problem of variable selection in sparse regression problems under the DP framework, [27, 45] , and [29] . Table 1 shows a clear comparison between those methods and our method. [27] proposed two algorithms under sparse regression setting. One of them is based on the exponential mechanism, which is known to be computationally inefficient due to exponentially large numbers of search queries. However, they do not analyze the algorithm under the model selection framework. Moreover, for the privacy analysis, they assume that the loss functions are bounded over the space of sparse vectors, which is rather restrictive in the linear regression setting. In comparison, our paper provides a solid model recovery guarantee (Theorem 3.5) for a similar exponential mechanism without using the bounded loss assumption. Furthermore, under a slightly stronger assumption, we design a computationally efficient MCMC algorithm that also enjoys desirable utility similar to the exponential mechanism (Theorem 4.3) under DP framework. The other algorithm in [27] is based on the resample-andaggregate framework [36, 42] . Although computationally efficient, this method requires sub-optimal β min condition compared to Theorem 3.5. In [45] , the authors introduced two concepts of stability for LASSO and proposed two PTR-based (propose-test-release) algorithms for variable selection. However, these methods have nontrivial probabilities of outputting the null (no result), which is undesirable in practice. Also, the support recovery probabilities for these methods do not approach 1 with a growing sample size even under the strong irrepresentability condition [57] on the design matrix. In [29] , the authors proposed to use the Akaike information criterion or Bayesian information criterion coupled with the exponential mechanism to choose the proper model. However, the runtime of this algorithm is exponential and also requires stronger β min condition. As mentioned earlier, in this paper, we show that our proposed MCMC algorithm is both computationally efficient and produces approximate DP estimates of γ * with a strong utility guarantee under a better β min condition. One may also apply sparse vector techniques (SVT) to choose important features [43] . In this case, each feature can be associated with an appropriate choice of score function, and then apply SVT to choose the relevant features. However, the choice of the score function in high-dimensional sparse regression cases remains unclear, and moreover, it is also known that the exponential mechanism enjoys better accuracy compared to SVT [30] under such an offline setting.  
differential privacy: 
 Differential privacy requires the output of a randomized procedure to be robust with respect to a small perturbation in the input dataset, i.e., an attacker can hardly recover the presence or absence of a particular individual in the dataset based on the output only. It is important to note that differential privacy is a property of the randomized procedure, rather than the output obtained.In this section, we will formalize the notion of differential privacy. Consider a dataset D := {z 1 , . . . , z n } ∈ Z n consisting of n datapoints in the sample space Z. A randomized algorithm A maps the dataset D to A(D) ∈ O, an output space. Thus, A(D) is a random variable on the output space O.For any two datasets D and D ′ , we say they are neighbors if |D∆D ′ | = 1. We can now formally introduce the definition of differential privacy. Definition 2.1 ((ε, δ)-DP, [9] ). Given the privacy parameters (ε, δ) ∈ R + × R + , a randomized algorithm A(•) is said to satisfy the (ε, δ)-DP property ifEQUATIONfor any measurable event K ∈ range(A) and for any pair of neighboring datasets D and D ′ .In the above definition, the probability is only with respect to the randomness of the algorithm A(•), and it does not impose any condition on the distribution of D or D ′ . If both ε and δ are small, then Definition 2.1 essentially entails that distribution of A(D) and A(D ′ ) are almost indistinguishable from each other for any choices of neighboring datasets D and D ′ . This guarantees strong privacy against an attacker by masking the presence or absence of a particular individual in the dataset. As a special case, when δ = 0, the notion of DP in Definition 2.1 is known as the pure differential privacy.For any DP procedure, a specific randomized procedure A must be designed that takes a database D ∈ Z n as input and returns an element of the output space O while satisfying the condition in (2) . Several approaches exist that are generic enough to be adaptable to different tasks, and which often serve as building blocks for more complex ones. A few popular examples include the Laplace mechanism [11] , Gaussian mechanism [10] , and Exponential mechanism [31] . We only provide the details of the last technique, since the other two techniques are out-of-scope for the methods and experiments in this paper.Intuitively, sensitivity quantifies the effect of any individual in the dataset on the outcome of the analysis. The score function u(•, •) is called data monotone if the addition of a data record can either increase (decrease) or remain the same with any outcome, e.g., u(α, D) ≤ u(α, D ∪ {z}). Next, we have the following result. Lemma 2.2 ( [8, 31] ). Exponential mechanism A E (D) that outputs samples from the probability distributionEQUATIONpreserves (2ε, 0)-differential privacy. If u(•, •) is data monotone, then we have (ε, 0)-differential privacy.In general, if the S is too large, the sampling from the distribution could be computationally inefficient. However, we show below that the special structure of the linear model (1) allows us to design an MCMC chain that can generate approximate samples efficiently from the distribution (3) for privately solving BSS under an appropriately chosen score function. 
best subset selection: 
 We briefly review the preliminaries of BSS, one of the most classical variable selection approaches.For a given sparsity level s, BSS solves for β best ( s) := arg min θ∈R p ,∥θ∥ 0 ≤ s ∥y -Xθ∥ β best ( s)] j ̸ = 0}. For a subset γ ⊆ [p], define the matrix X γ := (X j ; j ∈ γ). Let Φ γ := X γ (X ⊤ γ X γ ) -1 X ⊤γ be orthogonal projection operator onto the column space of X γ . Also, define the corresponding residual sum of squares (RSS) for model γ as L γ (y, X) := y ⊤ (I n -Φ γ )y. With this notation, the γ best ( s) can be alternatively written asEQUATIONLet X γ be the matrix comprised of only the columns of X with indices in γ, and Φ γ denotes the orthogonal projection matrix onto the column space of X γ . In addition, let Σ := n -1 X ⊤ X be the sample covariance matrix and for any two sets γ 1 , γ 2 ⊂ [p], Σ γ1,γ2 denotes the submatrix of Σ with row indices in γ 1 and column indices in γ 2 . Finally, define the collectionA s := {γ ⊂ [p] : γ ̸ = γ * , |γ| = s}, and for γ ∈ A s write Γ(γ) = Σ γ * \γ,γ * \γ -Σ γ * \γ,γ Σ -1 γ,γ Σ γ,γ * \γ . Then, it follows that β ⊤ γ * \γ Γ(γ)β γ * \γ is equal to the residualized signal strength n -1 ∥(I n -Φ γ )X γ * \γ β γ * \γ ∥ 2 2 . Therefore, β ⊤ γ * \γ Γ(γ)β γ * \γquantifies the separation between γ and the true model γ * . Ideally, a larger value of the quantity will help BSS to discriminate between γ * and any other candidate model γ. More details on this can be found in [40] . Now we are ready to introduce the identifiability margin that characterizes the model discriminative power of BSS.The discussion in Section 3 motivates us to define the following identifiablity margin:EQUATION)As mentioned earlier, the quantity m * ( s) captures the model discriminative power of BSS. To add more perspective, note that if the features are highly correlated among themselves then it is expected that m * ( s) is very close to 0. Hence, any candidate model γ is practically indistinguishable from the true model γ * which in turn makes the problem of exact model recovery harder. On the contrary, if the features are uncorrelated then m * ( s) becomes bounded away from 0 making the true model γ * easily recoverable. For example, [17] showed that under the knowledge of true sparsity, i.e., when s = s, the conditionEQUATIONis sufficient for BSS to achieve model consistency. One can also view m * (s) as a quantifier of the coupled effect of model correlation and signal strength. If we define the minimum and maximum eigenvalues over all models to be λ * = min γ∈As λ min (Γ(γ)) and λ * = max γ∈As λ max (Γ(γ)) respectively, then it follows that λ * β 2 min ≤ m * (s) ≤ λ * β 2 min . Therefore, it suffices to have β min ≳ σ{(log p)/(nλ * )} 1/2 in order to satisfy condition (6) . In this case, λ * captures the degree of model correlation, and β min is the minimum signal strength. Similar to our previous discussion, if there is high collinearity in the model, λ * will be typically small, and BSS needs a large value of β min to identify the true model γ * . On the other hand, if β min is too small for a given level of model correlation, i,e, if β min ≪ σ{(log p)/(nλ * )} 1/2 , then also BSS fails to achieve model consistency as it is hard to identify active features under the presence of weak signals [17, Theorem 2.1] . As we will see in the next section, the DP BSS algorithm also requires a margin condition similar to (6) to ensure model recovery, and this is indeed an indispensable condition as it is needed even in non-private case.In order to privatize the optimization problem in (4), we will adopt the exponential mechanism discussed in Section 2.2. In particular, for a tuning parameter K > 0, we consider the score function u K (γ; X, y) := -minθ∈R s :∥θ∥ 1 ≤K ∥y -X γ θ∥ 2 2 ,and for a given privacy budget ε > 0, we sample γ ∈ A s from the distributionEQUATIONAs we are concerned with the exact recovery γ * , from here on we assume s = s. The above algorithm is essentially the same as Algorithm 4 in [27] ; however, they do not introduce the extra ℓ 1 -constraint on the parameter space. Instead, their algorithm needs the loss-term (y -x ⊤ γ θ) 2 to be bounded by a constant for every possible choice of x, y, γ and θ. This assumption is not true in general for the squared error loss, and to remedy this issue, we introduce the extra ℓ 1 -constraint in the score function. This is a common strategy that is used to guarantee worst-case sensitivity bound and similar methods also have been adopted in [29, 6] to construct private estimators. Next, we present the following lemma that shows the data-monotonicity of the proposed score function. Lemma 3.1. The score function u K (γ; •) in (7) is data monotone.Therefore, Lemma 2.2 automatically guarantees that the above procedure is (ε, 0)-DP. However, in practice, we need an explicit form for ∆u K to carry out the sampling method, and it is also needed to analyze the utility guarantee of the exponential mechanism. To provide a concrete upper bound on the global sensitivity of u K (•; •), we make the following boundedness assumption on the database: Assumption 3.2. There exists positive constants r, x max such that sup y∈Y |y| ≤ r, sup x∈X ∥x∥ ∞ ≤ x max .Under this assumption, the following lemma provides an upper bound on the global sensitivity of the score function along with the DP guarantee. Lemma 3.3 (Sensitivity bound and DP). Under Assumption 3.2, the global sensitivity ∆u K is bounded by ∆ K := (r + x max K) 2 . Therefore, the exponential mechanism (7) with ∆u K replaced by ∆ K satisfies (ε, 0)-DP.The above lemma provides an upper bound on the global sensitivity of the score function rather than finding the exact value of it. However, to guarantee (ε, 0)-DP property of exponential mechanism, it suffices to use the upper bound of ∆u K in (7) . Now we will shift towards the utility analysis of the proposed exponential mechanism. First, we require some technical assumptions. Assumption 3.4. We assume the following hold:(a) There exists positive constants b max such that ∥β∥ 1 ≤ b max .(b) There exists positive constants κ -, κ + such thatEQUATIONfor all γ ∈ A s ∪ {γ * }.(c) The true sparsity level s follows the inequality s ≤ n/(log p).Assumption 3.4(a) tells that the true parameter β lies inside a ℓ 1 -ball. Similar boundedness assumptions are fairly standard in privacy literature [49, 29, 6] . Assumption 3.4(b) is a well-known assumption in the high-dimensional literature [54, 22, 32] which is known as the Sparse Riesz Condition (SRC). Finally, Assumption 3.4(c) essentially assumes that the s = o(n), i.e., sparsity grows with a sufficiently small rate compared to the sample size n. Theorem 3.5 (Utility gurantee). Let the conditions in Assumption 3.2 and Assumption 3.4 hold. SetK ≥ {(κ + /κ -)b max + (8x max /κ -)σ} √ s.Then, under the data generative model (1), there exist universal positive constants c 1 , C 1 such that wheneverEQUATIONwith probability at least 1 -c 1 p -2 we have π(γ * ) ≥ 1 -p -2 .Cost of privacy. Theorem 3.5 essentially says that whenever the identifiability margin is large enough, the exponential mechanism outputs the true model γ * with high probability. Note that ∆ K /σ 2 = Ω(s). In the low privacy regime, i.e., for ε > ∆ K /σ 2 we only require m * (s) ≳ σ 2 (log p)/n to achieve model consistency and this matches with the optimal rate for model consistency of non-private BSS. Note that, the margin condition does not depend at all on ε in this regime. In contrast, in a high privacy regime, i.e., for ε < ∆ K /σ 2 , Condition (9) essentially demands m * (s) ≳ σ 2 (s log p)/(nε) to achieve model consistency. Thus, in a high privacy regime, we pay an extra factor of (s/ε) in the margin requirement. Remark 3.6. The failure probability in Theorem 3.5 can be improved to O(p -M ) for any arbitrary integer M > 2. However, we have to pay a cost in the universal constant C 1 in terms of a multiplicative constant larger than 1. Remark 3.7. Under Assumption 3.4(b), it follows that λ * ≥ κ -. Therefore, it suffices to havemin j∈γ * β 2 j ≥ C1σ 2 κ- max 1, ∆ K /(εσ 2) log p n in order to hold condition (9) . Therefore, in highprivacy regime, our method requires min j∈γ * |β j | ≳ σ{(s log p)/(nεκ -)} 1/2 . In contrast, under the low-privacy regime, we retrieve the optimal requirement min j∈γ * |β j | ≳ σ{(log p)/(nκ -)} 1/2 . 
efficient sampling through mcmc: 
 In this section, we will propose an efficient sampling method to generate approximate samples from the distribution (7) . One of the challenges of sampling methods in high-dimension is their high computational complexity. For example, the distribution in (7) places mass on all p s subsets of [p], and it is practically infeasible to sample γ from the distribution as we have to essentially explore over an exponentially large space. This motivates us to resort to sampling techniques based on MCMC, through which we aim to obtain approximate samples from the distribution in (7) . Past works on MCMC algorithms for Bayesian variable selection can be divided into two main classes -Gibbs sampler [15, 24, 34] and Metropolis-Hastings [18, 28] . In this paper, we focus on a particular form of Metropolis-Hastings updates.In general terms, Metropolis-Hastings (MH) random walk is an iterative and local-move based method involving three steps:1. Given the current state γ, construct a neighborhood N (γ) of proposal states.Choose a new state γ ′ ∈ N (γ) according to some proposal distribution F(γ, •) over the neighborhood N (γ).3. Move to the new state γ ′ with probability R(γ, γ ′ ), and stay in the original state γ with probability 1 -R(γ, γ ′ ), where the acceptance probability is given byR(γ, γ ′ ) = min 1, π(γ ′ )F(γ ′ , γ) π(γ)F(γ, γ ′ ) ,where π(•) is same as in Equation ( 7).This procedure generates a Markov chain for any choice of the neighborhood structure N (γ) with the following transition probability:P MH (γ, γ ′ ) =    F(γ, γ ′ )R(γ, γ ′ ), if γ ′ ∈ N (γ), 1 -γ ′ ̸ =γ P MH (γ, γ ′ ), if γ ′ = γ, 0, otherwise.The specific form of Metropolis-Hastings update analyzed in this paper is obtained by following the double swap update scheme to update γ.Double swap update: Let γ ∈ A s ∪ {γ * } be the initial state. Choose an index pair (k, ℓ) ∈ γ × γ c uniformly at random. Construct the new state γ ′ by setting γ ′ = γ ∪ {ℓ} \ {k}.The above scheme can be viewed as a general MH update scheme when N (γ) is the collection of all models γ ′ which can be obtained by swapping two distinct coordinates of γ and γ c respectively. Thus, lettingd H (γ, γ ′ ) = |γ \ γ ′ | + |γ ′ \ γ| denote the Hamming distance between γ and γ ′ , the neighborhood is given by N (γ) = {γ ′ | d H (γ, γ ′ ) = 2, ∃ (k, ℓ) ∈ γ × γ c such that γ ′ = γ ∪ {ℓ} \ {k}}.With this definition, the transition matrix of the previously described Metropolis-Hastings scheme can be written as follows:EQUATIONLet C be a Markov chain on the discrete space S with a transition probability matrix P ∈ R |S |×|S | with stationary distribution ν. Throughout our discussion, we assume that C is reversible, i.e., it satisfies the balanced condition ν(γ)P(γ, γ ′ ) = ν(γ ′ )P(γ ′ , γ) for all γ, γ ′ ∈ S . Note that the previously described transition matrix P MH in (10) satisfies the reversibility condition. It is convenient to identify a reversible chain with a weighted undirected graph G on the vertex set S , where two vertices γ and γ ′ are connected if and only if the edge weight Q(γ, γ ′ ) := ν(γ)P(γ, γ ′ ) is strictly positive. For γ ∈ S and any subset S ⊆ S , we write P(γ, S) = γ ′ ∈S P(γ, γ ′ ). If γ is the initial state of the chain, then the total variation distance to the stationary distribution after t iterations is∆ γ (t) = P t (γ, •) -ν(•) TV := max S⊂S P t (γ, S) -ν(S) .The η-mixing time is given byEQUATIONwhich measures the number of iterations needed for the chain to be within distance η ∈ (0, 1) of the stationary distribution.Privacy of MCMC estimator: Now, we will show that once the MH chain in (10) has mixed with its stationary distribution π(•) defined in (7) , the model estimators at each iteration will enjoy approximate DP. To fix the notation, let γ t be the tth iteration of the MH chain in (10) . Then, we have the following useful lemma:Lemma 4.1. The model estimator γ τη is (ε, δ)-DP with δ = η(1 + e ε ).The above lemma shows that smaller η entails a better privacy guarantee for a fixed level ε as δ decreases with η. Therefore, allowing more mixing of the chain will provide better privacy protection. However, this raises a concern about how long a practitioner must wait until the chain archives η-mixing. In particular, it is important to understand how τ η scales in the difficulty parameters of the problem, for example, the dimension of the parameter space and sample size. In our case, we are interested in the covariate dimension p, sample size n, sparsity s, and the privacy parameter ε. In the next section, we will show that the chain with transition matrix (10) enjoys rapid mixing, meaning that the mixing time τ η grows at most at a polynomial rate in p, s and the sample size n.We now turn to develop sufficient conditions for MH scheme (10) to be rapidly mixing. To this end, we make a technical assumption on the design matrix. Essentially, the following assumption controls the amount of correlation between active features and spurious features. Assumption 4.2. For everyγ ′ ∈ A s \ {γ * }, there exists k / ∈ γ * ∪ γ ′ such that max j∈γ * \γ ′ X ⊤ j (I n -Φ γ ′ )X k ∥(I n -Φ γ ′ )X k ∥ 2 ≤ b -1 max (κ -C 1 σ 2 /2) log p,where C 1 is the same universal positive constant as in Theorem 3.5.First, note that Assumption 4.2 basically controls the length of the projection of the feature X j on the unit vectoru k := (I n -Φ γ ′ )X k / ∥(I n -Φ γ ′ )X k ∥ 2 .Therefore, the above inequality restricts the correlation between an active feature X j and the spurious scaled feature u k from being too large. To this end, we emphasize that stronger assumptions on model correlation (on top of the SRC condition) are common in literature for establishing the computational efficiency of Bayesian variable selection methods involving MH algorithm. For example, to show the computational efficiency MH algorithm under Zellner's g-prior, [51] assumesEQUATIONwhere s 0 (larger than s) is a specific tuning parameter of their algorithm that controls the model size.The assumption in the above display is akin to the well-known irrepresentability condition [57] which is a very strong assumption on the design. On a high level, at any given current state γ, Assumption 4.2 or Condition (12) helps to identify a good local move towards the true model γ * in the MH algorithm via deletion of the least influential covariate in γ. Now, we present our main result for the mixing time of MCMC. EQUATIONthere exist universal positive constants c 2 , C 2 such that the mixing time τ η of the MCMC chain (10) enjoys the following with probability at least 1 -c 2 p -2 :EQUATIONwhereΨ = r + (κ + /κ -)b max x max + (σ/κ -)x 2 max 2 .The main technical innovation in the theorem is the double swap updating scheme in the MCMC that allows us to leverage the canonical path ensemble construction argument [41] to prove the bound ( 14) on the mixing time. Essentially, we show that under Assumption 4.2, there exists a canonical path in a specially weighted graph corresponding to the MCMC random walk with low path congestion. The complete proof can be found in Appendix A.5.Regarding the statement of the theorem, note that the margin condition ( 13) is slightly stronger than the margin condition in Theorem 3.5. Under that condition, the above theorem shows that the η-mixing time of the MCMC algorithm designed for approximate sampling from the distribution (7) grows at a polynomial rate in (n, p, s). Recall that according to the previous definition (11) of the mixing time, Theorem 4.3 characterizes the worst-case mixing time, meaning the number of iterations when starting from the worst possible initialization. If we start with a good initial statefor example, the true model γ * would be an ideal though impractical choice -then we can remove the n term in the upper bound in (14) . Therefore, the bound in ( 14) can be thought of as the worst-case number of iterations required in the burn-in period of the MCMC algorithm. Furthermore, it is important to point out that Assumption 4.2 is only needed to ensure the "quick" mixing time of the MCMC chain. It is possible to relax this assumption, however, in that case, the MCMC chain is not guaranteed to mix under polynomial time. Nonetheless, given enough iterations, the chain will indeed converge to the distribution (7) as MH algorithm always generates an ergodic chain that eventually mixes to its stationary distribution.It is interesting to note that Theorem 4.3 suggests that in a large ε regime, the chain mixes slower compared to the small ε regime. The main reason for this is that Theorem 4.3 only relies on worst-case analysis. The intuition is the following: When ε is very large, then the target distribution is essentially fully concentrated on γ * (assuming the score for γ * is highest). Now, the current analysis of Theorem 4.3 does not assume any condition on the initial state of the MCMC chain. It treats the initial state γ 0 as if it is chosen in a completely random manner, i.e., it is the worst case. From this point of view, it is hard for a completely uninformative distribution to converge to a target distribution that is concentrated on a single subset (very informative), and resulting in a longer mixing time. Finally, Theorem 4.3 leads to the following corollary:Corollary 4.4. Let π t denote the distribution of the tth iterate γ t of the MCMC scheme (10) . Then, under the conditions of Theorem 3.5 and Theorem 4.3, there exists a universal constant c 3 > 0 such that for any fixed iterationt such that t ≥ C 2 ps 2 nεΨ -1 κ + b 2 max + log(1/η) , we have π t (γ * ) ≥ 1 -η -p -2 with probability at least 1 -c 3 p -2 .The above corollary is useful in the sense that it provides a quantitative choice of η that yields high utility of the estimator γ t . For example, if we set η = p -2 and ε = O(1), then for any t ≳ ps 2 (nε + log p) the resulting sample γ t will match γ * with probability 1 -c 3 p -2 . Remark 4.5. Similar to Remark 3.6, the failure probability in Theorem 4.3 and Corollary 4.4 can be improved to O(p -M ) for arbitrary large M > 2, but at the cost of paying higher values for the absolute constants C ′ 1 and C 2 . 
numerical experiments: 
 In this section, we will conduct some illustrative simulations. To compare the quality of the DP model estimator, we compare F-score [19] of the estimated model with that of the true model γ * and the BSS estimator. As the actual BSS is computationally infeasible, we use the adaptive best subset selection (ABESS) algorithm [58] as a computational surrogate to BSS. Throughout this section, we assume that the true sparsity s is known, i.e., we provide the knowledge of s to the algorithm. All codes are available at https://github.com/roysaptaumich/DP-BSS.Uniform design. We consider a random design matrix, formed by choosing each entry from the distribution Uniform(-1, 1) in i.i.d. fashion. In detail, we set n = 900, p = 2000, and the sparsity level s = 4. We generate the entries of the noise w independently from Uniform(-0.1, 0.1), and consider the linear model (1) . We choose the design vector β with true sparsity s = 4 and the support set γ * = {j : 1 ≤ j ≤ 4}. We set all the signal strength to be equal, taking the following two forms: (i) Strong signal: β j = 2{(s log p)/n} 1/2 , and (ii) Weak signal:β j = 2{(log p)/n} 1/2 for all j ∈ γ * .Under these setups, we consider the privacy parameter ε ∈ {0.5, 1, 3, 5, 10} which are acceptable choices of ε [35] . Moreover, similar (or larger) choices of ε are common in various applications including US census study [14] , socio-economic study [38] , and industrial applications [3, 52] . For the Metropolis-Hastings random walk, we vary K ∈ {0.5, 2, 3, 3.5} and initialize 10 independent Markov chains from random initializations and record the F-score of the last iteration. We use the CVXPY package [7, 2] for solving the ℓ 1 -constrained optimization problem in the updating step of MCMC. We also track the qualities of the model through its explanatory power for convergence diagnostics. In particular, we calculate the scale factor R γ := y ⊤ Φ γ y/ ∥y∥ 2 2 for each model update along the random walk and compare those with R γ best to heuristically gauge the quality of mixing. More details and a set of comprehensive plots can be found in Appendix D.1 where we also discuss more about the effect of ε and K on the utility. For K = 2, Table 2 shows that F-score increases as ε increases both in the cases of strong and weak signals. In fact, for ε ≥ 3, the performance of the algorithm is on par with the non-private BSS. This is consistent with the inflection phenomenon pointed out in Theorem 3.5 and Corollary 4.4. Furthermore, as expected, we see that for a fixed ε, the F-score is generally higher in the strong signal case. We also carry out experiments under independent Gaussian design. The details and more comprehensive discussion of the findings are deferred to Appendix D.2. In summary, in this case also, our algorithm enjoys greater utility under the strong signal case as shown in Table 2 .Computational resources and license information : All the experiments were performed in the Great Lakes cluster with 16 cores and 10 GB RAM. ABESS package is distributed under GNU General Public License, Version 3. CVXPY package is distributed under Apache License, Version 2.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: The Github Link is provided in the footnote of Page 9.Guidelines:• The answer NA means that paper does not include experiments requiring code. • At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable). • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
conclusion: 
 In this paper, we study the variable selection performance of BSS under the differential privacy constraint. In order to achieve (pure) differential privacy, we adopt the exponential mechanism and establish its high statistical utility guarantee in terms of exact model recovery. Furthermore, for computational efficiency, we design a MH random walk that provably mixes with the stationary distribution within a mixing time of the polynomial order in (n, p, s). We also show that the samples from the MH random walk enjoy approximate DP while retaining a high utility guarantee with experimental underpinnings. In summary, as discussed in Section 1.1, we establish both high utility and efficient computational guarantee for our model selection algorithm under privacy constraints, which is in sharp contrast with the previous works in DP model selection literature. Moreover, the proposed MCMC method is generic enough to adopt in other models beyond linear structures. For example, one can use this technique under the setup of generalized linear models with likelihood loss as the utility function. Therefore, our method can be used in diverse domains including medical studies to fast-track scientific discoveries and promote the practice of responsible AI.To this end, we also point out some of the open problems and future directions. One limitation, of our main result Theorem 3.5 is that it requires the condition min j∈γ * |β j | = Ω( (s log p)/n) in high-privacy regime. It is still an open question whether the extra √ s factor is necessary for model selection. Future research along this line could focus on solving BSS through DP mixed integer optimization (MIO). This would mean an important contribution in this field as commercial solvers like GUROBI or MOSEK would be capable of solving the BSS problems at an industrial scale with high computational efficiency using a general DP framework.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: We do the experiments on simulated data. Therefore, technically, we just have the testing step. We do specify clearly about this step in the simulation section.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [No] Justification: We do not include the error bar for the F-scores as the results are pretty robust across different repetitions.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: We have added the information about the computational resources in Section 5.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes]Justification: We point toward the potential positive impacts of our work in the conclusion sections. To the best of our knowledge, we could not think of any negative societal impact of our work.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: A.1 Proof of Lemma 3.1Recall that u K (γ; X, y) = -L γ,K (X, y) where L γ,K (X, y) := n i=1 (y i -x ⊤ i,γ β) 2 .Therefore, it suffices to show that L γ,K (•, •) is data monotone. Let D = {(x i , y i )} n i=1 and D ′ = D ∪ {x n+1 , y n+1 }. We defineβ n,γ := arg min β:∥β∥ 1 ≤K n i=1 (y i -x ⊤ i,γ β) 2 , β n+1,γ := arg min β:∥β∥ 1 ≤K n+1 i=1 (y i -x ⊤ i,γ β) 2 .Therefore, we have the following inequalities:L γ,K (D ′ ) = n+1 i=1 (y i -x ⊤ i,γ β n+1,γ ) 2 ≥ n i=1 (y i -x ⊤ i,γ β n+1,γ ) 2 ≥ n i=1 (y i -x ⊤ i,γ β n,γ ) 2 = L γ,K (D).The above inequalities conclude the proof.A.2 Proof of Sensitivity Bound (Lemma 3.3)Let (X, y) and ( X, ỹ) be two neighboring datasets with n and n + 1 observation respectively. For a subset γ ∈ A s ∪ {γ * }, consider the OLS estimators as follows:β γ,K := arg min θ:∥θ∥ 1 ≤K ∥y -X γ θ∥ 2 2, and βγ,K := arg min θ:∥θ∥1 ≤K ỹ -Xγ θ 2 2.From the definition of the score function u(γ; X, y), we haveu(γ; X, y) = - n i=1 (y i -x ⊤ i,γ β γ,K ) 2 u(γ; X, ỹ) = - n i=1 (y i -x ⊤ i,γ βγ,K ) 2 -(ỹ n+1 -x⊤ n+1,γ βγ,K ) 2 .By the property of the OLS estimators, we haveu(γ; X, y) -u(γ; X, ỹ) = n i=1 (y i -x ⊤ i,γ βγ,K ) 2 + (ỹ n+1 -x⊤ n+1,γ βγ,K ) 2 - n i=1 (y i -x ⊤ i,γ β γ,K ) 2 ≤ n i=1 (y i -x ⊤ i,γ β γ,K ) 2 + (ỹ n+1 -x⊤ n+1,γ β γ,K ) 2 - n i=1 (y i -x ⊤ i,γ β γ,K ) 2 = (ỹ n+1 -x⊤ n+1,γ β γ,K ) 2 ≤ (r + x max K) 2 .Similarly, we have u(γ; X, ỹ) -u(γ; X, y) ≤ (r + x max K) 2 . Next, the (ε, 0)-DP follows from Lemma 2.2. This finishes the proof.Consider the notation in Section B.1 and recall the event E K := ∩ γ:γ∈As {β γ,K = β γ,ols }. We will For notational brevity, we use L γ to denote L γ (X, y). Now, we restrict ourselves to the event E K . therefore we have L γ,K = L γ for all γ. To establish a lower bound π(γ * ), we make use of its specific form, thereby obtaining the following inequality:π(γ * ) = 1 1 + γ ′ ∈As exp - ε(L γ ′ -L γ * ) ∆u . Now we fix k ∈ [s],and and consider any γ ∈ A s,k . For any η ∈ [0, 1], note thatn -1 (L γ -L γ * ) = n -1 {y ⊤ (I n -Φ γ )y -y ⊤ (I n -Φ γ * )y} = n -1 (X γ * \γ β γ * \γ + w) ⊤ (I n -Φ γ )(X γ * \γ β γ * \γ + w) -w ⊤ (I n -Φ γ * )w = ηβ ⊤ γ * \γ Γ(γ)β γ * \γ + 2 -1 (1 -η)β ⊤ γ * \γ Γ(γ)β γ * \γ -2 n -1 (I n -Φ γ )X γ * \γ β γ * \γ ⊤ (-w) + 2 -1 (1 -η)β ⊤ γ * \γ Γ(γ)β γ * \γ -n -1 w ⊤ (Φ γ -Φ γ * )w.Consider the random variable Following the analysis of Theorem 2.1 in [17] , we haveP max γ∈A s,k 2n -1 {(I n -Φ γ )X γ * \γ β γ * \γ } ⊤ w ≥ 2 -1 (1 -η)β ⊤ γ * \γ Γ(γ)β γ * \γ ≤ 2e -6k log p ,and,P max γ∈A s,k n -1 w ⊤ (Φ γ -Φ γ * )w ≥ 2 -1 (1 -η)β ⊤ γ * \γ Γ(γ)β γ * \γ ≤ 4e -2k log p , whenever min γ∈A s,k β ⊤ γ * \γ Γ(γ)β γ * \γ k ≥ Cσ 2 log p n(1 -η)for large enough universal constant C > 0. Setting η = 1/2, we note that whenever m * (s) ≥ 2Cσ 2 {(log p)/n}, we getn -1 (L γ -L γ * ) ≥ 1 2 β ⊤ γ * \γ Γ(γ)β γ * \γ ≥ km * (s) 2 for all γ ∈ A s ,with probability at least 1 -2p -6 -4p -2 . Also, note that β ⊤ γ * \γ Γ(γ)β γ * \γ ≤ κ + sb 2 max . Hence, if we havem * (s) ≥ max 2C, 16∆u εσ 2 σ 2 log p n ,the following are true:γ ′ ∈As exp - ε(L γ ′ -L γ * ) ∆u ≤ γ ′ ∈As exp - nkεm * (s) 2∆u ≤ s k=1 p -s k s k exp - nkεm * (s) 2∆u ≤ s k=1 p 2k .p -4k ≤ p -2 .Therefore, we havemin γ∈As∪{γ * } π(γ) ≥ 1 1 + p -2 ≥ 1 -p -2with probability 1 -2p -6 -4p -2 . Now by the discussion in Section B.1, we haveEQUATION-)σx max . This finishes the proof.For clarity, we first specify some notations. Let γ D t denote the model update of MH chain run over dataset D. Let τ D η be the corresponding η-mixing time. Let D and D ′ be two neighboring datasets, and π D and π D ′ be the corresponding probability mass functions for the exponential mechanism. Then, we have the following:P(γ D τ D η = γ) ≤ π D (γ) + η ≤ e ε π D ′ (γ) + η ≤ e ε P(γ D ′ τ D ′ η = γ) + η(1 + e ε ).This finishes the proof.We again restrict ourselves to the event E K with K ≥ √ s ( κ+ κ-)b max + ( 8 κ-)σx max . For the proof, let P denote the transition matrix of the original Metropolis-Hastings sampler (10) . In this case, the state space is S = A s ∪ {γ * }. Now consider the transition matrix P = P/2 + I n /2, corresponding to the lazy version of the random walk that stays in its current position with a probability of at least 1/2. Due to the construction, the smallest eigenvalue of P is always non-negative, and the mixing time of the chain is completely determined by the second largest eigenvalue λ 2 of P. To this end, we define the spectral gap Gap(P) = 1 -λ 2 , and for any lazy Markov chain, we have the following sandwich relation [41, 50] EQUATION)Lower Bound on π(•) :To establish a lower bound on the target distribution in (7) , we make use of its specific form, thereby obtaining the following inequality:π(γ) = π(γ * ). π(γ) π(γ * ) = 1 1 + γ ′ ∈As exp - ε(L γ ′ -L γ * ) ∆u . exp - ε(L γ -L γ * ) ∆u .γ ∈ A s,k .Similar to the proof of Section A.3, we note that whenever m * (s) ≥ 2Cσ 2 {(log p)/n} for a large enough universal constant C > 0, we get3 2 β ⊤ γ * \γ Γ(γ)β γ * \γ ≥ n -1 (L γ -L γ * ) ≥ 1 2 β ⊤ γ * \γ Γ(γ)β γ * \γ ≥ km * (s) 2 for all γ ∈ A s , with probability at least 1 -2p -6 -4p -2 . Also, note that β ⊤ γ * \γ Γ(γ)β γ * \γ ≤ κ + sb 2 max . Hence, if we have m * (s) ≥ max 2C, 16∆u εσ 2 σ 2 log p n ,the following are true:γ ′ ∈As exp - ε(L γ ′ -L γ * ) ∆u ≤ γ ′ ∈As exp - nkεm * (s) 2∆u ≤ s k=1 p -s k s k exp - nkεm * (s) 2∆u ≤ s k=1 p 2k .p -4k ≤ p -2 ,and,exp - ε(L γ -L γ * ) ∆u ≥ exp - 3nεβ ⊤ γ * \γ Γ(γ)β γ * \γ 2∆u ≥ exp - 3nsεκ + b 2 max 2∆uCombining these two facts we haveEQUATIONwith probability 1 -2p -6 -4p -2 .Now it remains to prove a lower bound on the spectral gap Gap(P), and we do so via the canonical path argument [41] . We begin by describing the idea of a canonical path ensemble associated with a Markov chain. Given a Markov chain C with state space S , consider the weighted directed graph G(C) = (V, E) with vertex set V = S and the edge set E in which a ordered pair e = (γ, γ ′ ) is included as an edge with weight Q(e) = Q(γ, γ ′ ) = π(γ)P(γ, γ ′ ) iff P(γ, γ ′ ) > 0. A canonical path ensemble T corresponding to C is a collection of paths that contains, for each ordered pair (γ, γ ′ ) of distinct vertices, a unique simple path T γ,γ ′ connecting γ and γ ′ . We refer to any path in the ensemble T as a canonical path. [41] shows that for any reversible Markov chain and nay choice of a canonical path ensemble T , the spectral gap of P is lower bounded asEQUATIONwhere ℓ(T ) corresponds to the length of the longest path in the ensemble T , and the quantity ρ(T) := max e∈E 1 Q(e) (γ,γ ′ ):e∈T γ,γ ′ π(γ)π(γ ′) is known as the path congestion parameter. Thus, it boils down to the construction of a suitable canonical path ensemble T . Before going into further details, we introduce some working notations. For any two given paths T 1 and T 2 :• Their intersection T 1 ∩ T 2 denotes the collection of overlapping edges.• If T 2 ⊂ T 1 , then T 1 \ T 2 denotes the path obtained by removing all the edges of T 2 from T 1 .• We use T1 to denote the reverse of T 1 .• If the endpoint of T 1 is same as the starting point of T 2 , then T 1 ∪ T 2 denotes the path obtained by joining T 1 and T 2 at that point.We will now shift focus toward the construction of the canonical path ensemble. At a high level, our construction follows the same scheme as in [51] .First, we need to construct the canonical path T γ,γ * from any γ ∈ S to the true model γ * . To this end, we introduce the concept of memoryless paths. We call a set T M of canonical paths memoryless with respect to the central state γ * if 1. for any state γ ∈ S satisfying γ ̸ = γ * , there exists a unique simple path T γ,γ * in T M connecting γ and γ * ;2. for any intermediate state γ ∈ S on any path T γ,γ * ∈ T M , the unique path connecting γ and γ * is the sub-path of T γ,γ * starting from γ and ending at γ * .Intuitively, this memoryless property tells that for any intermediate step in any canonical path, the next step towards the central state does not depend on history. Specifically, the memoryless canonical path ensemble has the property that in order to specify the canonical path connecting any state γ ∈ S and the central state γ * , we only need to specify the next state from γ ∈ S \ {γ * }, i.e., we need a transition function G : S \ {γ * } → S that maps the current state γ to the next state. For simplicity, we define G(γ * ) = γ * to make S as the domain of G. For a more detailed discussion, we point the readers to Section 4 of [51] . We now state a useful lemma that is pivotal to the construction of the canonical path ensemble.Lemma A.1 ([51]). If a function G : S \ {γ * } → S satisfies the condition d H (G(γ), γ * ) < d H (γ, γ * ) for any state γ ∈ S \ {γ * }, then G is a valid transition map.Using the above lemma, we will now construct the memoryless set of canonical paths from any state γ ∈ S to γ * by explicitly specifying a transition map G. In particular, we consider the following transition function:• If γ ̸ = γ * , we define G(γ) to be γ ′ , whch is formed by replacing the least influential covariate in γ with most influential covariate in γ * \ γ. In notations, we have γ ′ j = γ j for all j / ∈ {j γ , k γ }, γ ′ jγ = 1 and γ ′ kγ = 0, where Based on this, we are now ready to construct the canonical path ensemble T . Specifically, due to memoryless property, two simple paths T γ,γ * and T γ ′ ,γ * share an identical subpath to γ * starting from their first common intermediate state. Let T γ∩γ ′ denote the common sub-path T γ∩γ * ∩ T γ ′ ∩γ * , and T γ\γ ′ := T γ,γ * \ T γ∩γ ′ denotes the remaining path of T γ,γ * after removing the segment T γ∩γ ′ . The path T γ ′ \γ is defined in a similar way. Then it follows that T γ\γ ′ and T γ ′ \γ have the same endpoint. Therefore, it is allowed to consider the path T γ\γ ′ ∪ Tγ ′ \γ .j γ := arg max j∈γ * \γ Φ γ∪{j} X γ * β γ * 2 2 and k γ := arg min k∈γ\γ * Φ γ∪{j} X γ * β γ * 2 2 -Φ γ∪{j}\{k} X γ * β γ *We call γ a precedent of γ ′ if γ ′ is on the canonical path T γ,γ * ∈ T , and a pair of states γ, γ ′ are adjacent if the canonical path T γ,γ ′ is e γ,γ ′ , the edge connecting γ and γ ′ . Next, for γ ∈ S , defineEQUATION)denote the set of all precedents. We denote by |T | the length of the path T . The following lemma provides some important properties of the previously constructed canonical path ensemble. Lemma A.2. For any distinct pair (γ, γ ′ ) ∈ S × S :EQUATION|T γ,γ ′ | ≤ 1 2 {d H (γ, γ * ) + d H (γ ′ , γ * )} ≤ 2s.(b) If γ and γ ′ are adjacent and γ is precedent of of γ ′ , then{(γ, γ′ ) | e γ,γ ′ ∈ T γ,γ ′ } ⊂ Λ(γ) × S .Proof. For the first claim, let us first assume that |T γ,γ * | = k, i.e., G k (γ) = γ * for the appropriate transition map G. Also, recall that |γ| = |γ * | = s. Hence, due to an elementary iterative argument, it follows that2s ≥ d H (γ, γ * ) = d H (G(γ), γ * ) + 2 = d H (G 2 (γ), γ * ) + 4 . . . = 2k. Also, note that |T γ,γ ′ | ≤ |T γ,γ * | + |T γ ′ ,γ * |.Hence, the claim follows using the previous inequality.For the second claim, note that for any pair (γ, γ′ ) such that T γ,γ ′ ∋ e γ,γ ′ , we have two possible options : (i) e γ,γ ′ ∈ T γ\γ ′ , or (ii) e γ,γ ′ ∈ T γ′ \γ . As γ is precedent of γ ′ , the only possibility that we have is e γ,γ ′ ∈ T γ\γ ′ . This shows that γ belongs to the path T γ,γ * and γ ∈ Λ(γ).According to Lemma A.2(b), the path congestion parameter ρ(T ) satisfiesEQUATION)where the setΓ * := {(γ, γ ′ ) ∈ S × S | T γ,γ ′ = e γ,γ ′ , γ ∈ Λ(γ ′ )}.Here we used the fact that the weight function Q satisfies the reversibility condition Q(γ, γ ′ ) = Q(γ ′ , γ) in order to restrict the range of the maximum to pairs (γ, γ ′ ) where γ ∈ Λ(γ ′ ).For the lazy form of the Metropolis-Hastings walk (10) , we haveQ(γ, γ ′ ) = π(γ)P(γ, γ ′ ) ≥ 1 ps π(γ) min 1, π(γ ′ ) π(γ) ≥ 1 ps min {π(γ), π(γ ′ )} .Substituting this bound in (19) , we getEQUATIONIn order to prove that ρ(T ) = O(ps) with high probability, it is sufficient to prove that the two terms inside the maximum are O(1). To this end, we introduce two useful lemmas. Lemma A.3. Consider the eventA n = max γ∈S ,ℓ / ∈γ w ⊤ (Φ γ∪{ℓ} -Φ γ )w ≤ 12σ 2 s log p Then we have P(A n ) ≥ 1 -p -2 .Proof. First note that w ⊤ (Φ γ∪{ℓ} -Φ γ )w = (h ⊤ γ,ℓ w) 2 for an appropriate unit vector h γ,ℓ depending only upon X γ and X ℓ . By Sub-gaussian tail inequality, we haveP (h ⊤ γ,ℓ w) 2 ≥ t ≤ 2e -t 2σ 2 .Setting t = 12σ 2 s log p and applying an union bound we getP max γ∈S ,ℓ / ∈γ (h ⊤ γ,ℓ w) 2 ≥ 12σ 2 s log p ≤ 2 p s (p -s)p -6s ≤ 2p -3s ≤ p -2 .Lemma A.4. Suppose that, in addition to the conditions in Theorem 4.3, the event A n holds. Then for all γ ̸ = γ * , we haveπ(γ) π(G(γ)) ≤ p -3 .Moreover, for all γ, π[Λ(γ)] π(γ) ≤ 2.Therefore, both Lemma A.3 and Lemma A.4 give ρ(T ) ≤ 2ps with probability 1 -p -2 . Lemma A.2(a) suggests that ℓ(T ) ≤ 2s. Therefore, Equation (17) shows that Gap(P) ≥ 1 4ps 2 with probability 1 -p -2 . Finally, combining ( 16) and ( 15), we get the following with 1 -8p -2τ η ≤ C 2 ps 2    nεκ + b 2 max r + ( κ+ κ-)b max x max + ( σ κ-)x 2 max 2 + log(1/η)    ,where C 2 > 0 is a universal constant. Finally, the proof is concluded by arguing thatP(E c K ) ≤ 2p -2 .Now we are ready to bound β γ,K 1 . Define the OLS estimator corresponding to the model γ asβ γ,ols = ( X ⊤ γ X γ n ) -1 X ⊤ γ X γ * β γ * n :=u1 + ( X ⊤ γ X γ n ) -1 X ⊤ γ w n :=u2.In this section, we will show that there exists a choice for K such that the event E K := ∩ γ:|γ|=s {β γ,ols = β γ,K } holds with high probability. By Holder's inequality we have∥u 1 ∥ 2 ≤ (X ⊤ γ X γ /n) -1 op X ⊤ γ X γ * /n op β γ * 2 ≤ ( κ+ κ-)b max .Hence, an application of Cauchy-Schwarz inequality directly yields that ∥u 1 ∥ 1 ≤ 2( κ+ κ-)√ sb max . Next, note that ∥u 2 ∥ 2 ≤ ( X ⊤ γ X γ n ) -1 2 X ⊤ γ w n 2 ≤ √ s ( X ⊤ γ X γ n ) -1 2 X ⊤ γ w n ∞ ≤ √ s ( X ⊤ γ X γ n ) -1 2 X ⊤ w n ∞ .Therefore, we get ∥u 2 ∥ 1 ≤ s κ-X ⊤ w/n ∞ . In order to upper bound the last term in the previous inequality, we define D i,j = X[i, j]w j for all (i, j) ∈ [s] × [n]. Using the sub-Gaussian property of w j , we have E(e λwj ) ≤ e λ 2 x 2 max σ 2 /2 . Therefore, due to Hoeffding's inequality, we haveP   1 n j∈[n] D i,j ≥ 8σx max log p n   ≤ 2p -4 . Note that X ⊤ w/n ∞ = max i∈[s] n -1 j∈[n] D i,j. Hence, by simple union-bound argument, it follows thatP max γ:|γ|=s X ⊤ γ w n ∞ ≥ 8σx max log p n ≤ 2p -4 ≤ 2p -2 .Thus, Assumption 3.4(c) yields that β γ,K2 1 ≤ s ( κ+ κ-)b max + ( 8 κ-)σx max 2 . Therefore, if K ≥ √ s ( κ+ κ-)b max + ( 8 κ-)σx max then P(E K ) ≥ 1 -2p -2 .Based on Theorem 4.3, we have ∥π t -π∥ TV ≤ eta with probability at least 1 -c 2 p -2 whenever t is sufficiently large. Also, by Theorem 3.5, we know π(γ * ) ≥ 1 -p -2 with probability 1 -c 1 p -2 . Therefore, we have π t (γ * ) ≥ 1 -η -p -2 with probability at least 1 -(c ! + c 2 )p -2 . This finishes the proof.Let j γ , k γ be the indices defined in the construction of G(γ). The we have γEQUATIONBy the form in (7) , we haveπ(γ) π(γ ′ ) = exp - y ⊤ (Φ γ ′ -Φ γ )y (∆u/ε) .To show that the above ration is O(1), it suffices to show that y ⊤ (Φ γ ′ -Φ γ )y is large. By simple algebra, it follows thaty ⊤ (Φ γ ′ -Φ γ )y = y ⊤ (Φ γ∪{jγ } -Φ γ )y -y ⊤ (Φ γ∪{jγ } -Φ γ ′ )y = ∥v 1 ∥ 2 2 + 2v ⊤ 1 w + w ⊤ (Φ γ∪{jγ } -Φ γ )w -∥v 2 ∥ 2 2 + 2v ⊤ 2 w + w ⊤ (Φ γ∪{jγ } -Φ γ ′ )w = ∥v 1 ∥ 2 2 + 2v ⊤ 1 (Φ γ∪{jγ } -Φ γ )w + w ⊤ (Φ γ∪{jγ } -Φ γ )w -∥v 2 ∥ 2 2 + 2v ⊤ 2 (Φ γ∪{jγ } -Φ γ ′ )w + w ⊤ (Φ γ∪{jγ } -Φ γ ′ )w ≥ ∥v 1 ∥ 2 (∥v 1 ∥ 2 -2 (Φ γ∪{jγ } -Φ γ )w 2 ) -∥v 2 ∥ 2 (∥v 2 ∥ 2 + 2 (Φ γ∪{jγ } -Φ γ ′ )w 2 ) -(Φ γ∪{jγ } -Φ γ ′ )w 2 2 .(22) Now, we recall the eventA n = max γ∈S ,ℓ / ∈γ w ⊤ (Φ γ∪{ℓ} -Φ γ )w ≤ 12σ 2 s log p . Let A 2 := nκ -m * (s) ≥ κ -C 0 σ 2 log p.Then for C 0 large enough so that κ -C 0 ≥ (128 × 12)s, Equation ( 22) leads to the following inequality under event A n :EQUATIONunder the margin condition of Theorem 4.3.From the previous part, the bound (23) implies that π(γ)/π(G(γ)) ≤ p -3 . For each γ ∈ Λ(γ), we have that γ ∈ T γ,γ ⊂ T γ,γ * . Let the path T γ,γ be γ 0 → γ 1 → . . . → γ k , where k = |T γ,γ| is the length of the path, and γ 0 = γ and γ k = γ are the two endpoints. Now note that {γ ℓ } ℓ≤k-1 ⊂ S , and (23) ensures that π(γ) π(γ) = k ℓ=1 π(γ ℓ-1 ) π(γ ℓ ) ≤ p -3k .Also, by Lemma A.2(a) we have k ∈ [s]. Now, we count the total number of sets in Λ(γ) for each k ∈ [s]. Recall that by the construction of the canonical path, we update the current state by adding a new influential covariate and deleting one unimportant one. Hence any state in S has at most sp adjacent precedents, implying that there could be at most s k p k distinct paths of length k. This entails that π(Λ(γ)) π(γ)≤ γ∈Λ(γ) π(γ) π(γ) ≤ s k=1 (ps) k p -3k ≤ s k=1 p -k ≤ 1 1 -1/p ≤ 2.Recall the definition of j γ and k γ . The first result in the following lemma shows that the gain in adding j γ to the current model γ is at least nκ -m * (s). The second result shows that the loss incurred by removing k γ from the model γ ∪ {j γ } is at most nκ -m * (s)/2. As a result, it follows that it is favorable to replace X kγ with the more influential feature X jγ in the current model γ.Lemma C.1. Under Assumption 3.4(b) and Assumption 4.2, the following hold for all γ ∈ A s :EQUATION(b) Φ γ∪{jγ } X γ * β γ * 2 2 -Φ γ∪{jγ }\{k} X γ * β γ * 2 2 ≤ nκ -m * (s)/2.Proof. For each ℓ ∈ γ * \ γ, we haveΦ γ∪{ℓ} X γ * β γ * 2 2 -Φ γ X γ * β γ * 2 2 = β ⊤ γ * X ⊤ γ * (Φ γ∪{ℓ} -Φ γ )X γ * β γ * = β ⊤ γ * X ⊤ γ * (I n -Φ γ )X ℓ X ⊤ ℓ (I n -Φ γ )X γ * β γ * X ⊤ ℓ (I n -Φ γ )X ℓ ≥ β ⊤ γ * \γ X ⊤ γ * \γ (I n -Φ γ )X ℓ X ⊤ ℓ (I n -Φ γ )X γ * \γ β γ * \γ n ,where the second equality simply follows from Gram-Schmidt orthogonal decomposition. By summing the preceding inequality over ℓ ∈ γ * \ γ, we getℓ∈γ * \γ Φ γ∪{ℓ} X γ * β γ * 2 2 -Φ γ X γ * β γ * 2 2 ≥ β ⊤ γ * \γ X ⊤ γ * \γ (I n -Φ γ )X γ * \γ X ⊤ γ * \γ (I n -Φ γ )X γ * \γ β γ * \γ n ≥ κ -β ⊤ γ * \γ X ⊤ γ * \γ (I n -Φ γ )X γ * \γ β γ * \γ ≥ nκ -|γ \ γ * | m * (s) = nκ -|γ * \ γ| m * (s).The last inequality follows from the fact that |γ| = |γ * | = s. Since j γ maximizes Φ γ∪{ℓ} X γ * β γ * 2 2 over all ℓ ∈ γ * \ γ, the preceding inequality implies thatΦ γ∪{jγ } X γ * β γ * 2 2 -Φ γ X γ * β γ * 2 2 ≥ nκ -m * (s).Similarly, to prove the second claim, first note that for any k ∈ γ \ γ * , we haveEQUATION2 over all possible k ∈ γ \ γ * , by Assumption 4.2 we haveΦ γ∪{jγ } X γ * β γ * 2 2 -Φ γ∪{jγ }\{k} X γ * β γ * 2 2 ≤ nκ -m * (s)/2.Under the setup of Section 5, we consider the privacy parameter ε ∈ {0.5, 1, 3, 5, 10}. For the Metropolis-Hastings random walk, we vary K ∈ {0.5, 2, 3, 3.5} and initialize 10 independent Markov chains from random initializations and record the F-score of the last iteration. We also track the qualities of the model through its explanatory power. In particular, we calculate the scale factor R γ := y ⊤ Φ γ y/ ∥y∥ Strong signal: Under this setup, note that the model estimate of ABESS exactly matches the true model. For ε ≥ 3 and K ≥ 2, Figure 1 shows that all the chains have identified a reasonably good estimate of the true model γ * within 50p iterations. This empirical phenomenon validates theoretical findings in Theorem 4.3. However, for larger values of K the performance is worse as the noise level is also large. On the other hand, for the case of K = 0.5, the performance is also worse due to too much shrinkage that results in a bad estimate of β. The mean F-score's also suggest the same phenomenon. For smaller values of ε, the performance is generally bad due to increased noise level. This is expected as higher privacy usually entails a worse performance in terms of utility.Weak signal: We perform the same experiments under a weak signal regime. As expected, both Figure 2 and Table 1 show that the performance of the proposed algorithm is generally inferior to that in the strong signal regime for K ≥ 2. However, note that our algorithm enjoys a better utility for K = 0.5. In fact, performance is as good as the non-private BSS for ε ≥ 3. This is not surprising as K = 0.5 closer to ∥β∥ 1 ≈ 0.7 in the weak signal case and results in better estimation for β. On the contrary, larger values of K inject more noise into the algorithm and the utility deteriorates.We consider an independent Gaussian design matrix, formed by sampling entries from identical independent standard normal distributions and normalized by the ℓ ∞ norm. Specifically, we set n = 900, p = 2000 with the sparsity level s = 4. Similar to the setup in Section 5, we generate entries with independent Uniform(-0.1, 0.1) noise w following the linear model (1) . We choose the design vector β with true sparsity s = 4 and the support set γ * = {j : 1 ≤ j ≤ 4}. All the signal strengths are set to be equal, taking the following two forms: (i) Strong signal: β j = 2{(s log p)/n} 1/2 , and (ii) Weak signal: β j = 2{(log p)/n} 1/2 for all j ∈ γ * . We consider the privacy parameter ε ∈ {0.5, 1, 3, 5, 10}. For the Metropolis-Hastings random walk, we vary K ∈ {0.5, 2, 3, 3.5} and initialize 10 independent Markov chains from random initialization and record the F-score of the last iteration.Strong signal: Under this setup, note that the model estimate of ABESS exactly matches the true model with F-score = 1. For the case of K = 0.5, Figure 3 shows the performance is better compared with settings K ≥ 2 when ε ≤ 5. However, when ε = 10, the performance is worse due to shrinkage of the estimate of β while the estimations in other settings are easier due to lower privacy requirements. For higher values of ε, the performance is generally strong because of the reduced noise level. This is expected since lower privacy typically leads to better utility performance. Notice that for ε = 10 and K = 2, we have a fairly accurate estimate of the true model γ * within 50p iterations.Weak signal: We conduct the same experiments under a weak signal regime. As expected, Figure 4 shows that the performance of the proposed algorithm is generally inferior to that in the strong signal regime for K ≥ 2. However, note that our algorithm enjoys a better utility for K = 0.5 when ε ≥ 3. In fact, performance is as good as the non-private BSS for ε = 10. This is not surprising as K = 0.5 closer to ∥β∥ 1 ≈ 0.7 in the weak signal case, leading to better estimation for β. In contrast, larger values of K introduce more noise into the algorithm and weaken the utility.• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 


Decision: The paper addresses the model selection problem in sparse linear regression, proposing a differentially private version of the best subset selection algorithm using the exponential mechanism. The authors show that their algorithm has a better privacy-utility (in terms of minimum value of coefficients required for identifiability) tradeoff than existing approaches. Additionally, they demonstrate that a simple MCMC algorithm for the problem converges efficiently under an additional condition (intuitively saying that spurious features have low correlation with true features). The main issue with the work is the required condition cannot be verified (without knowledge of true features) and privacy guarantees of the algorithm depend on it (since they depend on convergence). In addition, the empirical evaluation does not compare with any baselines. Overall, I recommend acceptance with some reservations.

Reviews: review: paper_summary: paper_summary: This paper studies the best subset selection (BSS) problem in high-dimensional sparse linear regression. The results of this paper are roughly as follows:
- First, adopt the exponential mechanism to design a DP BSS algorithm. The statistical/privacy guarantee of this approach can be derived based on standard techniques.
- Second, devise an efficient MCMC sampling algorithm that still enjoys approximate DP, while having a polynomial mixing time.
- Numerical experiments focus on a synthetic example, with a random design matrix, as well as a US census study and another socio-economic study.
Finally, the paper points out several open directions that would be plausible for future work.
strength_weakness: strength_weakness: - The simulation set up is somewhat restricted. Also, the main takeaways are not clear.

- The identifiablity margin assumption is not clearly explained; the intuition as to while this should be expected is not clear.

- In section 4, a comparison from the technical aspects is currently missing, and needs to be more clear to distinguish this work from prior works.
questions: questions: - What is the main technical innovation of this work as compared to prior work?
- I understand that this paper is focused on BSS in high-dimensional linear regression. What are the practical takeaways of the results beyond this setup? Could you discuss it in this context?
limitations: limitations: The authors discuss limitations of their work in Section 6.
paper_summary: paper_summary: This papers study the differentially private best subset selection (BSS) of sparse linear model selection.
strength_weakness: strength_weakness: - The comparisons with existing works are not comprehensive (see the comments below). 

- The assumptions on both the sensitivity and margin need a better justification. 

- The algorithm is based on computationally-intensive exponential sampling and it is not clear how practical the proposed MCMC improved version is. 

- Code is not released and variation of the experimental results are marked.
questions: questions: 1. As partially mentioned by the authors in the introduction, there are many existing more computationally-efficient sparse optimizations. So compared to privately select a subset of coordinates in the first place, one may also apply L_1-norm restriction or applying sparse vector techniques as alternative solutions.  I am wondering whether the authors could give a more comprehensive studies on those existing works with a more clear motivation on the advantage of BSS. 

2. Is it possible to remove the data-dependent assumption for bounded sensitivity (Assumption 3.2) but instead through clipping or truncation give an data-independent sensitivity bound? 

3. I think the MCMC of exponential sampling could be an independent and possibly more important contribution of this paper, which may provide an efficient version of exponential sampling methods. I am wondering whether there is any prior work on this topic and more rigorously how to determine the constant in Theorem 4.3 and Corollary 4.4.
limitations: limitations: Please see the above comments.
paper_summary: paper_summary: The paper addresses the model selection problem in sparse linear regression, proposing a differentially private version of the best subset selection algorithm using the exponential mechanism. The authors prove that the proposed private algorithm requires $O(\sigma^2 s\log p/n\epsilon)$ samples to identify the correct features with high probability. Additionally, they demonstrate that when spurious features have low correlation with true features, sampling through MCMC converges rapidly to the stationary distribution.
strength_weakness: strength_weakness: The rapid convergence of MCMC is achievable only when spurious features have low correlation with true features. However, practitioners typically do not know which features are true in advance, making it difficult to assess the strength of this correlation. Therefore, assumption 4.2 cannot be tested in practice.
questions: questions: Does the universal constant $C_1$ in theorem 3.5 depend on the exponent of the failure probability?
limitations: limitations: See weaknesses.
paper_summary: paper_summary: This paper addresses the challenge of model selection in high-dimensional sparse linear regression models under privacy constraints. It proposes a differentially private algorithm for best subset selection using the exponential mechanism, providing strong statistical utility guarantees under high-privacy regimes. The authors demonstrate that accurate model recovery is possible in low-privacy regimes, achieving minimax optimal conditions similar to non-private settings. To tackle computational complexity, the paper introduces a Metropolis-Hastings algorithm with polynomial mixing time under certain conditions, ensuring approximate differential privacy and good statistical utility. Experiments on simulated data illustrate the algorithm's effectiveness in identifying active features under reasonable privacy budgets.
strength_weakness: strength_weakness: 1. The high-probability events in Theorem 3.5, Theorem 4.3, and Corollary 4.4 could be further clarified. See question 1.
questions: questions: 1. Theorem 3.5, Theorem 4.3, and Corollary 4.4 are high-probability statements. Is it possible to simplify them by including the high probability (success probability) $1 - c_i p^{-2}$ into the conclusion? For example, can the failure probability $1 - c_2 p^{-2}$ be included in the total variation distance to the stationary distribution in the mixing time?
2. In the proof of the utility guarantee (Theorem 3.5), the authors cite the analysis of Guo et al. (2020) titled "Best subset selection is robust against design dependence." Could the authors cite the specific theorems?
3. What is the disadvantage of output perturbation to address the DP BSS problem compared to the exponential mechanism (the authors propose)? Can output perturbation achieve pure DP with a good utility guarantee?
limitations: limitations: Please check the questions above.

Comments: comment: Thank you for your work on the rebuttal. I will keep my rating.