Asymptotic Linearization of Robust Multivariate Estimates

Why Linearize a Robust Estimate?

A robust multivariate estimator is often defined implicitly: its estimate is the parameter value that makes a vector of sample estimating equations equal to zero. That definition can be useful computationally, but it does not by itself reveal how the estimate varies from sample to sample. Asymptotic linearization supplies that missing link. It approximates the estimator by an average of independent, mean-zero contributions, making its large-sample variability calculable.

The key is to distinguish two features of the estimating equation. Its derivative describes how strongly the equation responds when the parameter changes; its score covariance describes how much the equation varies across observations. These are different matrices, and both enter the resulting covariance.

Estimating Equations and Their Matrices

Definition: Let $X_1,\ldots,X_n$ be independent and identically distributed observations from a distribution $F$. Let $\theta\in\mathbb{R}^p$ be the parameter of interest and let $\psi(x,\theta)\in\mathbb{R}^p$ be an estimating function. Write $$ P_n\psi(\theta)=\frac{1}{n}\sum_{i=1}^n\psi(X_i,\theta), \qquad G(\theta,F)=E_F[\psi(X,\theta)]. $$ A population target $\theta_0$ satisfies $G(\theta_0,F)=0$. An estimator $\hat\theta$ solves $P_n\psi(\hat\theta)=0$, or solves it up to a negligible residual.

For the results below, assume $E_F\|\psi(X,\theta_0)\|^2<\infty$. Since the population estimating equation is zero at $\theta_0$, the variability matrix is

B=\operatorname{Var}_F\{\psi(X,\theta_0)\} =E_F[\psi(X,\theta_0)\psi(X,\theta_0)^T].

Suppose $G$ is differentiable at $\theta_0$. Define the sensitivity matrix, with a sign convention chosen to make it positive in many common location problems, by

A=-\left.\frac{\partial G(\theta,F)}{\partial\theta^T}\right|_{\theta=\theta_0}.

Both $A$ and $B$ are $p\times p$ matrices. In general, $A$ need not be symmetric. The transpose in the final covariance formula is therefore essential. Bounded estimating functions are common in robust procedures, but boundedness alone does not establish the regularity conditions or guarantee that $A$ is nonsingular.

Asymptotic Linearization

Theorem: Suppose $\theta_0$ is an interior point of the parameter space, $G(\theta_0,F)=0$, and $A$ is nonsingular. Assume $\hat\theta\xrightarrow{p}\theta_0$, and that $\psi(X_i,\theta)$ is differentiable in $\theta$ along the line segment from $\theta_0$ to $\hat\theta$, with the empirical Jacobian $P_n\dot\psi(\theta)=n^{-1}\sum_i\partial\psi(X_i,\theta)/\partial\theta^T$ defined there. Require the following local convergence condition: for every $\eta>0$, $$ \lim_{r\downarrow 0}\limsup_{n\to\infty} P\left( \sup_{\|\theta-\theta_0\|\le r} \|P_n\dot\psi(\theta)+A\|>\eta \right)=0. $$ Finally, suppose the equation residual satisfies $\|P_n\psi(\hat\theta)\|=o_p(n^{-1/2})$. Then $$ \hat\theta-\theta_0 =A^{-1}P_n\psi(\theta_0)+o_p(n^{-1/2}) =\frac{1}{n}\sum_{i=1}^n A^{-1}\psi(X_i,\theta_0)+o_p(n^{-1/2}). $$ If $B$ is as defined above, then $$ \sqrt n(\hat\theta-\theta_0)\xrightarrow{d} N_p(0,\ A^{-1}BA^{-T}). $$

The local convergence condition is stated directly for the empirical Jacobian near $\theta_0$. It is stronger than merely saying that a uniform law of large numbers holds for some Jacobian: it ensures that the Jacobian along the random segment to $\hat\theta$ approaches $-A$. In applications, it can follow from a suitable uniform law of large numbers together with continuity of the population Jacobian at $\theta_0$.

Proof. Put $\delta_n=\hat\theta-\theta_0$. By consistency, the segment $\theta_0+t\delta_n$, $0\le t\le1$, lies in an arbitrarily small neighborhood of $\theta_0$ with probability tending to one. The local condition therefore implies

J_n=\int_0^1 P_n\dot\psi(\theta_0+t\delta_n)\,dt \xrightarrow{p}-A.

The fundamental theorem of calculus applied along this segment gives

