Derivatives and residual distribution of regularized M-estimators with application to adaptive tuning
Year: 2021

Abstract:
This paper studies M-estimators with gradient-Lipschitz loss function regularized 1 with convex penalty in linear models with Gaussian design matrix and arbitrary 2 noise distribution. A practical example is the robust M-estimator constructed with 3 the Huber loss and the Elastic-Net penalty and the noise distribution has heavy-tails. 4 Our main contributions are three-fold. (i) We provide general formulae for the 5 derivatives of regularized M-estimators β̂(y,X) where differentiation is taken with 6 respect to both y andX; this reveals a simple differentiability structure shared by 7 all convex regularized M-estimators. (ii) Using these derivatives, we characterize 8 the distribution of the residual ri = yi−xi β̂ in the intermediate high-dimensional 9 regime where dimension and sample size are of the same order. (iii) Motivated 10 by the distribution of the residuals, we propose a novel adaptive criterion to select 11 tuning parameters of regularized M-estimators. The criterion approximates the 12 out-of-sample error up to an additive constant independent of the estimator, so 13 that minimizing the criterion provides a proxy for minimizing the out-of-sample 14 error. The proposed adaptive criterion does not require the knowledge of the 15 noise distribution or of the covariance of the design. Simulated data confirms the 16 theoretical findings, regarding both the distribution of the residuals and the success 17 of the criterion as a proxy of the out-of-sample error. Finally our results reveal 18 new relationships between the derivatives of β̂(y,X) and the effective degrees of 19 freedom of the M-estimator, which are of independent interest. 20

