In a recent article (Proc. Natl. Acad. Sci., 110(36), 14557-14562), El Karoui et al. study the distribution of robust regression estimators in the regime in which the number of parameters p is of the same order as the number of samples n. Using numerical simulations and 'highly plausible' heuristic arguments, they unveil a striking new phenomenon. Namely, the regression coefficients contain an extra Gaussian noise component that is not explained by classical concepts such as the Fisher information matrix. We show here that that this phenomenon can be characterized rigorously using techniques developed by the authors to analyze the Lasso estimator under high-dimensional asymptotics. We introduce an approximate message passing (AMP) algorithm to compute M-estimators and deploy state evolution to evaluate the operating characteristics of AMP and so also M-estimates. Our analysis clarifies that the 'extra Gaussian noise' encountered in this problem is fundamentally similar to phenomena already studied for regularized least squares in the setting n < p. Consider the traditional linear regression model with Y = (Y 1 , . . . , Y n ) T ∈ R n a vector of responses, X ∈ R n×p a known design matrix, θ 0 ∈ R p a vector of parameters, and W ∈ R n random noise having zero-mean components W = (W 1 , . . . , W n ) T i.i.d. with distribution F = F W having finite second moment 1 . We are interested in estimating θ 0 from observed data 2 (Y, X) using a traditional M-estimator, defined by a non-negative convex function ρ : R → R ≥0 : 2 We denote by X1, . . . , Xn the rows of X. We often omit the arguments Y , X as this dependency will hold throughout. Without loss of generality, we assume that the columns of X are normalized so that X ei 2 ≈ 1. (A more precise assumption will be formulated below.) Although this is a completely traditional problem, we consider it under high-dimensional asymptotics where the number of parameters p and the number of observations n are both tending to infinity, at the same rate. This is becoming a popular asymptotic model owing to the modern awareness of 'big data' and 'data deluge'; but also because it leads to entirely new phenomena. Classical statistical theory considered the situation where the number of regression parameters p is fixed and the number of samples n is tending to infinity. The asymptotic distribution was found by Huber [Hub73, Bic75] to be normal N(0, V) where the asymptotic variance matrix V is given by here ψ = ρ is the score function of the M-estimator and V (ψ, F ) = ( ψ 2 dF )/( ψ dF ) 2 the asymptotic variance functional of [Hub64] , and (X T X) the usual Gram matrix associated with the leastsquares problem. Importantly, it was found that for efficient estimation -i.e. the smallest possible asymptotic variance -the optimal M-estimator depended on the probability distribution F W of the errors W . Choosing ψ(x) = (log f W (x)) (with f W the density of W ), the asymptotic variance functional yields V (ψ, F W ) = 1/I(F W ), with I(F ) denoting the Fisher information. This achieves the fundamental limit on the accuracy of M-estimators [Hub73] . In modern statistical practice there is increasing interest in applications where the number of explanatory variables p is very large, and comparable to n. Examples of this new regime can be given, spanning bioinformatics, machine learning, imaging, and signal processing (a few research areas in the last domains include [LDSP08, Sca97, Ric05, Cha03] ). This paper considers the properties of M-estimators in the high-dimensional asymptotic n → ∞, n/p(n) → δ ∈ (1, ∞) In this regime, the asymptotic distribution of M-estimators no longer needs to obey the classical formula (3) in widespread use. We make a random-design assumption on the X's detailed below. We show that the asymptotic covariance matrix of the parameters is now of the form where V is still Huber's asymptotic variance functional, butΨ is the effective score function, which is different from ψ under high-dimensional asymptotics andF W is the effective error distribution, which is different from F W under high-dimensional asymptotics. In the limit δ → ∞, the effective score and the effective error distribution both tend to their classical counterparts, and one recovers The effective error distributionF W is a convolution of the noise distribution with an extra Gaussian noise component, not seen in the classical setting (here denotes convolution): The extra Gaussian noise depends in a complex way on ψ, F W , δ, which we characterize fully below in Corollary 4.2. Several important insights follow immediately: 2 1. Existing formulas are inadequate for confidence statements about M-estimates under high dimensional asymptotics, and will need to be systematically broadened. 2. Classical maximum likelihood estimates are inefficient under high-dimensional asymptotics. The idea dominating theoretical statistics since R.A. Fisher to use ψ = (− log f W ) as a scoring rule, does not yield the efficient estimator. 3. The usual Fisher Information bound is not necessarily attainable in the high-dimensional asymptotic, as I( F W ) < I(F W ). M-estimation in this high-dimensional asymptotic setting was considered in a recent article by El Karoui, Bean, Bickel, Lim, and Yu [EKBBL13], who studied the distribution of θ for Gaussian design matrices X. In short they observed empirically the basic phenomenon of extra Gaussian noise appearing in high-dimensional asymptotics and rendering classical inference incorrect. The dependence of the additional variance τ 2 * on δ, ψ and F was characterized by [EKBBL13] through a non-rigorous heuristics 3 that the authors describe as 'highly plausible and buttressed by simulations.' 4 (We refer to Section 5 for further discussion of related work.) In the present paper, we show that this important statistical phenomenon can be characterized rigorously, in a way that we think fully explains the main new concepts of extra Gaussian noise, effective noise and the effective score. Our proof strategy has three steps • Introduce an Approximate Message Passing (AMP) algorithm for M-estimation; an iterative procedure with the M-estimator as a fixed point, and having the effective score functionΨ as its score function at algorithm convergence. • Introduce State Evolution for calculating properties of the AMP algorithm iteration by iteration. We show that these calculations are exact at each iteration in the large-n limit where we freeze the iteration number and let n → ∞. At the center of the State Evolution calculation is precisely an extra Gaussian noise term that is tracked from iteration to iteration, and which is shown to converge to a nonzero noise level. In this way, State Evolution makes very explicit that AMP faces at each iteration and even in the limit, an effective noise that differs from the noise W by addition of an appreciable extra independent Gaussian noise. • Show that the AMP algorithm converges to the solution of the M-estimation problem in mean square, from which it follows that the asymptotic variance of the M-estimator is identical to the asymptotic variance of the AMP algorithm. More specifically, the asymptotic variance of the M-estimator is given by a formula involving the effective score function and the effective noise. 3 To the reader familiar with the mathematical theory of spin glasses, the argument of [EKBBL13] appears analogous to the cavity method from statistical physics [MPV87, MM09, Tal10] 4 After the first version of our manuscript was posted on ArXiv, Noureddine El Karoui announced an independent proof of related results, using a completely different approach. As it turns out, our formula for the asymptotic variance coincides with the one derived heuristically in [EKBBL13, Corollary 1] although our technique is remarkably different, and our proof provides a very clear understanding of the operational significance of the terms appearing in the asymptotic variance. It also allows explicit calculation of many other operating characteristics of the M-estimator, for example when used as an outlier detector 5 . At the heart of our analysis, we are simply applying an approach developed in [BM11, BM12] for rigorous analysis of solutions to convex optimization problems under high-dimensional asymptotics. That approach grew out of a series of earlier papers studying the compressed sensing problem [DMM09, DMM11, DJMM11, BM12]. From the perspective of this paper, those papers considered the same regression model (1) as here; however, they emphasized the challenging asymptotic regime where there are fewer observations than predictors, (i.e. n/p(n) → δ ∈ (0, 1)) so that even in the noiseless case, the equations Y = Xθ would be underdetermined. In the p > n setting, it became popular to use 1 -penalized least squares (Lasso, [Tib96, CD95] ). That series of papers considered the Lasso convex optimization problem in the case of X with iid N(0, 1/n) entries (just as here) and followed the same 3-step strategy we use here; namely, 1. Introducing an AMP algorithm; 2. Obtaining the asymptotic distribution of AMP by State Evolution; and 3. Showing that AMP agrees with the Lasso solution in the large-n limit. This procedure proved that the Lasso solution has the asymptotic distribution where σ 2 is the variance of the noise in the measurements, and τ 2 Lasso is the variance of an extra Gaussian noise, not appearing in the classical setting where p(n)/n → 0. The variance of this extra Gaussian noise was obtained by state evolution and shown to depend on the distribution of the coefficients being recovered, and on the noise level in a seemingly complicated way that can be characterized by a fixed-point relation, see [DMM11, BM12] . At the center of the rigorous analysis stand the papers [BM11, BM12] which analyze recurrences of the type used by AMP and establish the validity of State Evolution in considerable generality. Those same papers stand at the center of our analysis in this paper. Apart from allowing a simple treatment, this provides a unified understanding of the phenomenon of high-dimensional extra Gaussian noise. This paper introduces a new first-order algorithm for computing the M-estimator θ which is uniquely appropriate for the random-design case. This algorithm fits within the class of approximate message passing (AMP) algorithms introduced in [DMM09, BM11] (see also [Ran11] for extensions). This algorithm is of independent interest because of its low computational complexity. AMP has a deceptive simplicity. As an iterative procedure for convex optimization, it looks almost the same as the 'standard' application of simple fixed-stepsize gradent descent. However, it is 5 The slightly more general [EKBBL13, Result 1] covers heteroscedastic noise is not covered by the analysis of this paper, but should be provable by adapting our argument. 4 intended for use in the random-design setting, and it has an extra memory term (aka reaction term) that modifies the iteration in a profound and beneficial way. In the Lasso setting, AMP algorithms have been shown to have remarkable fast convergence properties [DMM09], far outperforming more complex-looking iterations like Nesterov and FISTA. In the present paper, AMP has an second important wrinkle -it solves a convex optimization problem associated to minimizing ρ with iterations based on gradient descent with an objective ρ bt which varies from one iteration to the next, as b t changes, but which does not tend to ρ in the limit. In the present paper, AMP is mainly used as a proof device, one component of the three-part strategy outlined earlier. However, a key benefit produced by the curious features of AMP is strong heuristic insight, which would not be available for a 'standard' gradient-descent algorithm. The AMP proof strategy makes visible the extra Gaussian noise appearing in the M-estimator θ. Elementary considerations show that such extra noise is present at iteration zero of AMP. State Evolution faithfully tracks the dynamics of this extra noise across iterations. State Evolution proves that the extra noise level does not go to zero asymptotically with increasing iterations, but instead that the extra noise level tends to a fixed nonzero value. Because AMP is solving the M-estimation problem, the M-estimator must be infected by this extra noise. The AMP algorithm and its State Evolution analysis shows that the extra noise in parameter θ t i at iteration t is due to cross-parameter estimation noise leakage, where errors in the estimation of all other parameters at the previous iteration (t − 1) cause extra noise to appear in θ t i . In the classical setting no such effect is visible. One could say that the central fact about the high-dimensional setting revealed here as well as in our earlier work [DMM09, DMM11, DJMM11, BM12], is that when there are so many parameters to estimate, one cannot really insulate the estimation of any one parameter from the errors in estimation of all the other parameters. 2 Approximate Message Passing (AMP) For the rest of the paper, we make the following smoothness assumption on ρ: Definition 2.1. We call the loss function ρ : R → R smooth if it is continuously differentiable, with absolutely continuous derivative ψ = ρ having an a.e. derivative ψ that is bounded: sup u∈R ψ (u) < ∞. Our assumption excludes some interesting cases, such as ρ(u) = |u|, but includes for instance the Huber loss 6 Associated to ρ, we introduce the family ρ b of regularizations of ρ: 6 We expect that the proof technique developed in this paper should be generalizable to a broader class of functions ρ, at the cost of additional technical complications. 5 in words, this is the min-convolution of the original loss with a square loss. Each ρ b has a corresponding score function The effective score of the M-estimator belongs to this family, for a particular choice of b, explained below. In the classical M-estimation literature [HR09], monotonicity and differentiability of the score function ψ is frequently useful; our assumptions on ρ guarantee these properties for the nominal score function ψ. The score family Ψ( · ; b) has such properties as well: for any b, Ψ( · ; b) is a strictly monotone increasing function; second, for any b > 0, Ψ( · ; b) is a contraction. With Ψ denoting differentiation with respect to the first variable, we have Ψ (z; b) ∈ (0, 1). For proof and further discussion, see Appendix A. Before proceeding, we give an example. Consider the Huber loss ρ H (z; λ), with score function ψ(z; λ) = min(max(−λ, z), λ). We have In particular the shape of each Ψ is similar to ψ, but the slope of the central part is now Ψ ( Our proposed approximate message passing (AMP) algorithm for the optimization problem (2) is iterative, starting at iteration 0 with an initial estimate θ 0 ∈ R p . At iteration t = 0, 1, 2, . . . it applies a simple procedure to update its estimate θ t ∈ R p , producing θ t+1 . The procedure involves three steps. Adjusted residuals. Using the current estimate θ t , we compute the vector of adjusted residuals where to the ordinary residuals Y − X θ t we here add the extra term 7 Ψ(R t−1 ; b t−1 ). Effective Score. We choose a scalar b t > 0, so that the effective score Ψ( · ; b t ) has empirical average slope p/n ∈ (0, 1). Setting δ = δ(n) = n/p > 1, we take any solution 8 (for instance the smallest solution) to 9 : 7 Here and below, given f : 8 This equation always admits at least one solution since b → Ψ (r; b) is continuous in b ≥ 0, with Ψ (r; 0) = 0 and (for ρ strictly convex) Ψ (r; ∞) = 1, cf. Proposition A.1. 9 Under this prescription, the sequence bt depends on the instance (Y, X). As explained in the next section, for the proof of our main result we will use a slightly different prescription, that is independent of the problem instance. Scoring. We apply the effective score function Ψ(R t ; b t ): The Scoring step of the AMP iteration (11) is similar to traditional iterative methods for Mestimation, compare [Bic75] . Indeed, using the traditional residual z t = Y − Xθ t , the traditional method of scoring at iteration t would read and one can see correspondences of individual terms to the method of scoring used in AMP. Of course the traditional term [ (10)), while the traditional term (X T X) −1 corresponds to AMP's implicit I p×p -which is appropriate in the present context because our random-design assumption below makes X T X behave approximately like the identity matrix. The next lemma explains the reason for using the effective score Ψ(·; b t ) in the AMP algorithm: this is what connects the AMP iteration to M-estimation (2). Lemma 2.2. Let ( θ * , R * , b * ) be a fixed point of the AMP iteration (9), (10), (11) having b * > 0. Then θ * is a minimizer of the problem (2). Viceversa, any minimizer θ * of the problem (2) corresponds to one (or more) AMP fixed points of the form ( θ * , R * , b * ). Proof. By differentiating Eq. (2), and omitting the arguments Y, X for simplicity from L(θ; Y, X), we get where as usual ρ is applied component-wise to vector arguments. The minimizers of L(θ) are all the vectors θ for which the right hand side vanishes. Consider then a fixed point ( θ * , R * , b * ), of the AMP iteration (9), (11). This satisfies the equations The first equation can be written as Using Proposition A.2 below, (16) implies that Ψ(R * ; b * ) = b * ρ (Y −X θ * ). Hence the second equation reads which coincides with the stationarity condition (13) for b * > 0. This concludes the proof.  To make the AMP algorithm concrete, we consider an example with n = 1000, p = 200, so δ = 5. For design matrix we let X i,j ∼ N(0, 1 n ), and we draw θ 0 a random vector of norm θ 0 2 = 6 √ p. For the distribution F = F W of errors, we use Huber's contaminated normal distribution CN(0.05, 10), so that F = 0.95Φ + 0.05H 10 , where H x denotes a unit atom at x. For the loss function, we use the Huber's ρ H (z; λ) with λ = 3. Starting the AMP algorithm with θ 0 = 0, we run 20 iterations. Separately, we solved the M-estimation problem using CVX, obtaining θ. Figure 1 (left panel) shows the progress of the AMP algorithm across iterations, presenting while Figure 1 (right panel) shows the progress of AMP in approaching the M-estimate θ, as measured by As is evident, the iterations converge rapidly, and they converge to the M-estimator, both in the sense of convergence of risks -measured here by RMSE( θ t ; θ 0 ) → RMSE( θ; θ 0 ) ≈ 1.6182 -and, more directly, in convergence of the estimates themselves: RMSE( θ t ; θ) → 0. Figure 2 (left panel) shows the process by which the effective score parameterb t is obtained at iteration t = 3, while the right panel shows howb t behaves across iterations. In fact it converges quickly towards a limit b ∞ ≈ 0.2710. Earlier we pointed to resemblances between AMP (11) and the traditional method of scoring for obtaining M-estimators (12). In reality the two approaches are very different: • The precise form of various terms in (9), (10) (11) is dictated by the statistical assumptions that we are making on the design X. In particular the memory terms are crucial for the state evolution analysis to hold. Several papers document this point [Mon12, Sch10, SSS10, Ran11, KMZ13]. • Under classical asymptotics, where p is fixed and n → ∞, it is sufficient to run a single step of such an algorithm [Bic75], in the high-dimensional setting it is necessary to iterate numerous times. The resulting analysis is considerably more complex because of correlations arising as the algorithm evolves. 3 State evolution description of AMP State Evolution is a method for computing the operating characteristics of the AMP iterates θ t and R t for arbitrary fixed t, under the high-dimensional asymptotic limit n, p → ∞, n/p → δ. In this section we initially describe a purely formal procedure which assumes that the AMP adjusted residuals R t = Y − X θ t + Ψ(R t ; b t ) really behave as W + τ t Z, with W the error distribution and Z an independent standard normal, for t = 0, 1, 2, . . .  one to calculate τ 2 t for each t = 0, 1, 2, 3, . . . , and to calculate the limit of τ t as t → ∞. Later in the section we present a rigorous result validating the method under the following random Gaussian design assumption. Definition 3.1. We say that a sequence of random design matrices {X(n)} n , with n → ∞ is a Gaussian design if each X = X(n) has dimensions n × p, and entries (X ij ) i∈[n],j∈ [p] that are i.i.d. N(0, 1/n). Further, p = p(n) is such that lim n→∞ n/p(n) = δ ∈ (0, ∞). Under the Gaussian design assumption, suppose that u is a vector in R p with norm u 2 . Then {E Xu 2 2 } = u 2 2 . Moreover, Xu is a Gaussian random vector with entries iid N(0, u 2 2 /n). It will be convenient to introduce for any estimatorθ the notation So initialize AMP with a deterministic estimate θ 0 , and take R −1 = 0. Then the initial residual is R 1 = Y − X θ 0 = W + X(θ 0 − θ 0 ). The terms W and X(θ 0 − θ 0 ) are independent, and X(θ 0 − θ 0 ) is Gaussian with variance Hence, when AMP is started this way, we see that the adjusted residuals initially contain an extra Gaussian noise of variance τ 2 0 = MSE( θ 0 , θ 0 )/δ. Assuming the adjusted residuals continue, at later iterations, to behave as W + τ t Z with Z an independent standard normal, we now calculate τ 2 t for each t = 1, 2, 3, . . . , and eventually identify the limit of τ t as t → ∞. For a given τ > 0, δ = n/p and noise distribution F W , define the variance map where W ∼ F W , and, independently, Z ∼ N(0, 1). In this display, the reader can see that extra Gaussian noise of variance τ 2 is being added to the underlying noise W , and V measures the δ-scaled variance of the resulting output. Evidently for Under our assumptions for Ψ, for each given specification (τ ; δ, F W ) of the ingredients besides b that go into V, there is (as clarified by Lemma A.3) a well-defined value b = b(τ ; δ, F W ) giving the smallest solution b ≥ 0 to Definition 3.2. State Evolution is an iterative process for computing the scalars {τ 2 t } t≥0 , starting from an initial condition τ 2 0 ∈ R ≥0 following DefiningṼ(τ 2 ) = V(τ 2 , b(τ )), we see that the evolution of τ t follows the iterations of the mapṼ. In particular, we make these observations: •Ṽ(τ 2 ) is a continuous, nondecreasing function of τ . •Ṽ(τ 2 ) < τ 2 as τ → ∞. Figure 3 , left panel, considers the case where W again follows the Huber's contaminated normal distribution CN(0.05, 10) and ψ is the standard Huber estimator with parameter λ = 3. The ratio n/p = δ = 2, and the parameter vector has θ 0 2 2 /p = 6 2 . It displays the functionṼ (τ 2 ) as a function of τ . Evidently, there is a stable fixed point τ * = τ * (δ, F W ), i.e. a point obeyingṼ(τ 2 * ) = τ 2 * , such that τ 2 →Ṽ(τ 2 ) has a derivative less than 1 at τ 2 * . We conclude that τ t evolves under state evolution to a nonzero limit. Figure 3 , right panel, shows how τ 2 t evolves to the fixed point near 0.472 starting from τ 2 0 = 2.056. State Evolution offers a formal 10 procedure for predicting operating characteristics of the AMP iteration at any fixed iteration t or in the limit t → ∞. Nater in this section, we will provide rigorous validation of these predictions. 10 By formal, we mean a rule-based procedure which we can follow to get a prediction, without any guarantees that the prediction is correct. Call the tuple S = (τ ; b, δ, F ) a state; in running the AMP algorithm we assume that the algorithm is initialized with θ 0 so that τ 2 0 = MSE( θ 0 , θ 0 )/δ, so that AMP starts in state S = (τ 0 ; b 0 , δ, F ), and visits S 1 = (τ 1 ; b 1 , δ, F ), S 2 = (τ 2 ; b 2 , δ, F ), . . . ; eventually AMP visits states arbitrarily close to the equilibrium state S * = (τ * ; b * , δ, F ). SE predictions of operating characteristics are provided by two rules assigning predictions to certain classes of observables, based on the state that AMP is in. Definition 3.3. The state evolution formalism assigns predictions E to two types of observables under specific states. Observables Involving θ − θ 0 . Given a univariate test function ξ : R → R, assign the predicted value for p −1 i∈p ξ( θ i − θ 0,i ) under state S by the rule where expectation on the right hand side is with respect to Z ∼ N(0, 1). Observables involving Residual, Error. Let R denote some coordinate of the adjusted residual for AMP in state S and W the same coordinate of the underlying error. Given a bivariate test function ξ 2 : R 2 → R, assign the prediction of where Z ∼ N(0, 1) and W ∼ F W is independent of Z. The two most important predictions of operating characteristics are undoubtedly: • MSE at iteration t. We let S t = (τ t , b(τ t ), δ, F W ) denote the state of AMP at iteration t, and predict • MSE at convergence. With τ * > 0 the limit of τ t , let S * = (τ * , b(τ * ), δ, F W ) denote the state of AMP at convergence. and predict Other predictions might also be of interest. Thus, concerning the mean absolute error MAE( θ t , θ 0 ) = θ t − θ 0 1 /p, state evolution predicts MAE ≈ 2δτ 2 t /π. Concerning functions of (R, W ), consider the ordinary residuals Y − X θ * at AMP convergence. These residuals will of course in general not have the distribution of the errors W . Setting η(z; b) = z − Ψ(z; b), we have Y − X θ * = η(R; b * ). State evolution predicts that the ordinary residuals will have the same distribution as η(W +τ * Z; b * ). Continuing with our running example, we again consider the case of contaminated normal data W ∼ CN(0.05, 10) and Huber ρ with λ = 3. If we start AMP with the all-zero estimate θ 0 = 0, then since θ 0 2 = 6 √ p we start SE with τ 0 = 2.056. Again in our running example, these predictions can be tested empirically. For illustration, we conducted a very small experiment, generating 10 independent realizations of the running model at n = 1000 and p = 200, and comparing the actual evolutions of observables during AMP iterations with the predicted evolutions. Figure 5 shows that the predictions from SE are very close to the averages across realizations. State Evolution cannot evolve so that τ 2 t → 0; under minimal regularity, it always exceeds a specific nonzero noise level. Lemma 3.4. Suppose that F W has a well-defined Fisher information I(F W ). Then for any t > 0 . δ and, if ξ G denotes the score function for location of G, then .  From convexity and translation-invariance of Fisher Information We can sharpen this bound one step further. It will be convenient to write I(X) for the Fisher information of distribution F X . Lemma 3.5. Proof. Barron and Madiman [MB07] give the inequality I(W + τ Z) ≤ x 2 I(W ) + (1 − x) 2 I(τ Z), valid for any x ∈ (0, 1). By calculus, we know that for a, b > 0, Setting a = I(W ) and b = I(τ Z) = τ −2 , and dividing both numerator and denominator by b, we are done. Revisit the argument of Lemma 3.4; the inequality τ 2 t ≥ 1/(δI(F W )) shows that if t > 0, then τ 2 t I(W ) ≥ 1/δ. Using this in the previous Lemma, This yields a 'one-step' improvement: Corollary 3.6. Suppose that F W has a well-defined Fisher Information I(F W ). Then for any t > 1 We can iterate this argument across many steps, obtaining that, for every t > k, We obtain immediately: Corollary 3.7. Suppose that F W has a well-defined Fisher information I(F W ). Then for every accumulation point τ * of State Evolution . The predictions of state evolution can be validated in the large-system limit n, p → ∞, under the random Gaussian design assumption of Definition 3.1. We impose regularity conditions on the observables whose behavior we attempt to predict: In particular, ξ(x) = x 2 is pseudo-Lipschitz. Recall also the definition of MSE in equation (18). For a sequence of estimatorsθ, define the per-coordinate asymptotic mean squared error (AMSE) as the following large-system limit: when the indicated limit exists. The following result validates the predictions of State Evolution for pseudo-Lipschitz observables. Our proof is deferred to Appendix B. Theorem 3.9. Assume that the loss function ρ is convex and smooth, that the sequence of matrices {X(n)} n is a standard Gaussian design, and that θ 0 , θ 0 are deterministic sequences such that AMSE(θ 0 , θ 0 ) = δτ 2 0 . Further assume that F W has finite second moment and let {τ 2 t } t≥0 be the state evolution sequence with initial condition τ 2 0 . Let { θ t , R t } t≥0 be the AMP trajectory with parameters b t as per Eq. (19). Let ξ : R → R, ξ 2 : R × R → R be pseudo-Lipschitz functions. Then, for any t > 0, we have, for In particular, we may take ξ(x) = x 2 and obtain for the AMP iteration in full agreement with the predictions of state evolution in Definition 3.3. The key step for characterizing the distribution of the M-estimator θ, cf. Eq. (2), is to prove that the AMP iterates θ t converge to θ. We will prove that this is indeed the case, at least in the limit n, p → ∞, and for suitable initial conditions 11 . Throughout this section, we shall assume that ρ is strongly convex, i.e. that inf x∈R ρ (x) > 0. This corresponds to assuming inf x∈R ψ (x) > 0, which is rather natural from the point of view of robust statistics since it ensures uniqueness of the M estimator 12 . The key step is to establish the following high-dimensional convergence result. 11 We expect convergence for arbitrary initial conditions (as long as they are independent of (W, X)), but proving this claim is not needed for our main goal, and we leave it for future study. Proving this claim would require showing convergence of the state evolution recursion (20). 12 The Huber estimator is not covered by the result of this section; although we expect our approach to apply in such generality. We focus here on the strongly convex case to avoid un-necessary complications. Theorem 4.1. (Convergence of AMP to the M-Estimator.) Assume the same setting as in Theorem 3.9, and further assume that ρ is strongly convex and that δ > 1. Let (τ * , b * ) be a solution of the two equations and assume that AMSE( θ 0 , θ 0 ) = δτ 2 * . Then From this and Theorem 3.9, the desired characterization of θ immediately follows. To tie back to the introduction, we prove formula (4): Corollary 4.2. (Asymptotic Variance Formula under High-Dimensional Asymptotics.) Assume the setting of Theorem 3.9, and further assume that ρ is strongly convex and δ > 1. The asymptotic variance of θ obeys lim n,p→∞ where Ave i∈[p] denotes the average across indices i, V (ψ, F ) denotes the usual Huber asymptotic variance formula for M-estimates -V (ψ, F ) = ( ψ 2 dF )/( ψ dF ) 2 -and the effective scoreΨ is while the effective noise distributionF isF Here (τ * , b * ) are the unique solutions of the equations (24)-(25). Proof. By symmetry, Ave i∈ [p] Var( θ i ) = EMSE( θ, θ 0 ). Theorem 4.1 and State Evolution show that AMSE( θ, θ 0 ) = δτ 2 * . By (24)- (25) V Recall that the traditional information bound for M-estimators is V (Ψ, F ) ≥ 1 I(F W ) , and that this is achievable under p fixed, n → ∞ asymptotics. Considering the formula forF we see that because τ * > 0, such an asymptotic variance is not achievable under high-dimensional asymptotics. We now make this effect more visible. Combining Corollary 4.2 with Corollary 3.7's lower bound on the equilibrium noise τ * reachable by State Evolution, we have the following.  . In this inequality, the effect of the high-dimensional asymptotics parameter δ is extremely clear; it shows that the classical information bound is not achievable when δ = n/p < ∞, There is always an inflation in variance at least by (1 − 1/δ) −1 = n n−p . Moreover, the inflation completely blows up as δ → 1. Corollary 4.4. Assume the setting of Theorem 3.9, and further assume that ρ is strongly convex and δ > 1. Then for any pseudo-Lipschitz function ξ : R → R, we have, for Z ∼ N(0, 1) In particular, the solution of Eqs. (24), (25) is necessarily unique. Among other applications, this result can be used to bound the suboptimality of AMP after a fixed number of iterations. Combining Theorems 3.9 and 4.1 gives: Corollary 4.5. Assume the same setting as in Theorem 3.9, and further assume that ρ is strongly convex and δ > 1. Then the almost sure limits AMSE( θ t ; θ 0 ) and AMSE( θ; θ 0 ) exist, and obey Theorem 4.1 extends to cover general Gaussian matrices X with i.i.d. rows. Definition 4.6. We say that a sequence of random design matrices {X(n)} n , with n → ∞, is a general Gaussian design if each X = X(n) has dimensions n × p, and rows (X i ) i∈ [n] that are i.i.d. N(0, Σ/n), where Σ = Σ(n) ∈ R p×p is a strictly positive definite matrix. Further, p = p(n) is such that lim n→∞ n/p(n) = δ ∈ (0, ∞). Notice that, if X is a general Gaussian design, then XΣ −1/2 is a standard Gaussian design. The following then follows from Corollary 4.7 together with a simple change of variables argument, cf. [EKBBL13, Lemma 1]. Corollary 4.7. Assume the same setting as in Theorem 3.9, but with {X(n)} n≥0 being a general Gaussian design with covariance Σ, and further assume that ρ is strongly convex and δ > 1. There is a scalar random variable T n so that where Z ∼ N(0, I p×p ) and we have the almost-sure limit lim n→∞ T n = a.s. τ * , where τ * solves Eqs. (24), (25). This result coincides with Corollary 1 in [EKBBL13] apart from a factor √ n in the random part of Eq. (30) that arises because of a difference in the normalization of X. Several generalizations of the present proof technique should be possible, and would be of interest. We list a few in order of increasing difficulty: 19 1. Generalize the i.i.d. Gaussian rows model for X by allowing different rows to be randomly scaled copies of a common X ∼ N(0, Σ/n). This is the setting of [EKBBL13, Result 1]. 2. Remove the smoothness and strong convexity assumptions on ρ. 3. Add a regularization term to the objective function L(θ) cf. Eq. (2), of the form p i=1 J(θ i ), with J : R → R a convex penalty. For 1 penalty and 2 loss, this reduces to the Lasso, studied in [BM12]. 4. Generalize the present results to non-Gaussian designs. We expect -for instance-that they should hold universally across matrices X with i.i.d. entries (under suitable moment conditions). A similar universality result was established in [BLM12] for compressed sensing. Let us mention that alternative proof techniques would be worth exploring as well. In particular, Shcherbina and Tirozzi [ST03] define a statistical mechanics model with energy function that is analogous to the loss L(θ), cf. Eq. (2), and Talagrand [Tal10, Chapter 3] proves further results on the same model. While this treatment focuses on estimating a certain partition function, in the case of strongly convex ρ it should be possible to extract properties of the minimizer from a 'zero-temperature' limit. Finally, Rangan [Ran11] considers a similar regression model to the one studied here using approximate message passing algorithms, albeit from a Bayesian point of view. The reader might have noticed many analogies between the analysis in the last pages and earlier work on estimation in the underdetermined regime n < p using the Lasso [DMM09, DMM11, DJMM11, BM12]. Most specifically, the central tool in our proof of the correctness of State Evolution is a set of lemmas and theorems about analysis of recursive systems that were developed to understand the Lasso. That the same machinery directly gives results in robust regression -see for example our proof of correctness of State Evolution in Appendix B below -might seem particularly unexpected. In this section we briefly point out that the two problems are so closely linked that phenomena which appear in one situation are bound to appear in the other. In a very strong sense, solving an M-estimation problem with p < n is the very same thing as solving a related penalized regression problem inp >ñ. Given a convex function J : R → R, define the ρ function We then have the M-Estimation problem This problem has p < n and is generically a determined problem. We now construct a corresponding underdetermined problem with the 'same' solution. Setñ = n − p,p = n. We soon will 20 construct a vector/matrix pair ( Y ∈ Rñ, X ∈ Rñ ×p ) obeyingñ <p, where Y and X are related to Y and X in a specific way. With this pair we pose the J-penalized least squares problem with solution β( Y ; X), say. Here is the specific pair that links (M J ) with (L J ). We let X be a matrix with orthonormal rows such that XX = 0, i.e. finally, we set Y = XY . Of special interest is the case J(x) = λ |x| in which case (L J ) of (33) defines the Lasso estimator. Then ρ J (x) = ρ H (x; λ) is the Huber loss and (M J ) of (32) defines the Huber M-estimate. Indeed, in that case (L J ) is more classically presented as while (M J ) is more classically presented as (Huber λ ) min In this special case, our general result from the next section implies the following: Proposition 6.1. With problem instances (Y, X) and ( Y , X) related as above, the optimal values of the Lasso problem (Lasso λ ) and the Huber problem (Huber λ ) are identical. The solutions of the two problems are in one-one-relation. In particular, we have In a sense the Lasso problem solution β is finding the outliers in Y ; once the solution is known, the solution of the M-estimation problem is simply a least squares regression on adjusted data Y adj ≡ (Y −β) with outliers removed. We will now show that the problem (32) is dual to (33) under or special choice of ( Y , X), via (34). Notation. For x ∈ R n , we denote by ∂ρ(x) the subgradient of the convex function n i=1 ρ(x i ), at x. Analogously, for z ∈ Rp, we denote by ∂J(z) the subgradient of the convex function p i=1 J(z i ), at z. Proposition 6.2. Assume that ρ( · ) = ρ J ( · ), that X has orthonormal rows with null( X) = image(X), and finally that Y = X Y . Then the solutions of the regularized least squares problem (33) are in one-to-one correspondence with the solutions of the robust regression problem (2), via the mappings Proof. 'Differentiating' Eq. (31) it is easy to see that u ∈ ∂ρ(x) if and only if u ∈ ∂J(x − u) . First assume θ is a minimizer of problem (32). This happens if and only if there exists u ∈ R n such that We then claim that β ≡ Y − X θ − u is a minimizer of Eq. (33). Indeed where the last identity follows since, by Eq. (34), null(X T ) = image( X T ), and hence u ∈ image( X T ) by Eq. (41). Using again Eqs. (41) and (40), we deduce that u ∈ ∂J( β), i.e. which is the stationarity condition for the problem (33). Viceversa a similar argument shows that, given β that minimizes Eq. (33), and θ ≡ (X T X) −1 X T (Y − β) is a minimizer of the robust regression problem (32). The last section raises the possibility that the phenomena found in this paper for M-estimation in the p < n case are actually isomorphic to those found in our previous work on penalized regression in the p > n case; [DMM09, DMM11, DJMM11, BM12]. Here we merely content ourselves with sketching a few similarities. To be definite, consider robust regression using the Huber loss [Hub64, HR09] ρ(x) = x 2 /2 for |x| ≤ λ and ρ(x) = λ|x| − λ 2 /2 otherwise. In this case it is easy to see that  In order to make contact with the Lasso, recall the definition of soft thresholding operator η(x; α) = sign(x) (|x| − α) + . We have the relationship 22 Letting c t ≡ b t /(1 + b t ), the state evolution equation (20), then reads This is very close to the state evolution equation in compressed sensing for reconstructing a sparse signal whose entries have distribution F W , from an underdetermined number of linear measurements; indeed in that setting we have the state evolution recursion [DMM09, DMM11, DJMM11, BM12]. The connection is quite suggestive: while in compressed sensing we look for the few non-zero coefficients in the signal, in robust regression we try to identify the few outliers contaminating the linear relation. A similar duality was already pointed out in [DT09], although in a specific setting.
A closely related work REF studied the same questions using a different technique known as approximate message passing (AMP), again assuming no regularization.
High Dimensional Robust M-Estimation: Asymptotic Variance via Approximate Message Passing
ArXiv
ArXiv
Mathematics Computer Science