P_n\psi(\hat\theta) =P_n\psi(\theta_0)+J_n\delta_n.

Let $e_n=P_n\psi(\hat\theta)$. The multivariate central limit theorem, using $E_F[\psi(X,\theta_0)]=0$ and its finite second moment, implies $P_n\psi(\theta_0)=O_p(n^{-1/2})$. Also $e_n=o_p(n^{-1/2})$. Since $J_n\xrightarrow{p}-A$ and $A$ is nonsingular, $J_n$ is invertible with probability tending to one and $J_n^{-1}=O_p(1)$. Rearranging the displayed equation yields $\delta_n=J_n^{-1}\{e_n-P_n\psi(\theta_0)\}=O_p(n^{-1/2})$.

Now write $J_n=-A+R_n$, where $R_n=o_p(1)$. Substitution gives $e_n=P_n\psi(\theta_0)-A\delta_n+R_n\delta_n$. Because $R_n\delta_n=o_p(n^{-1/2})$, rearranging proves the asserted linear representation. Finally, the summands $A^{-1}\psi(X_i,\theta_0)$ are independent, identically distributed, have mean zero and covariance $A^{-1}BA^{-T}$. The multivariate central limit theorem applied to their normalized sum gives the stated limit. $\square$

The covariance $A^{-1}BA^{-T}$ is called the sandwich covariance: the sensitivity matrix and its transpose surround the variability matrix. It is the asymptotic covariance of $\sqrt n(\hat\theta-\theta_0)$; the approximate covariance of $\hat\theta$ itself is $A^{-1}BA^{-T}/n$.

Worked Covariance Calculations

Worked Example: A Bounded-Score Location Estimate

Let $X=(X_1,X_2)^T$, with independent coordinates. Each coordinate takes values $1$ and $-1$, each with probability $0.4$, and $5$ and $-5$, each with probability $0.1$. Estimate a two-dimensional location $\theta$ using the componentwise score $\psi_j(X,\theta)=\phi_2(X_j-\theta_j)$, where $\phi_2(u)=\max(-2,\min(u,2))$.

Symmetry gives $G(0,F)=0$, so $\theta_0=0$ is a target. Near zero, the observations at $\pm1$ are in the linear part of the score, while those at $\pm5$ remain clipped. Thus the derivative of each score with respect to its location coordinate is $-1$ for the central observations and zero for the outer observations. The central probability is $0.8$, so

$$ A=0.8I_2. $$

For either coordinate the score takes values $1,-1,2,-2$ with probabilities $0.4,0.4,0.1,0.1$. Its mean is zero and its second moment is $0.4(1^2)+0.4((-1)^2)+0.1(2^2)+0.1((-2)^2)=1.6$. Independence makes the cross-covariance zero. Hence $B=1.6I_2$, and

A^{-1}BA^{-T} =(1/0.8)(1.6)(1/0.8)I_2 =2.5I_2.

The estimator therefore has approximate covariance $2.5I_2/n$. The calculation uses score variability $B$, not the covariance of the raw observations: the score clips the outer observations before they contribute to the estimating equation.

Worked Example: Why the Sandwich Uses a Transpose

Use the same distribution and componentwise score as above, but premultiply the score by $C=\begin{pmatrix}1&1\\0&1\end{pmatrix}$, defining $\psi_C=C\psi$. The root is unchanged because $C$ is nonsingular, but the sensitivity is not symmetric.

Here $A_C=CA=0.8C=\begin{pmatrix}0.8&0.8\\0&0.8\end{pmatrix}$, and $B_C=CBC^T=1.6CC^T=\begin{pmatrix}3.2&1.6\\1.6&1.6\end{pmatrix}$. Since $A_C^{-1}=1.25C^{-1}$ and $A_C^{-T}=1.25C^{-T}$, direct substitution gives

A_C^{-1}B_CA_C^{-T} =(1.25C^{-1})(1.6CC^T)(1.25C^{-T}) =2.5I_2.

The result agrees with the first calculation, as it must: premultiplying the estimating equation by a fixed nonsingular matrix does not change its solution. Omitting the transpose, or using $A_C^{-1}B_CA_C^{-1}$, would not in general produce the covariance of that solution.

Influence Function Under Point Contamination

