High Dimensional Robust M-Estimation: Asymptotic Variance via Approximate Message Passing

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 techniques that were 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


M-Estimation under high dimensional asymptotics
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 : where u, v = m i=1 u i v i is the standard scalar product in R m , and θ is chosen arbitrarily if there is multiple minimizers.
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.

Extra Gaussian noise due to high-dimensional asymptotics
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 (3) 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 and FW 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 V (ψ, F W ). The effective error distribution FW is a convolution of the noise distribution with an extra Gaussian noise component, not seen in the classical setting (here denotes convolution): 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 heuristics3 that the authors describe as 'highly plausible and buttressed by simulations.'4(We refer to Section 5 for further discussion of related work.)

Proof Strategy: Approximate Message Passing
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.
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 detector5 .

Underlying tools
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.

The role of AMP
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 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 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.

Approximate Message Passing (AMP) 2.1 A family of score functions
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 Associated to ρ, we introduce the family ρ b of regularizations of ρ: 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 Ψ (

AMP algorithm
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 term7 Ψ(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 solution8 (for instance the smallest solution) to9 : 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.

Relation to M-estimation
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).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 which coincides with the stationarity condition (13) for b * > 0. This concludes the proof.

Example
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 parameter bt is obtained at iteration t = 3, while the right panel shows how bt behaves across iterations.In fact it converges quickly towards a limit b ∞ ≈ 0.2710.

Contrast to iterative M-estimation
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.

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, . . . .The variable τ 2 t thus quantifies the extra Gaussian noise supposedly present in the adjusted residuals of AMP; we show how this ansatz allows 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

Initialization of the extra variance
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 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 )/δ.

Evolution of the extra Gaussian variance to its ultimate limit
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 b > 0, 0 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 τ .
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 τ .

Predicting operating characteristics from State Evolution
State Evolution offers a formal10 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.
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( • 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 State evolution predicts that the ordinary residuals will have the same distribution as η(W +τ * Z; b * ).

Example of State Evolution predictions
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.

A lower bound on State Evolution
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 .
From convexity and translation-invariance of Fisher Information I(G) = I(F W N(0, τ 2 )) < I(F W ). Then τ 2 t = Ṽ(τ 2 t−1 ) ≥ 1/(δI(F W )). 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.
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

Correctness of State Evolution predictions
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.

Convergence and characterization of M-estimators
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 conditions11 .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.
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.
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 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 distribution F is 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) Recall that the traditional information bound for M-estimators is , and that this is achievable under p fixed, n → ∞ asymptotics.Considering the formula for F 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.
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.

Discussion
Several generalizations of the present proof technique should be possible, and would be of interest.We list a few in order of increasing difficulty: 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.

Duality between robust regression and regularized least squares
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.

Duality of optimization problems
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 in p > ñ.Given a convex function J : R → R, define the ρ function We then have the M-Estimation problem (M J ) min 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 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 .

The Lasso-Huber connection
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.

