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

M-Estimation under high dimensional asymptotics

Consider the traditional linear regression model

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 ${\sf N}(0,{\mathbf{V}})$ where the asymptotic variance matrix ${\mathbf{V}}$ is given by

here $\psi=\rho^{\prime}$ is the score function of the M-estimator and $V(\psi,F)=(\int\psi^{2}{\rm d}F)/(\int\psi^{\prime}{\rm d}F)^{2}$ the asymptotic variance functional of [Hub64], and $({\mathbf{X}}^{{\sf T}}{\mathbf{X}})$ the usual Gram matrix associated with the least-squares 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 $\psi(x)=(\log f_{W}(x))^{\prime}$ (with $f_{W}$ the density of $W$), the asymptotic variance functional yields $V(\psi,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\to\infty$, $n/p(n)\to\delta\in(1,\infty)$ 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 ${\mathbf{X}}$’s detailed below. We show that the asymptotic covariance matrix of the parameters is now of the form

The extra Gaussian noise depends in a complex way on $\psi$, $F_{W}$, $\delta$, which we characterize fully below in Corollary 4.2.

Several important insights follow immediately:

Existing formulas are inadequate for confidence statements about M-estimates under high dimensional asymptotics, and will need to be systematically broadened.

Classical maximum likelihood estimates are inefficient under high-dimensional asymptotics. The idea dominating theoretical statistics since R.A. Fisher to use $\psi=(-\log f_{W})^{\prime}$ as a scoring rule, does not yield the efficient estimator.

The usual Fisher Information bound is not necessarily attainable in the high-dimensional asymptotic, as $I(\widetilde{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 $\widehat{\theta}$ for Gaussian design matrices ${\mathbf{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 $\tau_{*}^{2}$ on $\delta$, $\psi$ and $F$ was characterized by [EKBBL13] through a non-rigorous heuristics 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] that the authors describe as ‘highly plausible and buttressed by simulations.’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. (We refer to Section 5 for further discussion of related work.)

2 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 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\to\infty$.

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 detectorThe 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..

3 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.

where $\sigma^{2}$ is the variance of the noise in the measurements, and $\tau_{\rm Lasso}^{2}$ is the variance of an extra Gaussian noise, not appearing in the classical setting where $p(n)/n\to 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.

4 The role of AMP

This paper introduces a new first-order algorithm for computing the M-estimator $\widehat{\theta}$ 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 $\rho$ with iterations based on gradient descent with an objective $\rho_{b_{t}}$ which varies from one iteration to the next, as $b_{t}$ changes, but which does not tend to $\rho$ 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 $\widehat{\theta}$. 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 $\widehat{\theta}_{i}^{t}$ 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 $\widehat{\theta}_{i}^{t}$. 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)

For the rest of the paper, we make the following smoothness assumption on $\rho$:

Our assumption excludes some interesting cases, such as $\rho(u)=|u|$, but includes for instance the Huber loss We expect that the proof technique developed in this paper should be generalizable to a broader class of functions $\rho$, at the cost of additional technical complications.

Associated to $\rho$, we introduce the family $\rho_{b}$ of regularizations of $\rho$:

in words, this is the min-convolution of the original loss with a square loss. Each $\rho_{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 $\psi$ is frequently useful; our assumptions on $\rho$ guarantee these properties for the nominal score function $\psi$. The score family $\Psi(\,\cdot\,;b)$ has such properties as well: for any $b$, $\Psi(\,\cdot\,;b)$ is a strictly monotone increasing function; second, for any $b>0$, $\Psi(\,\cdot\,;b)$ is a contraction. With $\Psi^{\prime}$ denoting differentiation with respect to the first variable, we have $\Psi^{\prime}(z;b)\in(0,1)$. For proof and further discussion, see Appendix A.

Before proceeding, we give an example. Consider the Huber loss $\rho_{{\rm H}}(z;\lambda)$, with score function $\psi(z;\lambda)=\min(\max(-\lambda,z),\lambda)$. We have

In particular the shape of each $\Psi$ is similar to $\psi$, but the slope of the central part is now $\|\Psi^{\prime}(\,\cdot\,;b)\|_{\infty}=\frac{b}{1+b}<1$.

2 AMP algorithm

We choose a scalar $b_{t}>0$, so that the effective score $\Psi(\,\cdot\,;b_{t})$ has empirical average slope $p/n\in(0,1)$. Setting $\delta=\delta(n)=n/p>1$, we take any solutionThis equation always admits at least one solution since $b\mapsto\Psi^{\prime}(r;b)$ is continuous in $b\geq 0$, with $\Psi^{\prime}(r;0)=0$ and (for $\rho$ strictly convex) $\Psi^{\prime}(r;\infty)=1$, cf. Proposition A.1. (for instance the smallest solution) to Under this prescription, the sequence $b_{t}$ depends on the instance $(Y,{\mathbf{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.:

We apply the effective score function $\Psi(R^{t};b_{t})$:

The Scoring step of the AMP iteration (11) is similar to traditional iterative methods for M-estimation, compare [Bic75]. Indeed, using the traditional residual $z^{t}=Y-{\mathbf{X}}\theta^{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 $[\sum_{i=1}^{n}\psi^{\prime}(z^{t}_{i})/n]^{-1}$ corresponds to AMP’s $[\sum_{i=1}^{n}\Psi^{\prime}(R^{t}_{i};b_{t})/n]^{-1}\equiv\delta$ (because of step (10)), while the traditional term $({\mathbf{X}}^{{\sf T}}{\mathbf{X}})^{-1}$ corresponds to AMP’s implicit ${\rm I}_{p\times p}$ – which is appropriate in the present context because our random-design assumption below makes ${\mathbf{X}}^{{\sf T}}{\mathbf{X}}$ behave approximately like the identity matrix.

3 Relation to M-estimation

The next lemma explains the reason for using the effective score $\Psi(\cdot;b_{t})$ in the AMP algorithm: this is what connects the AMP iteration to M-estimation (2).

Let $(\widehat{\theta}_{*},R_{*},b_{*})$ be a fixed point of the AMP iteration (9), (10), (11) having $b_{*}>0$. Then $\widehat{\theta}_{*}$ is a minimizer of the problem (2). Viceversa, any minimizer $\widehat{\theta}_{*}$ of the problem (2) corresponds to one (or more) AMP fixed points of the form $(\widehat{\theta}_{*},R_{*},b_{*})$.

By differentiating Eq. (2), and omitting the arguments $Y,{\mathbf{X}}$ for simplicity from ${\cal L}(\theta;Y,{\mathbf{X}})$, we get

where as usual $\rho^{\prime}$ is applied component-wise to vector arguments. The minimizers of ${\cal L}(\theta)$ are all the vectors $\theta$ for which the right hand side vanishes.

Consider then a fixed point $(\widehat{\theta}_{*},R_{*},b_{*})$, of the AMP iteration (9), (11). This satisfies the equations

Using Proposition A.2 below, (16) implies that $\Psi(R_{*};b_{*})=b_{*}\rho^{\prime}(Y-{\mathbf{X}}\widehat{\theta}_{*})$. Hence the second equation reads

which coincides with the stationarity condition (13) for $b_{*}>0$. This concludes the proof. ∎

4 Example

To make the AMP algorithm concrete, we consider an example with $n=1000$, $p=200$, so $\delta=5$. For design matrix we let $X_{i,j}\sim{\sf N}(0,\frac{1}{n})$, and we draw $\theta_{0}$ a random vector of norm $\|\theta_{0}\|_{2}=6\sqrt{p}$. For the distribution $F=F_{W}$ of errors, we use Huber’s contaminated normal distribution ${\sf CN}(0.05,10)$, so that $F=0.95\Phi+0.05H_{10}$, where $H_{x}$ denotes a unit atom at $x$. For the loss function, we use the Huber’s $\rho_{{\rm H}}(z;\lambda)$ with $\lambda=3$. Starting the AMP algorithm with $\widehat{\theta}^{0}=0$, we run 20 iterations.

Separately, we solved the M-estimation problem using CVX, obtaining $\widehat{\theta}$.

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 $\widehat{\theta}$, 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 ${\rm RMSE}(\widehat{\theta}^{t};\theta_{0})\to{\rm RMSE}(\widehat{\theta};\theta_{0})\approx 1.6182$ - and, more directly, in convergence of the estimates themselves: ${\rm RMSE}(\widehat{\theta}^{t};\widehat{\theta})\to 0$.

Figure 2 (left panel) shows the process by which the effective score parameter $\hat{b}_{t}$ is obtained at iteration $t=3$, while the right panel shows how $\hat{b}_{t}$ behaves across iterations. In fact it converges quickly towards a limit $b_{\infty}\approx 0.2710$.

5 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 ${\mathbf{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\to\infty$, 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 $\widehat{\theta}^{t}$ and $R^{t}$ for arbitrary fixed $t$, under the high-dimensional asymptotic limit $n,p\to\infty$, $n/p\to\delta$.

In this section we initially describe a purely formal procedure which assumes that the AMP adjusted residuals $R^{t}=Y-{\mathbf{X}}\widehat{\theta}^{t}+\Psi(R^{t};b_{t})$ really behave as $W+\tau_{t}Z$, with $W$ the error distribution and $Z$ an independent standard normal, for $t=0,1,2,\dots$. The variable $\tau_{t}^{2}$ thus quantifies the extra Gaussian noise supposedly present in the adjusted residuals of AMP; we show how this ansatz allows one to calculate $\tau_{t}^{2}$ for each $t=0,1,2,3,\dots$, and to calculate the limit of $\tau_{t}$ as $t\to\infty$. Later in the section we present a rigorous result validating the method under the following random Gaussian design assumption.

We say that a sequence of random design matrices $\{{\mathbf{X}}(n)\}_{n}$, with $n\to\infty$ is a Gaussian design if each ${\mathbf{X}}={\mathbf{X}}(n)$ has dimensions $n\times p$, and entries $(X_{ij})_{i\in[n],j\in[p]}$ that are i.i.d. ${\sf N}(0,1/n)$. Further, $p=p(n)$ is such that $\lim_{n\to\infty}n/p(n)=\delta\in(0,\infty)$.

So initialize AMP with a deterministic estimate $\widehat{\theta}^{0}$, and take $R^{-1}=0$. Then the initial residual is $R^{1}=Y-{\mathbf{X}}\widehat{\theta}^{0}=W+{\mathbf{X}}(\theta_{0}-\widehat{\theta}^{0})$. The terms $W$ and ${\mathbf{X}}(\theta_{0}-\widehat{\theta}^{0})$ are independent, and ${\mathbf{X}}(\theta_{0}-\widehat{\theta}^{0})$ is Gaussian with variance $\tau_{0}^{2}=\|\widehat{\theta}^{0}-\theta_{0}\|_{2}^{2}/n={\rm MSE}(\widehat{\theta}^{0},\theta_{0})/\delta$. Consider some fixed coordinate $R^{1}(i)$ of $R^{1}$. Then

Hence, when AMP is started this way, we see that the adjusted residuals initially contain an extra Gaussian noise of variance $\tau_{0}^{2}={\rm MSE}(\widehat{\theta}^{0},\theta_{0})/\delta$.

2 Evolution of the extra Gaussian variance to its ultimate limit

Assuming the adjusted residuals continue, at later iterations, to behave as $W+\tau_{t}\,Z$ with $Z$ an independent standard normal, we now calculate $\tau_{t}^{2}$ for each $t=1,2,3,\dots$, and eventually identify the limit of $\tau_{t}$ as $t\to\infty$.

For a given $\tau>0$, $\delta=n/p$ and noise distribution $F_{W}$, define the variance map

where $W\sim F_{W}$, and, independently, $Z\sim{\sf N}(0,1)$. In this display, the reader can see that extra Gaussian noise of variance $\tau^{2}$ is being added to the underlying noise $W$, and ${\cal V}$ measures the $\delta$-scaled variance of the resulting output. Evidently for $b>0$, $0\leq{\cal V}(\tau^{2},b)\cdot\delta\leq({\rm Var}(W)+\tau^{2})\cdot\delta$.

Under our assumptions for $\Psi$, for each given specification $(\tau;\delta,F_{W})$ of the ingredients besides $b$ that go into ${\cal V}$, there is (as clarified by Lemma A.3) a well-defined value $b=b(\tau;\delta,F_{W})$ giving the smallest solution $b\geq 0$ to

3 Predicting operating characteristics from State Evolution

State Evolution offers a formalBy formal, we mean a rule-based procedure which we can follow to get a prediction, without any guarantees that the prediction is correct. procedure for predicting operating characteristics of the AMP iteration at any fixed iteration $t$ or in the limit $t\to\infty$. Nater in this section, we will provide rigorous validation of these predictions.

Call the tuple $S=(\tau;b,\delta,F)$ a state; in running the AMP algorithm we assume that the algorithm is initialized with $\widehat{\theta}^{0}$ so that $\tau_{0}^{2}={\rm MSE}(\widehat{\theta}^{0},\theta_{0})/\delta$, so that AMP starts in state $S=(\tau_{0};b^{0},\delta,F)$, and visits $S_{1}=(\tau_{1};b^{1},\delta,F)$, $S_{2}=(\tau_{2};b^{2},\delta,F)$, …; eventually AMP visits states arbitrarily close to the equilibrium state $S_{*}=(\tau_{*};b^{*},\delta,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.

The state evolution formalism assigns predictions ${\cal E}$ to two types of observables under specific states.

where expectation on the right hand side is with respect to $Z\sim{\sf N}(0,1)$.

where $Z\sim{\sf N}(0,1)$ and $W\sim F_{W}$ is independent of $Z$.

The two most important predictions of operating characteristics are undoubtedly:

${\rm MSE}$ at iteration $t$. We let $S_{t}=(\tau_{t},b(\tau_{t}),\delta,F_{W})$ denote the state of AMP at iteration $t$, and predict

${\rm MSE}$ at convergence. With $\tau_{*}>0$ the limit of $\tau_{t}$, let $S_{*}=(\tau_{*},b(\tau_{*}),\delta,F_{W})$ denote the state of AMP at convergence. and predict

Other predictions might also be of interest. Thus, concerning the mean absolute error ${\rm MAE}(\widehat{\theta}^{t},\theta_{0})=\|\widehat{\theta}^{t}-\theta_{0}\|_{1}/p$, state evolution predicts ${\rm MAE}\approx\sqrt{2\delta\tau_{t}^{2}/\pi}$. Concerning functions of $(R,W)$, consider the ordinary residuals $Y-X\widehat{\theta}^{*}$ at AMP convergence. These residuals will of course in general not have the distribution of the errors $W$. Setting $\eta(z;b)=z-\Psi(z;b)$, we have $Y-X\widehat{\theta}^{*}=\eta(R;b_{*})$. State evolution predicts that the ordinary residuals will have the same distribution as $\eta(W+\tau_{*}Z;b_{*})$.

4 Example of State Evolution predictions

Continuing with our running example, we again consider the case of contaminated normal data $W\sim{\sf CN}(0.05,10)$ and Huber $\rho$ with $\lambda=3$. If we start AMP with the all-zero estimate $\widehat{\theta}^{0}=0$, then since $\|\theta_{0}\|_{2}=6\sqrt{p}$ we start SE with $\tau_{0}=2.056$. Figure 4 presents predictions by state evolution for the MSE (left panel) and for the mean absolute error MAE.

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.

5 A lower bound on State Evolution

State Evolution cannot evolve so that $\tau_{t}^{2}\to 0$; under minimal regularity, it always exceeds a specific nonzero noise level.

Suppose that $F_{W}$ has a well-defined Fisher information $I(F_{W})$. Then for any $t>0$

We can sharpen this bound one step further. It will be convenient to write $I(X)$ for the Fisher information of distribution $F_{X}$.

Barron and Madiman [MB07] give the inequality $I(W+\tau Z)\leq x^{2}I(W)+(1-x)^{2}I(\tau Z)$, valid for any $x\in(0,1)$. By calculus, we know that for $a,b>0$,

Setting $a=I(W)$ and $b=I(\tau Z)=\tau^{-2}$, and dividing both numerator and denominator by $b$, we are done. ∎

Revisit the argument of Lemma 3.4; the inequality $\tau_{t}^{2}\geq 1/(\delta I(F_{W}))$ shows that if $t>0$, then $\tau_{t}^{2}I(W)\geq 1/\delta$. Using this in the previous Lemma,

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$,

Suppose that $F_{W}$ has a well-defined Fisher information $I(F_{W})$. Then for every accumulation point $\tau_{*}$ of State Evolution

6 Correctness of State Evolution predictions

The predictions of state evolution can be validated in the large-system limit $n,p\to\infty$, 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, $\xi(x)=x^{2}$ is pseudo-Lipschitz.

The following result validates the predictions of State Evolution for pseudo-Lipschitz observables. Our proof is deferred to Appendix B.

Assume that the loss function $\rho$ is convex and smooth, that the sequence of matrices $\{{\mathbf{X}}(n)\}_{n}$ is a standard Gaussian design, and that $\theta_{0}$, $\widehat{\theta}^{0}$ are deterministic sequences such that ${\rm AMSE}(\theta_{0}$, $\widehat{\theta}^{0})=\delta\tau_{0}^{2}$. Further assume that $F_{W}$ has finite second moment and let $\{\tau_{t}^{2}\}_{t\geq 0}$ be the state evolution sequence with initial condition $\tau_{0}^{2}$. Let $\{\widehat{\theta}^{t},R^{t}\}_{t\geq 0}$ be the AMP trajectory with parameters $b_{t}$ as per Eq. (19).

In particular, we may take $\xi(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 $\widehat{\theta}$, cf. Eq. (2), is to prove that the AMP iterates $\widehat{\theta}^{t}$ converge to $\widehat{\theta}$. We will prove that this is indeed the case, at least in the limit $n,p\to\infty$, and for suitable initial conditionsWe expect convergence for arbitrary initial conditions (as long as they are independent of $(W,{\mathbf{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)..

The key step is to establish the following high-dimensional convergence result.

(Convergence of AMP to the M-Estimator.) Assume the same setting as in Theorem 3.9, and further assume that $\rho$ is strongly convex and that $\delta>1$.

Let $(\tau_{*},b_{*})$ be a solution of the two equations

and assume that ${\rm AMSE}(\widehat{\theta}^{0},\theta_{0})=\delta\tau_{*}^{2}$. Then

From this and Theorem 3.9, the desired characterization of $\widehat{\theta}$ immediately follows.

To tie back to the introduction, we prove formula (4):

(Asymptotic Variance Formula under High-Dimensional Asymptotics.) Assume the setting of Theorem 3.9, and further assume that $\rho$ is strongly convex and $\delta>1$. The asymptotic variance of $\widehat{\theta}$ obeys

Here $(\tau_{*},b_{*})$ are the unique solutions of the equations (24)-(25).

Information Bound under High-Dimensional Asymptotics:

In this inequality, the effect of the high-dimensional asymptotics parameter $\delta$ is extremely clear; it shows that the classical information bound is not achievable when $\delta=n/p<\infty$, There is always an inflation in variance at least by $(1-1/\delta)^{-1}=\frac{n}{n-p}$. Moreover, the inflation completely blows up as $\delta\to 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:

Assume the same setting as in Theorem 3.9, and further assume that $\rho$ is strongly convex and $\delta>1$. Then the almost sure limits ${\rm AMSE}(\widehat{\theta}^{t};\theta_{0})$ and ${\rm AMSE}(\widehat{\theta};\theta_{0})$ exist, and obey

Theorem 4.1 extends to cover general Gaussian matrices ${\mathbf{X}}$ with i.i.d. rows.

Notice that, if ${\mathbf{X}}$ is a general Gaussian design, then ${\mathbf{X}}\Sigma^{-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].

Assume the same setting as in Theorem 3.9, but with $\{{\mathbf{X}}(n)\}_{n\geq 0}$ being a general Gaussian design with covariance $\Sigma$, and further assume that $\rho$ is strongly convex and $\delta>1$. There is a scalar random variable $T_{n}$ so that

where ${\mathbf{Z}}\sim{\sf N}(0,{\rm I}_{p\times p})$ and we have the almost-sure limit $\lim_{n\to\infty}\ T_{n}=_{a.s.}\tau_{*}$, where $\tau_{*}$ solves Eqs. (24), (25).

This result coincides with Corollary 1 in [EKBBL13] apart from a factor $\sqrt{n}$ in the random part of Eq. (30) that arises because of a difference in the normalization of ${\mathbf{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:

Generalize the i.i.d. Gaussian rows model for ${\mathbf{X}}$ by allowing different rows to be randomly scaled copies of a common $X\sim{\sf N}(0,\Sigma/n)$. This is the setting of [EKBBL13, Result 1].

Remove the smoothness and strong convexity assumptions on $\rho$.

Generalize the present results to non-Gaussian designs. We expect –for instance– that they should hold universally across matrices ${\mathbf{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 ${\cal L}(\theta)$, 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 $\rho$ 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.

with solution $\widehat{\beta}(\widetilde{Y};\widetilde{\mathbf{X}})$, say.

Here is the specific pair that links $(M_{J})$ with $(L_{J})$. We let $\widetilde{\mathbf{X}}$ be a matrix with orthonormal rows such that $\widetilde{\mathbf{X}}{\mathbf{X}}=0$, i.e.

finally, we set $\widetilde{Y}=\widetilde{\mathbf{X}}Y$.

Of special interest is the case $J(x)=\lambda\,|x|$ in which case $(L_{J})$ of (33) defines the Lasso estimator. Then $\rho_{J}(x)=\rho_{{\rm H}}(x;\lambda)$ 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

In this special case, our general result from the next section implies the following:

With problem instances $(Y,X)$ and $(\widetilde{Y},\widetilde{\mathbf{X}})$ related as above, the optimal values of the Lasso problem $(\mbox{Lasso}_{\lambda})$ and the Huber problem $(\mbox{Huber}_{\lambda})$ are identical. The solutions of the two problems are in one-one-relation. In particular, we have

In a sense the Lasso problem solution $\widehat{\beta}$ 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_{\rm adj}\equiv(Y-\hat{\beta})$ with outliers removed.

1.2 General duality result

We will now show that the problem (32) is dual to (33) under or special choice of $(\widetilde{Y},\widetilde{\mathbf{X}})$, via (34).

Assume that $\rho(\,\cdot\,)=\rho_{J}(\,\cdot\,)$, that $\widetilde{\mathbf{X}}$ has orthonormal rows with ${\rm null}(\widetilde{\mathbf{X}})={\rm image}({\mathbf{X}})$, and finally that $\widetilde{Y}=\widetilde{\mathbf{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

‘Differentiating’ Eq. (31) it is easy to see that

We then claim that $\widehat{\beta}\equiv Y-{\mathbf{X}}\widehat{\theta}-u$ is a minimizer of Eq. (33). Indeed

where the last identity follows since, by Eq. (34), ${\rm null}({\mathbf{X}}^{{\sf T}})={\rm image}(\widetilde{\mathbf{X}}^{{\sf T}})$, and hence $u\in{\rm image}(\widetilde{\mathbf{X}}^{{\sf T}})$ by Eq. (41). Using again Eqs. (41) and (40), we deduce that $u\in\partial J(\widehat{\beta})$, i.e.

which is the stationarity condition for the problem (33).

Viceversa a similar argument shows that, given $\widehat{\beta}$ that minimizes Eq. (33), and $\widehat{\theta}\equiv({\mathbf{X}}^{{\sf T}}{\mathbf{X}})^{-1}{\mathbf{X}}^{{\sf T}}(Y-\widehat{\beta})$ is a minimizer of the robust regression problem (32). ∎

2 Comparison to AMP in the p>n𝑝𝑛p>n case

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] $\rho(x)=x^{2}/2$ for $|x|\leq\lambda$ and $\rho(x)=\lambda|x|-\lambda^{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 $\eta(x;\alpha)={\rm sign}(x)\,(|x|-\alpha)_{+}$. We have the relationship

Letting $c_{t}\equiv 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.

Acknowledgements

This work was partially supported by the NSF CAREER award CCF-0743978, the NSF grant DMS-0806211, and the grants AFOSR/DARPA FA9550-12-1-0411 and FA9550-13-1-0036.

Appendix A Properties of the functions 𝖯𝗋𝗈𝗑𝖯𝗋𝗈𝗑{\sf Prox}, ΨΨ\Psi

Since, for $b>0$, $x\mapsto\rho(x)+(x-z)^{2}/(2b)$ is differentiable and strongly convex, $x={\sf Prox}(z;b)$ is uniquely determined by setting to zero the first derivative:

The claim then follows from the Implicit Function theorem. ∎

and hence $\Psi$ is differentiable, with partial derivatives

In particular, for any fixed $b$, $z\mapsto\Psi(z;b)$ is strictly increasing and Lipschitz continuous, with

Using again the stationarity condition (53) that holds for $x={\sf Prox}(z;b)$, we have

which is our first claim. The other claims immediately follow by calculus. ∎

Finally, we prove that Eq. (19) that defines $b_{t}$ as a function of $\tau_{t}$ always has at least one solution.

Then for any $a\in(0,1)$, the set of solutions

It follows immediately from the continuity properties of $\Psi$ that $b\mapsto G(b)$ is continuous. The claim follows by proving that $\lim_{b\to 0}G(b)=0$ and $\lim_{b\to\infty}G(b)=1$.

By Proposition A.2 equation (56) $0\leq\Psi^{\prime}(z;b)\leq 1$. The limit $b\to 0$ follows from dominated convergence since, by the upper bound in (56) $\lim_{b\to 0}\Psi^{\prime}(z;b)=0$ for each $z$.

In order to obtain the limit as $b\to\infty$, note that by Stein Lemma:

Appendix 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 $\widehat{\theta}^{{\rm cen},0}=\widehat{\theta}^{0}-\theta_{0}$, and change data by letting $Y^{{\rm cen}}=Y-{\mathbf{X}}\theta_{0}\equiv W$. Applying the AMP recursion in these new coordinates gives the new trajectory $\widehat{\theta}^{{\rm cen},t}=\widehat{\theta}^{t}-\theta_{0}$ for all $t$, and $R^{{\rm cen},t}=R^{t}$ for all $t$.

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 $\{{\vartheta}^{t},S^{t}\}$ that approximates closely the iterates $\{\widehat{\theta}^{{\rm cen},t},R^{{\rm cen},t}\}$ defined by (64),(65). The new sequence is defined by letting ${\vartheta}^{0}=\widehat{\theta}^{0}-\theta_{0}$ and, for all $t\geq 0$

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.

Under the assumptions of Theorem 3.9, we have, for any fixed $t\geq 0$,

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

where $b^{t}$, $h^{t}$, $q^{t}$ and $m^{t}$ are vectors and $\lambda_{t}$ and $\xi_{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\times N$ random matrix with iid Gaussian entries. The scalars $\xi_{t}=\langle g^{\prime}_{t}(b^{t},w)\rangle$ and $\lambda_{t}=\frac{1}{\delta}\langle f^{\prime}_{t}(h^{t})\rangle$, where $\langle\cdot\rangle$ 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

and the state evolution recursion involves the pair of variables

Table 1 sets up a ‘dictionary’ of correspondences between this paper and [BM11].

We get exact correspondence between the two systems, provided we identify $\delta\Psi(W+S^{t};b_{t})$ with $m^{t}=g_{t}(b^{t};w)$ and $-\delta h^{t}$ with $f_{t}(h^{t})$. One has, in particular, that $\lambda_{t}=\frac{1}{\delta}\langle f^{\prime}_{t}(h^{t})\rangle=-1$, and that $\xi_{t}=\langle g^{\prime}_{t}(b^{t},w)\rangle=\langle\delta\Psi^{\prime}(W+S^{t};b_{t})\rangle=q_{t}$.

Theorem 3.9 now follows from the equivalence of the last two recursions – i.e. equivalence of (64)-(65) with (66)-(67).

Under the assumptions of Theorem 3.9, we have, for any fixed $t\geq 0$,

Comparing the first of these equations with Eq. (64), and using triangular inequality, we get

Comparing analogously Eq. (65) and (75), we obtain

Iterating the upper bounds (77), (79), and using the fact that ${\vartheta}^{0}=\widehat{\theta}^{0}$, we conclude that there exists a constant $A=A(\delta)<\infty$ such that

Finally, it is a standard result in random matrix theory [AGZ09] that $\lim_{n\to\infty}\|{\mathbf{X}}\|_{2}=C(\delta)<\infty$. 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).

Appendix C Proof that AMP converges to the M-estimator (Theorem 4.1)

Notice first of all that, by construction, $\tau_{t}^{2}=\tau_{*}^{2}$, $b_{t}=b_{*}$ for all $t$.

Given $\delta$, $\rho$ as in the statement of the theorem and $\tau_{*}$, $b_{*}$ a solution of the fixed point equation (24), (25), we define the doubly infinite matrix ${\Gamma}=({\Gamma}_{t,s})_{t,s\geq 0}$ by letting, recursively for $t,s\geq 0$

Notice that, in particular, ${\Gamma}_{s,t}={\Gamma}_{t,s}$ for all $s,t\geq 0$ and ${\Gamma}_{t,t}=\tau_{*}^{2}$ for all $t$.

The significance of these quantities is clarified by the following result.

Under the hypotheses of Theorem 3.9, further assume that $\tau_{*}^{2}$ and ${\Gamma}$ are defined as above. Then, for any $t,s\geq 0$,

As a special case of the latter result, we have

The following lemma provides information about the asymptotic behavior of ${\Gamma}_{t,s}$. Its proof is deferred to Section C.2.

Let $\tau_{*}$, ${\Gamma}$ be defined as above for $\delta>1$. Then

(The case $h>1$ follows from $h=1$ by the triangle inequality.)

We are now ready to prove Theorem 4.1. Recall that ${\cal L}(\theta)={\cal L}(\theta;Y,{\mathbf{X}})$ denotes the loss function defined in Eq. (2), and that its gradient and Hessian are given by

In particular, letting $\sigma_{\rm min}({\mathbf{X}})$ denote the minimum non-zero singular value of ${\mathbf{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 $\delta>1$ such that for any $n\geq n_{0}$,

As a consequence, with probability at least $1-e^{-c_{1}n}$, we have

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 $\Psi(z;b_{*})=z-{\sf Prox}(z;b_{*})$, this can be rewritten as

By Eq. (11), and recalling that $\Psi(z;b)=b\,\rho^{\prime}({\sf Prox}(z;b))$, we have

Hence, using Eqs (91) and (91), and recalling that $\lim_{n\to\infty}\|{\mathbf{X}}\|_{2}<\infty$ almost surely [AGZ09], we get

This is equivalent to the claim (99) since $\nabla_{\theta}{\cal L}(\theta)=-{\mathbf{X}}\rho^{\prime}(Y-{\mathbf{X}}\theta)$.

First of all note that, due to Lemma B.2, it is sufficient to prove that

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].

Here ${\mathbf{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}\equiv{\Gamma}_{t,t+1}/\tau_{*}^{2}$. We then have the recursion

${\sf H}(q)$ is strictly convex for $q\in$.

In order to prove $(i)$, note that, for $q=1$, $Z_{1}=Z_{2}\equiv Z\sim{\sf N}(0,1)$ and hence

which is equal to $1$ since $b_{*}$, $\tau_{*}$ satisfy Eq. (24).

In order to prove $(ii)$, $(iii)$, define

Then we have the spectral representation (for $t=\log(1/q)$)

Because of the remarks $(i)$-$(iii)$ just proven, it follows that $\lim_{t\to\infty}q_{t}=1$ (and hence $\lim_{t\to\infty}{\Gamma}_{t,t+1}=\tau_{*}^{2}$) if and only if ${\sf H}^{\prime}(1)\leq 1$. A simple calculation yields

where $Z\sim{\sf N}(0,1)$. Recalling that $\Psi^{\prime}(z;b)\in(0,1)$, we have $(\Psi^{\prime})^{2}\leq\Psi^{\prime}$ and so

where the last identity follows because $(\tau_{*},b_{*})$ solve Eq. (25). This finishes the proof.

References