The same sensitivity matrix describes a local response to contamination. For a point $z$, define $F_\varepsilon=(1-\varepsilon)F+\varepsilon\Delta_z$, where $\Delta_z$ places all its probability at $z$. Suppose $G(\theta,F)$ is differentiable at $\theta_0$, $A$ is nonsingular, and a solution path $\theta_\varepsilon$ of $G(\theta_\varepsilon,F_\varepsilon)=0$ is right-differentiable at zero with $\theta_{\varepsilon=0}=\theta_0$. In addition, require $\psi(z,\theta)$ to be continuous at $\theta_0$. This last condition ensures that the contribution from the contaminating point has the required limit as the parameter moves.

Result: Under these conditions, the influence function, defined as the right derivative of the target along this contamination path, is $$ \operatorname{IF}(z;T,F) =\left.\frac{d\theta_\varepsilon}{d\varepsilon}\right|_{\varepsilon=0+} =A^{-1}\psi(z,\theta_0). $$

To verify the formula, write the contaminated equation as $0=(1-\varepsilon)G(\theta_\varepsilon,F)+\varepsilon\psi(z,\theta_\varepsilon)$. Let $q=\left.d\theta_\varepsilon/d\varepsilon\right|_{0+}$. Differentiability of $G$ gives $G(\theta_\varepsilon,F)=-A(\theta_\varepsilon-\theta_0)+o(\|\theta_\varepsilon-\theta_0\|)$, and the assumed path differentiability gives $(\theta_\varepsilon-\theta_0)/\varepsilon\to q$. Continuity at the contaminating point gives $\psi(z,\theta_\varepsilon)\to\psi(z,\theta_0)$. Divide the equation by $\varepsilon$ and take the right limit: $0=-Aq+\psi(z,\theta_0)$. Nonsingularity of $A$ yields $q=A^{-1}\psi(z,\theta_0)$. $\square$

Worked Example: Influence of a Specific Contaminating Point

Return to the bounded-score model above and take $z=(5,1)^T$. At $\theta_0=0$, the score is $\psi(z,0)=(2,1)^T$, while $A=0.8I_2$. Therefore

\operatorname{IF}(z;T,F) =(0.8I_2)^{-1} \begin{pmatrix}2\\1\end{pmatrix} = \begin{pmatrix}2.5\\1.25\end{pmatrix}.

For a direct check, the first coordinate of the population equation has derivative $-0.8$ at zero, and the contaminating score in that coordinate is $2$. Differentiating the contaminated equation gives $-0.8q_1+2=0$, so $q_1=2.5$. The second coordinate gives $-0.8q_2+1=0$, so $q_2=1.25$. The clipped score is continuous at this point, as required for the contamination calculation.

Using the Expansion and Avoiding Common Errors

The linear representation identifies the observation-level contribution $A^{-1}\psi(X_i,\theta_0)$. It explains why a robust procedure can limit the direct effect of extreme observations through a bounded score, while still having sampling uncertainty determined by both $A$ and $B$. In practice, estimates of these matrices can be substituted into the sandwich formula when their consistency is justified for the model and estimating procedure at hand.

Using only the score covariance. $B$ describes variability of the estimating equation, not of the estimator. The estimator covariance also requires the sensitivity adjustment on both sides.
Assuming differentiability at one point proves the expansion. A Taylor argument involves the Jacobian along a segment ending at the random estimate. A local convergence condition, such as the one stated above, is needed to control that matrix.
Differentiating contamination without continuity. The contaminating term is $\psi(z,\theta_\varepsilon)$. Its limit is $\psi(z,\theta_0)$ only when suitable continuity or another explicit justification is available.
Ignoring the transpose. The sandwich is $A^{-1}BA^{-T}$, including when $A$ is nonsymmetric.

Key takeaway: Under suitable local regularity, a robust estimating-equation solution is asymptotically an average of transformed scores. The sensitivity $A$ transforms the score, while the variability matrix $B$ records its sampling variation; together they determine the sandwich covariance and the point-contamination influence function.

Check Your Understanding

Use the definitions and conditions in this tutorial to answer the following.

For an estimating equation with mean-zero score, what does $B$ measure, and what does $A$ measure?
Why is a local condition on the empirical Jacobian needed in the proof of the asymptotic linear representation?
If $A$ is nonsymmetric, write the asymptotic covariance of $\sqrt n(\hat\theta-\theta_0)$ with the correct matrix order.
In the bounded-score example, explain why the outer observations contribute to $B$ even though they contribute zero to the local derivative defining $A$.
What condition on the score at a contaminating point justifies replacing $\psi(z,\theta_\varepsilon)$ by $\psi(z,\theta_0)$ in the influence-function derivation?