General duality result
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 ∈ R p, we denote by ∂J(z) the subgradient of the convex function p i=1 J(z i ), at z. 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 ].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.
B Proof of correctness of State Evolution (Theorem 3.9) We will show correctness of State Evolution for the AMP algorithm using analytically defined b t .Namely, we suppose that with b t defined recursively as the smallest positive solution of the second equation in this system: For analysis purposes, we consider a recursion equivalent to the AMP recursion, in which the data are recentered and the recursion is recast around recentered variables.We change the initial condition of the AMP iteration by letting θ cen,0 = θ 0 − θ 0 , and change data by letting Y cen = Y − Xθ 0 ≡ W . Applying the AMP recursion in these new coordinates gives the new trajectory θ cen,t = θ t − θ 0 for all t, and R cen,t = R t for all t.
The new trajectory follows the recursion In this form, the recursion can be reduced to a recursion studied in [BM11], for which State Evolution has been proven correct.The reduction is to introduce a recursion generating iterates {ϑ t , S t } that approximates closely the iterates { θ cen,t , R cen,t } defined by ( 64),(65).The new sequence is defined by letting ϑ 0 = θ 0 − θ 0 and, for all t ≥ 0 where The only difference between this recursion and the previous one cf.Eqs.(64), (65), lies in the new coefficient q t , which was identically equal to 1 in the previous recursion.The benefit of this specific recursion is that we already know that State Evolution is correct.
Lemma B.1.Under the assumptions of Theorem 3.9, we have, for any fixed t ≥ 0, Proof.This is an immediate application of Theorem 2 in [BM11].That Theorem considers general recursions which include (66)-(67) as a special case, and corresponding state evolution equations, and shows the correctness of state evolution, in the process establishing conclusions of the precise form shown in the conclusion of this lemma.So it is simply a matter of establishing the correspondence of variables.
In the original notation of [BM11], the generalized AMP recursions studied are b t = Aq t − λ t m t−1 (71) where b t , h t , q t and m t are vectors and λ t and ξ t scalars.In addition, the vectors q t = f t (h t ) and m t = g t (b t , w) are produced by element wise applications of nonlinearities f t and g t , the latter involving the random vector w.Here A is a rectangular n × N random matrix with iid Gaussian entries.The scalars ξ t = g t (b t , w) and λ t = 1 δ f t (h t ) , where • denotes an empirical mean over the entries in a vector.and the iteration takes m −1 = 0.For state evolution, the Theorem 2 assumes the sequence of initial conditions q 0 obeys σ 2 0 = lim and the state evolution recursion involves the pair of variables Table 1 sets up a 'dictionary' of correspondences between this paper and [BM11]. (66) Table 1: Correspondences between terms in the recursions of this paper, (66)-(67), and the recursions (71)-(72), analyzed in [BM11].
We get exact correspondence between the two systems, provided we identify δΨ(W + S t ; b t ) with m t = g t (b t ; w) and −δh t with f t (h t ).One has, in particular, that λ t = 1 δ f t (h t ) = −1, and that In the [BM11] general study of state evolution, there are two state variables τ t and σ t .However, when applied here, the distinction vanishes.The variable called τ 2 t corresponds with δ 2 E{Ψ(W + σ t Z) 2 } while the variable σ 2 t corresponds with E(ϑ t ) 2 .Using facts about τ 2 t in this paper, we have the identities σ 2 t = δτ 2 t and τ 2 t = δτ 2 t .Equations ( 69)-(70) now follow from Theorem 2 of [BM11].
Lemma B.2.Under the assumptions of Theorem 3.9, we have, for any fixed t ≥ 0,