text: This paper studies M-estimators with gradient-Lipschitz loss function regularized1 with convex penalty in linear models with Gaussian design matrix and arbitrary2 noise distribution. A practical example is the robust M-estimator constructed with3 the Huber loss and the Elastic-Net penalty and the noise distribution has heavy-tails.4 Our main contributions are three-fold. (i) We provide general formulae for the5 derivatives of regularized M-estimators β̂(y,X) where differentiation is taken with6 respect to both y andX; this reveals a simple differentiability structure shared by7 all convex regularized M-estimators. (ii) Using these derivatives, we characterize8 the distribution of the residual ri = yi−x>i β̂ in the intermediate high-dimensional9 regime where dimension and sample size are of the same order. (iii) Motivated10 by the distribution of the residuals, we propose a novel adaptive criterion to select11 tuning parameters of regularized M-estimators. The criterion approximates the12 out-of-sample error up to an additive constant independent of the estimator, so13 that minimizing the criterion provides a proxy for minimizing the out-of-sample14 error. The proposed adaptive criterion does not require the knowledge of the15 noise distribution or of the covariance of the design. Simulated data confirms the16 theoretical findings, regarding both the distribution of the residuals and the success17 of the criterion as a proxy of the out-of-sample error. Finally our results reveal18 new relationships between the derivatives of β̂(y,X) and the effective degrees of19 freedom of the M-estimator, which are of independent interest.20
1 Introduction21
This paper studies properties of robust estimators in linear models y = Xβ∗ + ε with response22 y ∈ Rn, unknown regression vector β∗ whereX is a design matrix with n rows x1, ...,xn, each row23 xi being a high-dimensional feature vector in Rp with covariance Σ. Throughout, let β̂ = β̂(y,X)24 be a regularized M -estimator given as a solution of the convex minimization problem25
β̂(y,X) = argminb∈Rp 1 n ∑n i=1 ρ(yi − x>i b) + g(b) (1)
where ρ : R → R is a convex data-fitting loss function and g : Rp → R a convex penalty. We26 may write β̂ρ,g(y,X) for (1) to emphasize the dependence on the loss-penalty pair (ρ, g); if the27 argument (y,X) is dropped then β̂ is implicitly understood at the observed that (y,X). Typical28 examples of losses include the square loss ρ(u) = u2/2, the Huber loss H(u) = ∫ |u| 0 min(1, t)dt29 or its scaled version ρ = Λ2H(u/Λ) for some tuning parameter Λ > 0, while typical examples of30 penalty functions include the Elastic-Net g(b) = λ‖b‖1 + µ‖b‖2/2 for tuning parameters λ, µ ≥ 0.31 The paper introduces the following criterion to select a loss-penalty pair (ρ, g) with small out-of-32 sample error ‖Σ1/2(β̂ − β∗)‖2: for a given set of candidate loss-penalty pairs {(ρ, g)} and the33
Submitted to 35th Conference on Neural Information Processing Systems (NeurIPS 2021). Do not distribute.
corresponding M -estimator β̂ρ,g in (1), select the pair (ρ, g) that minimizes the criterion34
Crit(ρ, g) = ∥∥∥r+ d̂f tr[V ] ψ ( r )∥∥∥ 2 with    r = y −Xβ̂ρ,g ∈ Rn, d̂f = tr[X(∂/∂y)β̂ρ,g] ∈ R, V = diag{ψ′(r)}(In −X(∂/∂y)β̂ρ,g) ∈ Rn×n (2)
where tr[·] is the trace, ψ : R → R is the derivative of ρ, ψ′ the derivative of ψ and we extend ψ35 and ψ′ to functions Rn → Rn by componentwise application of the univariate function of the same36 symbol. Above, (∂/∂y)β̂ρ,g ∈ Rp×n denotes the Jacobian of (1) with respect to y for X fixed,37 at the observed data (y,X). As we will see while studying particular examples, for pairs (ρ, g)38 commonly used in robust high-dimensional statistics such as the square loss, Huber loss with the39 `1-penalty or Elastic-Net penalty, the ratio d̂f/ tr[V ] in (2) admits simple, closed-form expressions40 and can be computed at a negligible computational cost once β̂ρ,g(y,X) itself has been computed.41 The criterion (2) has an appealing adaptivity property: it does not require any knowledge of the noise42 ε or its distribution, nor any knowledge of the covariance Σ of the design.43
1.1 Contributions44
1. The end goal of paper is to provide theoretical justification and theoretical guarantees for the45 criterion (2) in the high-dimensional regime where the ratio p/n has a finite limit and X has46 anisotropic Gaussian distribution. The theoretical results will justify the approximation47
∥∥r + ( d̂f/ tr[V ] ) ψ ( r )∥∥2/n ≈ ‖ε‖2/n+ ‖Σ1/2(β̂ − β∗)‖2. (3)
Figure 1 illustrates the accuracy of (3) on simulated data. To study the criterion (2) and derive the48 approximation (3), we develop novel results of independent interest regarding M -estimators in (1):49
2. The paper derives general formula for the derivatives (∂/∂yi)β̂ and (∂/∂xij)β̂. This sheds light50 on the differentiability structure ofM -estimators for general loss-penalty pairs: for any ρ, g with g51 strongly convex, there exists Â ∈ Rp×p depending on (y,X) such that for almost every (y,X),52
(∂/∂yi)β̂(y,X) = ÂX >eiψ ′(ri), (∂/∂xij)β̂(y,X) = Âejψ(ri)− ÂX>eiψ′(ri)β̂j ,
for ri = yi − x>i β̂, ∀i ∈ [n],j ∈ [p] where ej ∈ Rp and ei ∈ Rn are canonical basis vectors.53
3. The paper obtains a stochastic representation for the residual yi−x>i β̂ for some fixed i = 1, ..., n,54 extending some results of [10] on unregularized M -estimators to penalized ones as in (1). In55 short, for each i = 1, ..., n the i-th residual satisfies ri = yi − x>i β̂56
ri + (d̂f/ trV )ψ(ri) ≈ εi + Zi‖Σ1/2(β̂ − β∗)‖ (4) where Zi ∼ N(0, 1) is independent of εi. This stochastic representation is the motivation for57 the criterion (2) as the amplitude of the normal part in the right-hand side is proportional to the58 out-of-sample error ‖Σ1/2(β̂ − β∗)‖ that we wish to minimize, while the variance of the noise59 εi does not depend on the choice of (ρ, g).60
Simulated data in Figure 2 confirms that the stochastic representation for the i-th residual ri =61 yi − x>i β̂ is accurate. Our working assumption throughout the paper is the following.62 Assumption 1.1. For constants γ, µ > 0 independent of n, p we have p/n ≤ γ, the loss ρ : R→ R63 is convex with a unique minimizer at 0, continuously differentiable and its derivative ψ = ρ′ is64 1-Lipschitz. The design matrix X has iid N(0,Σ) rows for some invertible covariance Σ and the65 noise ε is independent of X with continuous distribution. The penalty g : Rp → R is µ-strongly66 convex w.r.t. Σ in the sense that b 7→ g(b)− (µ/2)b>Σb is convex in b ∈ Rp.67
Throughout the paper, we consider a sequence (say, indexed by n) of regression problems with p,68 β∗, Σ and the loss-penalty pair (ρ, g) depending implicitly on n. For some deterministic sequence69 (an), the stochastically bounded notation OP (an) in this context may hide constants depending on70 γ, µ only, that is, OP (an) denotes a sequence of random variables Wn such that for any ε > 0 there71 exists K depending on (ε, γ, µ) satisfying P(|Wn| ≥ Kan) ≤ ε.72 Since Assumption 1.1 requires p/n ≤ γ, the Bolzano-Weierstrass theorem lets us extract a subse-73 quence of regression problems such that p/n→ γ′ along this subsequence, for some constant γ. This74 is the asymptotic regime we have in mind throughout the paper, although our results do not require a75 specific limit for the ratio p/n. For some results, we will require the following additional assumption76 which is satisfied by robust loss functions and penalty that shrink towards 0.77
Assumption 1.2. The penalty is minimized at 0, that is, g(0) = minb∈Rp g(b); the loss is Lipschitz as78 in |ψ| ≤M for some constant M independent of n, p; the signal is bounded as in ‖Σ1/2β∗‖2 ≤M .79
1.2 Related works80
The context of the present work is the study of M -estimators in the regime pn has a finite limit. This81 literature pioneered in [2, 10, 9, 15] typically describes the subtle behavior of β̂ in this regime by82 solving a system of nonlinear equations. This system typically depends on a prior distribution for the83 components of β∗, and either depends on the covariance Σ [7] or assume Σ = Ip [2, 16, 6, among84 many others]. Solutions to the nonlinear system are a powerful tool to understand β̂ in theory, e.g.,85 to characterize the deterministic limit of ‖Σ1/2(β̂ − β∗)‖, see e.g., the general results in [6] for the86 square loss and [16] for general loss-penalty pairs. However, since the system and its solution depend87 on unobservable quantity (Σ and prior on β∗), the system solution is not directly usable for practical88 purposes such as parameter tuning.89
The present work distinguishes itself from most of this literature as the goal is to describe the behavior90 of β̂ using observable quantities that only depend on the data (y,X) (and not unobservable ones such91 as Σ or a prior distribution on β∗ that appear in the aforementioned nonlinear system of equations).92 As we will see this view lets us perform adaptive tuning of parameters in a fully adaptive manner93 using the criterion (2). The criterion (2) appeared in previous works for the square loss only: [1, 12]94 studied (2) for the Lasso with Σ = Ip and [3, Section 3] for the square loss and general penalty95 (note that for the square loss ρ(u) = u2/2, (2) reduces to n2‖r‖2/(n− d̂f)2 due to ψ(u) = u and96 tr[V ] = n − d̂f. The property ψ(u) = u of the square loss hides the subtle interplay between97 r, ψ(r), d̂f and tr[V ] in (2) for ρ different than the square loss). A criterion different from (2) is98 studied in [12, 3] to estimate the out-of-sample error. That criterion has the drawback to require the99 knowledge of Σ, unlike (2) which is fully adaptive.100
This work leverages probabilistic results on functions of standard normal random variables [4][3,101 §6, §7] which are consequences of Stein’s formula [14]. Consequently, the main limitation of our102 work is that it currently requires Gaussian design for the probabilistic results (on the other hand, the103 differentiability result (5) is deterministic and does not rely on any probabilistic assumption).104
2 Differentiability of regularized M-estimators105
The first step towards the study of the criterion (2) is to justify the almost sure existence of the106 derivatives of β̂ that appear in (2) through the scalar scalar d̂f and the matrix V in (2). Although the107 criterion (2) only involves the derivatives of β̂(y,X) with respect to y for a fixedX , the proof of108
our results rely on the interplay between the derivatives with respect to y and with respect toX: this109 differentiability structure of M -estimators is the content of the following result.110
Theorem 2.1. Let Assumption 1.1 be fulfilled. For almost every (y,X) the map (y,X) 7→ β̂(y,X)111 is differentiable at (y,X) and there exists a matrix Â ∈ Rp×p with ‖Σ1/2ÂΣ1/2‖op ≤ (nµ)−1 s.t.112
(∂/∂yi)β̂(y,X) = ÂX >eiψ ′(ri),
(∂/∂xij)β̂(y,X) = Âejψ(ri)− ÂX>eiψ′(ri)β̂j , where ri = yi − x>i β̂, (5)
ei ∈ Rn, ej ∈ Rp are canonical basis vectors , ψ := ρ′ and ψ′ denote the derivatives. Furthermore,113
d̂f = tr[X(∂/∂y)β̂] = tr[XÂX diag{ψ′(r)}], (6)
V = diag{ψ′(r)}(In −X(∂/∂y)β̂) = diag{ψ′(r)} − diag{ψ′(r)}XÂX diag{ψ′(r)}. (7)
satisfy 0 ≤ d̂f ≤ n and 0 ≤ tr[V ] ≤ n.114
Since the same matrix Â appears in both the derivatives with respect to yi and to xij , (5) provides115 relationship between (∂/∂yi)β̂ and (∂/∂xij)β̂, for instance (∂/∂xij)β̂ = Âejψ(ri)− β̂j(∂/∂yi)β̂.116 Although the matrix Â is not explicit for arbitrary loss-penalty pair, closed-form expressions are117 available for particular examples such as the Elastic-Net penalty as discussed in Section 6.118 Remark 2.1. For the square loss ρ(u) = u2/2, the differentiability formulae (5) reduce to119
(∂/∂yl)β̂(y,X) = ÂX >el, (∂/∂xij)β̂(y,X) = Âej(yi − x>i β̂)− ÂX >eiβ̂j
for most every (y,X) and some matrix Â ∈ Rp×p depending on (y,X), since in this case ψ′ = 1.120
In the simple case where g is twice continuously differentiable, (5) follows [4] with121
Â = ( X> diag{ψ′(r)}X + n∇2g(β̂) )−1 (8)
by differentiating the KKT conditions X>ψ(y −Xβ̂) = n∇g(β̂). To illustrate why this is true, provided that β̂(y,X) is differentiable, if (y(t),X(t)) are smooth perturbations of (y,X) with (y(0),X(0)) = (y,X) and ddt (y(t),X(t))|t=0 = (ẏ, Ẋ), differentiation of X(t)
>ψ(y(t) − X(t)β̂(y(t),X(t))) = n∇g(β̂(y(t),X(t))) at t = 0 and the chain rule yields
Ẋ>ψ(r)−X> diag{ψ′(r)}(ẏ − Ẋβ̂(y,X)) = Â−1 ddt β̂(y(t),X(t)) ∣∣ t=0
with Â in (8). This gives (5) if the penalty g is twice-differentiable. Theorem 2.1 reveals that for122 arbitrary convex penalty functions including non-differentiable ones, the differentiability structure123 (5) always holds, as in the case of twice differentiable penalty g, even for penalty functions such as124 g(b) = µ‖b‖2/2 + λ‖mat(b)‖nuc where mat : Rp → Rd1×d2 is a linear isomorphism to the space of125 d1 × d2 matrices and ‖ · ‖nuc is the nuclear norm: in this case by Theorem 2.1 there exists a matrix126 Â ∈ Rp×p such that (5) holds although no closed-form expression for Â is known.127 The representation (5) is a powerful tool as it provides explicit derivatives of quantities of interest128 such as r = y −Xβ̂, ‖ψ(r)‖2 or ‖Σ1/2(β̂ − β∗)‖2. These explicit derivatives can then be used in129 probabilistic identities and inequalities that involve derivatives, for instance Stein’s formulae [14],130 the Gaussian Poincaré inequalty [5, Theorem 3.20], or normal approximations [8, 4].131 Remark 2.2. Similar derivative formulae hold if an intercept is included in the minimization, as in132
( β̂0(y,X), β̂(y,X) ) = argmin b0∈R,b∈Rp 1 n n∑
i=1
ρ(yi − b0 − x>i b) + g(b) (9)
Let Assumption 1.1 be fulfilled, and assume further ‖ψ′(r)‖2 > 0 with r := y− 1nβ̂0−x>i β̂ where133 1n = (1, ..., 1)
> ∈ Rn. For almost every (y,X) the map (y,X) 7→ β̂(y,X) is differentiable at134 (y,X), and there exists Â ∈ Rp×p depending on (y,X) with ‖Σ1/2ÂΣ1/2‖op ≤ (nµ)−1 such that135
136
(∂/∂yi)β̂(y,X) = ÂX >Ψ′ei, (∂/∂xij)β̂(y,X) = Âejψ(ri)− ÂX>Ψ′eiβ̂j , (10)
where ei ∈ Rn, ej ∈ Rp are canonical basis vectors, ψ = ρ′ and Ψ′ := diag{ψ′(r)} −137 ψ′(r)ψ′(r)>/ ∑ i∈[n] ψ ′(ri).138
3 Distribution of individual residuals139
We now turn to the distribution of a single residual ri = yi − x>i β̂ for some fixed observation140 i ∈ {1, ..., n} (for instance, fix i = 1). By leveraging the differentiability structure (5) and the normal141 approximation from [4], the following result provides a clear picture of the distribution of ri.142 Theorem 3.1. Let Assumption 1.1 be fulfilled and let Â ∈ Rp×p be given by Theorem 2.1. Then for143 every i = 1, ..., n there exists Zi ∼ N(0, 1) such that144 ∣∣∣ ( ri+tr[ΣÂ]ψ(ri) ) − ( εi+‖Σ1/2(β̂−β∗)‖Zi
)∣∣∣ ≤ OP (n−1/4)(|ψ(εi)|+‖Σ1/2(β̂−β∗)‖) (11) Furthermore, if εi has a fixed distribution F , there exists a bivariate variable (ε̃ni , Z̃ n i ) converging in145 distribution to the product measure F ⊗N(0, 1) such that146 ri + tr[ΣÂ]ψ(ri) = ε̃ n i + ‖Σ
1/2(β̂ − β∗)‖Z̃ni . (12) If εi has a fixed distribution F and Assumption 1.2 holds then |ψ(εi)|+ ‖Σ1/2(β̂ − β∗)‖ = OP (1).147
Theorem 3.1 is a formal statement regarding the informal normal approximation148
ζi := ri + tr[ΣÂ]ψ(ri)− εi ‖Σ1/2(β̂ − β∗)‖ ≈ N(0, 1). (13)
Simulations in Figure 2 confirm the normality of ζi for the Huber loss with Elastic-Net penalty149 and four combinations of tuning parameters. For the square loss ρ(u) = u2/2, because ψ(u) = u,150 asymptotic normality of the residuals hold in the following form.151 Theorem 3.2. Let Assumption 1.1 hold with ρ(u) = u2/2 and ε ∼ N(0, σ2In). Then for i = 1,152
(σ2 + ‖Σ1/2(β̂ − β∗)‖2)−1/2(1 + tr[ΣÂ])(yi − x>i β̂)→d N(0, 1) as n→ +∞. (14)
It is informative to provide a sketch of the proof of Theorem 3.1 explain the appearance of ψ(ri) and153 tr[ΣÂ] in the normal approximation results (11) and (13). A variant of the normal approximation154 of [4] proved in the supplement states that for a differentiable function f : Rq → Rq \ {0} and155 z ∼ N(0, Iq), there exists Z ∼ N(0, 1) such hat156
E [∣∣∣ f(z) >z − ∑q k=1(∂/∂zk)fk(z)
‖f(z)‖ − Z
∣∣∣ 2] ≤ C1E [∑q k=1 ‖(∂/∂zk)f(z)‖2 ‖f(z)‖2 ] . (15)
Some technical hurdles aside, the proof sketch is the following: Apply the previous display to q = p,157 z = Σ−1/2xi conditionally on (ε, (xl)l∈[n]\{i}) and to f(z) = Σ
1/2(β̂ − β∗) in the simple case158 where β∗ = 0 (this amounts to performing a change of variable by translation of β̂ to β̂ − β∗). Then159 the right-hand side of the previous display is negligible in probability compared to Z, and in the160 left-hand side f(z)>z = x>i (β̂ − β ∗) and ∑q k=1(∂/∂zk)fk(z) ≈ tr[ΣÂ]ψ(ri) as the second term161 in (5) is negligible. This completes the sketch of the proof of (13).162
Proximal operator representation. From the above asymptotic normality results, a stochastic representation for the i-th residual ri = yi − x>i β̂ can be obtained as follows: With prox[tρ](u) the proximal operator of x 7→ tρ(x) defined as the unique solution z ∈ R of equation z + tψ(z) = u,
ri = yi − x>i β̂ = prox[t̂ρ] ( ε̃ni + ‖Σ 1/2(β̂ − β∗)‖Z̃ni )
with t̂ = tr[ΣÂ]
where (ε̃ni , Z̃ n i ) converges in distribution to product measure F ⊗N(0, 1) where F is the law of εi.163
4 A proxy of the out-of-sample error if Σ is known164
The approximations of the previous sections for ri + tr[ΣÂ]ψ(ri) and the fact that εi is independent165 of Zi ∼ N(0, 1) in (11) suggest that (ri + tr[ΣÂ]ψ(ri))2 ≈ ε2i + ‖Σ
1/2(β̂ − β∗)‖2Z2i ; and166 averaging over {1, ..., n} one can hope for the approximation ‖r + tr[ΣÂ]ψ(r)‖2/n ≈ ‖ε‖2/n+167 ‖Σ1/2(β̂ − β∗)‖2. The following result makes this heuristic precise.168
Theorem 4.1. Let Assumption 1.1 be fulfilled and Â be given by Theorem 2.1. Then169
‖Σ1/2(β̂ − β∗)‖2 + ‖ε‖2/n = ∥∥r + tr[ΣÂ]ψ(r) ∥∥2/n+OP (n−1/2) Rem,
where Rem := ‖Σ1/2(β̂−β∗)‖2 + 1n‖ψ(r)‖ 2 + (‖Σ1/2(β̂−β∗)‖2 + 1n‖ψ(r)‖ 2)1/2‖ 1√ n ε‖. Thus
‖Σ1/2(β̂ − β∗)‖2 + ‖ε‖2/n = (1 +OP (n−1/2)) ∥∥r + tr[ΣÂ]ψ(r) ∥∥2/n.
Theorem 4.1 provides a first candidate, ∥∥r + tr[ΣÂ]ψ(r) ∥∥2/n to estimate170
‖Σ1/2(β̂ − β∗)‖2 + ‖ε‖2/n. (16)
Estimation of (16) is useful as ‖ε‖2/n is independent of the choice of the estimator β̂ and in particular171 independent of the chosen loss-penalty pair in (1). Given two or more estimators (1), choosing the172 one with smallest ∥∥r + tr[ΣÂ]ψ(r) ∥∥2 is thus a good proxy for minimizing the out-of-sample error.173
Corollary 4.2. Let β̂, β̃ be two M -estimators (1) Assumption 1.1 with loss-penalty pair (ρ, g) and174 (ρ̃, g̃) respectively. Assume that both satisfy Assumption 1.1 and let ψ = ρ′ and ψ̃ = ρ̃′. Let175 r = y −Xβ̂, r̃ = y −Xβ̃ be the residuals, Â, Ã be the corresponding matrices of size p × p176 given by Theorem 2.1. Further assume that both estimators satisfy Assumption 1.2 and that ε has iid177 coordinates independent with E[|εi|1+q] ≤M for constants q ∈ (0, 1),M > 0 independent of n, p.178 Let Ω = {‖XΣ−1/2‖op ≤ 2 √ n+ √ p} ∩ {‖ε‖2 ≤ n2/(1+q)}. Then for any η > 0 independent of179 n, p there exists C(γ, µ, η, q,M) > 0 depending only on {γ, µ, η, q,M} such that180
P ( ‖Σ1/2(β̂ − β∗)‖2 − ‖Σ1/2(β̃ − β∗)‖2 > η, ‖r + tr[ΣÂ]ψ(r)‖2 ≤ ‖r̃ + tr[ΣÃ]ψ̃(r̃)‖2 )
≤ C(γ, µ, η, q,M)n−q/(1+q) + P(Ωc)→ 0.
Provided that the noise random variables εi have at least 1 + q moments, Corollary 4.2 implies181 that with probability approaching one given two M -estimators β̂ and β̃, choosing the estimator182 corresponding to the smallest criteria among ‖r + tr[ΣÂ]r‖2 and ‖r̃ + tr[ΣÃ]r̃‖2 leads to the183 smallest out-of-sample error, up to any small constant η > 0. This allows noise random variables εi184 with infinite variance. A similar result can be obtained to select among K different M -estimators (1).185
Corollary 4.3. As in Corollary 4.2, assume E[|εi|1+q] ≤M and let β̂1, ..., β̂K be M -estimators of186 the form (1) with loss-penalty pair (ρk, gk) satisfying Assumptions 1.1 and 1.2. For each k = 1, ...,K,187 let rk = y−Xβ̂k be the residuals and Âk be the corresponding matrix of size p×p from Theorem 2.1.188 Let k̂ ∈ argmink=1,...,K ‖rk + tr[ΣÂk]ψk(rk)‖ where ψk = ρ′k. Then if (γ, µ, η, q,M) are189 constants independent of n, p190
P ( ‖Σ1/2(β̂k̂ − β ∗)‖2 > mink=1,...,K ‖Σ1/2(β̂k − β ∗)‖2 + η ) → 0 if K = o(nq/(1+q)).
Given K different loss-penalty pairs and the corresponding M -estimators in (1), minimizing the191 criterion ‖r + tr[ΣÂ]r‖ thus provably selects a loss-penalty pair that leads to an optimal out-192 of-sample error, up to an arbitrary small constant η > 0 independent of n, p. The requirement193 K = o(nq/(1+q)) means that the cardinality of the collection of M -estimators to select from should194 grow more slowly than a power of n. This is typically satisfied for default tuning parameter grids in195 popular libraries (e.g., sklearn.linear_model.Lasso [13]) with tuning parameters evenly spaced196 in a log-scale that consequently have cardinality logarithmic in the parameter range. The major197 drawback of the criterion ‖r + tr[ΣÂ]r‖ is the dependence through tr[ΣÂ] on the covariance Σ198 of the design, which is typically unknown. The next section introduces an estimator of tr[ΣÂ] that199 does not require the knowledge of Σ.200

heading: 5 Degrees of freedom and estimating tr[ΣÂ] without the knowledge of Σ201
text: This section focuses on estimating tr[ΣÂ]. The matrix Â from Theorem 2.1 can estimated from the data (y,X) in the sense that Â is a measurable function of (y,X) (thanks to the observation that derivatives are limits, and limits of measurable functions are again measurable). The difficulty is thus to estimate tr[ΣÂ] without the knowledge of Σ. To illustrate this difficulty, consider Ridge regression with square loss ρ(u) = u2/2 and penalty g(b) = τ‖b‖2/2. Then β̂(y,X) = (X>X + τnIp) −1X>y and Â in Theorem 2.1 is given explicitly by Â = (X>X + τnIp)−1 and
tr[ΣÂ] = tr[(G>G+ nτΣ−1)−1], whereG = XΣ−1/2.
Above,G is a random matrix with iid N(0, 1) entries the value of tr[ΣÂ] is highly dependent on the202 spectrum of Σ−1. In this particular case, the limit of tr[(G>G+ nτΣ−1)−1] can be obtained using203 random matrix theory [11] as the limiting behavior of the Stieltjes transform of G>G/n + τΣ−1204 and its spectral distribution is known; however the limit of the spetral distribution depends on the205 spectrum of τΣ−1. This is not desirable here as we wish to construct estimators that require no206 knowledge on Σ. For more involved loss-penalty pairs such as the Elastic-Net in Example 6.1, such207 random matrix theory results do not apply as tr[ΣÂ] depends on the random support of β̂.208
Instead, we do not rely on known random matrix theory results. With the matrix Â ∈ Rp×p given by209 Theorem 2.1, our proposal to estimate tr[ΣÂ] is the ratio d̂f/ tr[V ] with d̂f and V in (6)-(7). Both210 the scalar d̂f and the matrix V ∈ Rn×n are observable; in particular they do not depend on Σ.211 Theorem 5.1. Let Assumption 1.1 be fulfilled and Â be given by Theorem 2.1. Then212
E[| tr[ΣÂ] tr[V ]/n− d̂f/n|] ≤ C2(γ, µ)n−1/2. (17)
Simulations in Figure 3 and Table 1 confirm that the approximation tr[ΣÂ] ≈ d̂f/ tr[V ] is accurate213 for the Huber loss with Elastic-Net penalty. For the square loss, ψ′ = 1 and tr[V ] = n− d̂f so that214 (17) becomes E|(1− d̂f/n)(1 + tr[ΣÂ])− 1| ≤ C3(γ, µ)n−1/2 and the following result holds.215 Corollary 5.2. Let Assumption 1.1 be fulfilled with ρ(u) = u2/2 and ε ∼ N(0, σ2In). Then216 (1−d̂f/n)(1+tr[ΣÂ])→P 1 and the normality (14) holds with 1+tr[ΣÂ] replaced by (1−d̂f/n)−1.217
For general loss ρ, the criterion (2) replaces tr[ΣÂ] by d̂f/ tr[V ] in the proxy of the out-of-sample218 error ‖r+ tr[ΣÂ]ψ(r)‖2 studied in the previous section. Thanks to (17), this replacement preserves219 the good properties of ‖r + tr[ΣÂ]ψ(r)‖2 proved in Corollaries 4.2 and 4.3.220 Theorem 5.3. For k = 1, ...,K, let (ρk, gk) be a loss-penalty pair satisfying Assumptions 1.1 and 1.2 with ψk = ρ′k, let β̂k, rk, Âk be the corresponding M -estimator residual vector and matrix of size p × p given by Theorem 2.1 as in Corollary 4.3 and let d̂fk = tr[XAkX> diag{ψ′k(rk)}] and V k = diag{ψ′k(rk)}(In −XAkX
> diag{ψ′k(rk)}). For a small constant η > 0 independent of n, p, say η = 0.05, define
k̂ ∈ argmin k=1,...,K
∥∥∥rk + d̂fk
tr[V k] ψk(rk)
∥∥∥ 2
subject to 1
n
n∑
i=1
ψ′k(rki) ≥ η.
0 .0
0 3
2
0 .0
0 4
7
0 .0
0 7
1
0 .0
1 1
0 .0
1 6
0 .0
2 4
0 .0
3 6
0 .0
5 4
0 .0
8 1
0 .1
2
0 .1
8
0 .2
7
0 .4
1
λ
1e-10
1e-07
1e-05
0.001
0.01
0.032
0.063
0.1
τ
n̂/n
0.0
0.2
0.4
0.6
0.8
1.0
0 .0
0 3
2
0 .0
0 4
7
0 .0
0 7
1
0 .0
1 1
0 .0
1 6
0 .0
2 4
0 .0
3 6
0 .0
5 4
0 .0
8 1
0 .1
2
0 .1
8
0 .2
7
0 .4
1
λ
1e-10
1e-07
1e-05
0.001
0.01
0.032
0.063
0.1
τ
d̂f/n
0.0
0.2
0.4
0.6
0.8
1.0
0 .0
0 3
2
0 .0
0 4
7
0 .0
0 7
1
0 .0
1 1
0 .0
1 6
0 .0
2 4
0 .0
3 6
0 .0
5 4
0 .0
8 1
0 .1
2
0 .1
8
0 .2
7
0 .4
1
λ
1e-10
1e-07
1e-05
0.001
0.01
0.032
0.063 0.1 τ
tr[V ]/n
0.0
0.2
0.4
0.6
0.8
1.0
If εi has 1 + q moments in the sense that E[|εi|1+q] ≤ M for constants q ∈ (0, 1),M > 0. If221 (M, q, η, µ, γ) and η̃ > 0 are independent of n, p then222
P ( ‖Σ1/2(β̂k̂−β
∗)‖ > min k=1,...,K: 1n ∑n i=1 ψ ′ k(rki)≥η ‖Σ1/2(β̂k−β ∗)‖+η̃
) → 0 if K = o(nq/(1+q)).
Figure 1 illustrates on simulations the success of the criterion (2) over a grid of tuning parameters223 for M -estimators with the Huber loss and Elastic-Net penalty. The criterion (2) is thus successful224 at selecting a M -estimator with smallest out-of-sample error up to an additive constant η̃, among225 those M -estimators indexed in {1, ...,K} that are such that 1n ∑n i=1 ψ ′ k(rki) ≥ η. On the one hand226
it is unclear to us whether the restriction 1n ∑n i=1 ψ ′ k(rki) ≥ η; on the other hand there is a practical227
meaning in excluding M -estimators with small 1n ∑n i=1 ψ ′ k(rki): For the Huber loss H(u) := u 2/2228
for |u| ≤ 1 and |u| − 1/2 for |u| ≥ 1 the quantity 1n ∑n i=1 ψ ′ k(rki) is the number of of data points229 in {1, ..., n} such that the residual yi − x>i β̂k fall within the quadratic regime of the loss function.230 Observations i ∈ {1, ..., n} that fall in the linear regime of the loss are excluded from the fit, in the231 sense that for some iwith rki = yi−x>i β̂k > 1, replacing yi by ỹi = yi+1000 (or any positive value)232 does not change the M -estimator solution β̂k (this can be seen from the KKT conditions directly,233 or by integration the derivative with respect to yi in (5)). Thus the constraint 1n ∑n i=1 ψ ′ k(rki) ≥ η234 requires that at most a constant fraction of the observations are excluded from the fit (or equivalently,235 at least a constant fraction of the n observations participate in the fit). For scaled versions of the236 Huber loss, ρk(u) = a2H(a−1u) for some a > 0, the value n̂ = 1n ∑n i=1 ψ ′ k(rki) again counts237 the number of residuals falling in the quadratic regime of the loss, i.e., the number of observations238 participating in the fit. The heatmaps of Figure 3 illustrate n̂ in a simulation for a wide range of239 parameters. Similarly, for smooth robust loss functions such as ρk(u) = √ 1 + u2, the constraint240
1 n ∑n i=1 ψ ′ k(rki) ≥ η requires that at most a constant fraction of the n observations are such that241
ψ′k(rki) < η/2, i.e., such that the second derivative ψ ′ k is too small (and the loss ρk too flat).242
Theorems 2.1, 3.2, 4.1 and 5.1 provide our general results applicable to a single regularized M -243 estimator (1) while corollaries such as Theorem 5.3 are obtained using the union bound. The next244
section specializes our results and notation to the Huber loss with Elastic-Net penalty and details the245 simulation setup used in the figures.246
6 Example and simulation setting: Huber loss with Elastic-Net penalty247
In simulations and in the example below, we focus on the loss-penalty pair248 ρ(u; Λ) = Λ2H(Λ−1u), g(b;λ, τ) = λ‖b‖1 + (τ/2)‖b‖22 (18) for tuning parameters Λ, λ, τ ≥ 0 where H(u) := u2/2 for |u| ≤ 1 and |u| − 1/2 for |u| ≥ 1.249 Example 6.1. With (ρ, g) in (18), matrix Â in (5) matrix V in (7) and d̂f in (6) we have250
ÂŜ,Ŝ = (X > Ŝ diag{ψ′(r)}X Ŝ + nτI p̂) −1, Ai,j = 0 if i 6∈ Ŝ or j 6∈ Ŝ, V = diag{ψ′(r)} − diag{ψ′(r)}X Ŝ(X > Ŝ diag{ψ′(r)}X Ŝ + nτI p̂) −1X> Ŝ diag{ψ′(r)}, d̂f = tr[X Ŝ(X > Ŝ diag{ψ′(r)}X Ŝ + nτI p̂) −1X> Ŝ diag{ψ′(r)}], (19)
where Ŝ is the active set {j ∈ [p] : β̂j 6= 0} and p̂ is the size of Ŝ; X Ŝ is the submatrix of X251 selecting columns with index in Ŝ and ÂŜ,Ŝ is the submatrix of Â with entries indexed in Ŝ × Ŝ.252
The identities (19) are proved in [3, §2.6]. Simulations in Figures 1 to 3 and Table 1 illustrate typical253 values for d̂f, tr[V ], tr[ΣÂ], the out-of-sample error and the criterion (2), n̂ = ∑n i=1 ψ
′(ri) and254 p̂ = |Ŝ| under anisotropic Gaussian design and heavy-tailed εi. The simulation setup is as follows.255 Data Generation Process. Simulation data are generated from a linear model y = Xβ∗ + ε with256 anisotropic Gaussian design Σ and heavy-tail noise vector ε. The design matrix X has n = 1001257 rows and p = 1000 columns. Each row of X is i.i.d. N(0,Σ), with the same Σ across all258 repetitions, generated once by Σ = R>R/(2p) withR ∈ R2p×p being a Rademacher matrix with259 i.i.d. entries P(Rij = ±1) = 12 . The true signal vector β
∗ ∈ Rp has its first 100 coordinates set to260 p1/2/100 = √ 10/10 and the rest 900 coordinates set to 0. The noise vector ε ∈ Rn has i.i.d. entries261 from the t-distribution with 2 degrees of freedom (so that Var[εi] =∞, i.e., εi is heavy-tailed).262 Estimation Process. Each dataset (y,X) is fitted by a Huber Elastic-Net estimator with263 loss-penalty pair in (18). We focus on 2d heatmaps with respect to the two penalty parame-264 ters (λ, τ) of the penalty; to this end the Huber loss parameter Λ is set to Λ = 0.054n1/2265 and a grid for (λ, τ) in then set so that d̂f/n varies on the grid from 0 to 1 (cf. the mid-266 dle heatmap in Figure 3). The Elastic-Net penalty g(b;λ, τ) = λ‖b‖1 + (τ/2)‖b‖22 is used267 with (λ, τ) ∈ {(0.036, 10−10), (0.054, 0.01), (0.036, 0.01), (0.024, 0.1)} in Figure 2 and Table 1,268 (λ, τ) ∈ [0.0032, 0.41]×{10−10} in Figure 3, and (λ, τ) ∈ [0.0032, 0.041]×[10−10, 0.1] in Figure 1.269 More simulation results are provided in the supplementary materials.270
References271 [1] Mohsen Bayati, Murat A Erdogdu, and Andrea Montanari. Estimating lasso risk and noise level.272 In Advances in Neural Information Processing Systems, pages 944–952, 2013.273
[2] Mohsen Bayati and Andrea Montanari. The lasso risk for gaussian matrices. IEEE Transactions274 on Information Theory, 58(4):1997–2017, 2012.275
[3] Pierre C Bellec. Out-of-sample error estimate for robust m-estimators with convex penalty.276 arXiv:2008.11840, 2020.277
[4] Pierre C Bellec and Cun-Hui Zhang. Second order poincare inequalities and de-biasing arbitrary278 convex regularizers when p/n→ γ. arXiv:1912.11943, 2019.279
[5] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. Concentration inequalities: A280 nonasymptotic theory of independence. Oxford University Press, 2013.281
[6] Michael Celentano and Andrea Montanari. Fundamental barriers to high-dimensional regression282 with convex penalties. arXiv preprint arXiv:1903.10603, 2019.283
[7] Michael Celentano, Andrea Montanari, and Yuting Wei. The lasso with general gaussian designs284 with applications to hypothesis testing. arXiv preprint arXiv:2007.13716, 2020.285
[8] Sourav Chatterjee. Fluctuations of eigenvalues and second order poincaré inequalities. Proba-286 bility Theory and Related Fields, 143(1):1–40, 2009.287
[9] David Donoho and Andrea Montanari. High dimensional robust m-estimation: Asymptotic288 variance via approximate message passing. Probability Theory and Related Fields, 166(3-289 4):935–969, 2016.290
[10] Noureddine El Karoui, Derek Bean, Peter J Bickel, Chinghway Lim, and Bin Yu. On robust291 regression with high-dimensional predictors. Proceedings of the National Academy of Sciences,292 110(36):14557–14562, 2013.293
[11] Vladimir A Marčenko and Leonid Andreevich Pastur. Distribution of eigenvalues for some sets294 of random matrices. Mathematics of the USSR-Sbornik, 1(4):457, 1967.295
[12] Léo Miolane and Andrea Montanari. The distribution of the lasso: Uniform control over sparse296 balls and adaptive parameter tuning. arXiv preprint arXiv:1811.01212, 2018.297
[13] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion,298 Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Scikit-299 learn: Machine learning in python. the Journal of machine Learning research, 12:2825–2830,300 2011.301
[14] Charles M Stein. Estimation of the mean of a multivariate normal distribution. The annals of302 Statistics, pages 1135–1151, 1981.303
[15] Mihailo Stojnic. A framework to characterize performance of lasso algorithms. arXiv preprint304 arXiv:1303.7291, 2013.305
[16] Christos Thrampoulidis, Ehsan Abbasi, and Babak Hassibi. Precise error analysis of regularized306 m-estimators in high dimensions. IEEE Transactions on Information Theory, 64(8):5592–5628,307 2018.308
Checklist309
1. For all authors...310 (a) Do the main claims made in the abstract and introduction accurately reflect the paper’s311 contributions and scope? [Yes]312 (b) Did you describe the limitations of your work? [Yes] See Assumptions 1.1 and 1.2 and313 the limitations mentioned in Section 1.2314 (c) Did you discuss any potential negative societal impacts of your work? [N/A]315 (d) Have you read the ethics review guidelines and ensured that your paper conforms to316 them? [Yes]317 2. If you are including theoretical results...318
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Assump-319 tions 1.1 and 1.2.320
(b) Did you include complete proofs of all theoretical results? [Yes] See supplementary321 material.322
3. If you ran experiments...323 (a) Did you include the code, data, and instructions needed to reproduce the main exper-324 imental results (either in the supplemental material or as a URL)? [Yes] See code in325 supplementary material.326
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they327 were chosen)? [Yes] The code is also provided.328
(c) Did you report error bars (e.g., with respect to the random seed after running experi-329 ments multiple times)? [Yes]330
(d) Did you include the total amount of compute and the type of resources used (e.g., type331 of GPUs, internal cluster, or cloud provider)? [Yes] It takes an Amazon EC2 server332 approximately 40 hours to generate all our simulation results. This is also mentioned333 in supplementary.334
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...335 (a) If your work uses existing assets, did you cite the creators? [N/A] Simulations are336 implemented using Python.337 (b) Did you mention the license of the assets? [N/A]338 (c) Did you include any new assets either in the supplemental material or as a URL? [N/A]339
340
(d) Did you discuss whether and how consent was obtained from people whose data you’re341 using/curating? [N/A]342
(e) Did you discuss whether the data you are using/curating contains personally identifiable343 information or offensive content? [N/A] Simulated data only.344
5. If you used crowdsourcing or conducted research with human subjects...345 (a) Did you include the full text of instructions given to participants and screenshots, if346 applicable? [N/A]347 (b) Did you describe any potential participant risks, with links to Institutional Review348 Board (IRB) approvals, if applicable? [N/A]349 (c) Did you include the estimated hourly wage paid to participants and the total amount350 spent on participant compensation? [N/A]351

Review summaries:
This paper studies the derivatives of regularized M-estimators and distribution of the residual.  These results give rises to a new novel adaptive criterion for tuning loss and penalty function. Numerical experiments are conducted on various models to validate the theoretical result.

The paper is dedicated to the analysis of robust linear regression, where the proposed estimation is based on sparse M-estimation. The statistical criterion corresponds to a penalised M-estimation problem, where both the non-penalised loss function and penalty function are assumed convex. The contributions of the paper can be summarised as follows: - the authors provide a criterion for the optimal selection of the non-penalised loss and penalty function so that some out-of-sample distance, say D, is minimal.  - the authors quantify the behaviour of the penalised estimator \hat\beta using quantities depending on the observations only (idest no unknown quantities such as true sparse support) in a regime where p/n has finite limit. - some explicit derivative formulas of the penalised estimator are provided.  - the distribution of the estimated residuals is provided, which allows for a selection of the tuning parameter for regularization. - the authors derive an approximation formula of the distance D.  

The authors study M-estimators with convex penalties. They give algebraic forms for derivatives of the estimator with respect to response $y$ and the design matrix $X$. They describe the distribution of residuals or certain functions of residuals. Given a loss function and penalty (with tuning parameters) they propose a criterion that can effectively approximate the out-of-sample error. This criterion can be computed without the knowledge of the covariance matrix $\Sigma$ (rows of $X \sim N(0,\Sigma)$). They give finite sample bounds that show how fast the criterion converges to the out-of-sample error.   

Finding the right parameter when doing robust estimation is a recurring problem in the robust community. This article propose to choose the parameters in a robust regularized linear regression model by minimizing a completely data-driven criterion. This criterion is shown to be a good approximation to the out of sample error, this is also illustrated on practical experiments.

This paper studies M-estimators in linear models, with convex and gradient-Lipschitz loss function and a convex penalty.   Main contributions:  1. The authors provide a criterion to select tuning parameters (the loss function and the penalty function) of M-estimators, such that minimizing the criterion provides a proxy for minimizing the out-of-sample error, since the criterion is an approximation of the out-of-sample error up to a constant shift. One advantage of this criterion is that it does not require the prior knowledge of the covariance matrix of the design and the noise distribution.  2. The authors provide simulation results that confirm theoretical findings.