B.1 Proof of Lemma B.2 (Equivalence of recursions)
Throughout this proof, we will drop the superscript 'cen' from R cen,t and θ cen,t .Define S t + ≡ W + S t , whence Comparing the first of these equations with Eq. (64), and using triangular inequality, we get where the last inequality follows since Ψ( • ; b) : R → R is Lipschitz continuous with Lipschitz constant at most 1, cf Proposition A.2. Comparing analogously Eq. ( 65) and (75), we obtain Iterating the upper bounds (77), (79), and using the fact that ϑ 0 = θ 0 , we conclude that there exists By Lemma B.1, we have, almost surely and where the second identity follows from Lemma B.1 and, in the third, we used the definition of b t .
(Note that we are applying here Lemma B.1 to ξ( • ) = Ψ ( • ; b t ) which is bounded and non-negative but not necessarily continuous.However, since W + τ t Z has a density for every τ t > 0, the limit holds by a standard weak convergence argument, approximating ξ by simple functions.Namely, we construct a sequence of simple functions ξ such that ξ (t) ≤ ξ(t) ≤ ξ (t) + (1/ ) for all t, and apply Lemma B.1 -which implies weak convergence of the empirical distribution of {W i + S t i }-to ξ .)Finally, it is a standard result in random matrix theory [AGZ09] that lim n→∞ X 2 = C(δ) < ∞.Hence, by taking the limit of Eq. (80) we get, almost surely, The norm R t − S t + 2 is then controlled using Eq.(77).
We are now ready to prove Theorem 4.1.Recall that L(θ) = L(θ; Y, X) denotes the loss function defined in Eq. (2), and that its gradient and Hessian are given by In particular, letting σ min (X) denote the minimum non-zero singular value of X, we have Using the hypothesis of strong convexity and standard concentration of measure for the singular values of Wishart matrices [Ver12], these exists constants c 0 , c 1 , n 0 > 0 for δ > 1 such that for any n ≥ n 0 , As a consequence, with probability at least 1 − e −c 1 n , we have Hence using Cauchy-Schwartz The last step of the proof consists in showing that, almost surely In order to prove this claim, reconsider Eq. ( 9), for time t + 1, with b t = b * .Using the fact that Ψ(z; b * ) = z − Prox(z; b * ), this can be rewritten as By Eq. ( 11), and recalling that Ψ(z; b) = b ρ (Prox(z; b)), we have where the last identity followed by Eq. (100).Using the triangle inequality and noting that, by the smoothness assumption C ≡ sup z∈R ρ (z) < ∞, we get Note that a similar statement is proved in [BM12, Theorem 4.2] for characterizing the Lasso estimator.While the same argument can be followed here, we outline an alternative argument that is based on a reduction to the setting of [JM12].
The matrix multiplication in Eq. ( 109) operates in the natural way over (R q ) N , namely we identified A with the Kronecker product A ⊗ I q×q .Explicitly, Eq. (109) reads A ij f j (z t j ; t) − B t f i (z t−1 i ; t − 1) .
Finally B t ∈ R q×q is given by The recursion (109) is characterized in [JM12, Theorem 1], which establishes -for instance-that, for ξ : R q → R pseudo-Lipschitz, we have, almost surely, lim Here Z t is a Gaussian random vector whose covariance is fully specified in [JM12].The proof of the lemma is finished by comparing the expressions in [JM12] for the covariance wit the ones in the statement of the lemma.

C.2 Proof of Lemma C.2
First of all we introduce the notation q t ≡ Γ t,t+1 /τ 2 * .We then have the recursion where expectation E q is with respect to the centered Gaussian vector (Z 1 , Z 2 ) with E q {Z 2 1 } = E q {Z 2 2 } = 1 and E q {Z 1 Z 2 } = q, independent of W ∼ F W .We claim that: (i) H(1) = 1; (ii) H(q) is increasing for q ∈ [0, 1]; (iii) H(q) is strictly convex for q ∈ [0, 1].In order to prove (i), note that, for q = 1, Z 1 = Z 2 ≡ Z ∼ N(0, 1) and hence which is equal to 1 since b * , τ * satisfy Eq. (24).In order to prove (ii), (iii), define H(q) ≡ E q {h W (Z 1 )h W (Z 2 )|W } , We will prove that H is strictly increasing and convex for any W , whence claims (ii) and (iii) follow by linearity.The argument is the same as in [BM12, Lemma C.1] Let {X t } t≥0 be the stationary Ornstein-Uhlenbeck process with covariance E(X 0 X t ) = e −t , and denote by E expectation with respect to X. Then Then we have the spectral representation (for t = log(1/q)) whence the claim follows since c = 0 for some ≥ 2 as long as h W (x) is non-linear.Because of the remarks (i)-(iii) just proven, it follows that lim t→∞ q t = 1 (and hence lim t→∞ Γ t,t+1 = τ 2 * ) if and only if H (1) ≤ 1.A simple calculation yields where Z ∼ N(0, 1).Recalling that Ψ (z; b) ∈ (0, 1), we have (Ψ ) 2 ≤ Ψ and so where the last identity follows because (τ * , b * ) solve Eq. ( 25).This finishes the proof.

CFigure 2 :
Figure 1: Left Panel: RMSE of AMP versus iteration (black curve), and its convergence to RMSE of M-estimation (constant green curve).Right Panel: Discrepancy of AMP from M-estimate, versus iteration.
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 .
C.1 Proof of Lemma C.1First of all note that, due to Lemma B.2, it is sufficient to prove that Eξ 2 (Z t , Z s ) .