Radial M-Estimation Of Multivariate Location

Why Weight Residuals by Their Distance?

The sample mean responds to a residual in direct proportion to its length: a point ten times farther from a proposed center contributes ten times as much to the mean equation. A radial M-estimator changes that relationship by assigning a weight according to the residual’s distance from the center. “Radial” means that the weight depends on distance, not direction. This preserves symmetry across directions while limiting how strongly distant observations can pull the estimate.

We study location with a fixed, known or previously chosen scatter matrix. Jointly estimating location and scatter is a different problem. Let $x_1,\ldots,x_n$ be observations in $\mathbb{R}^p$, and let $S$ be a fixed symmetric positive-definite $p\times p$ matrix. The distance from observation $x_i$ to a proposed location $\mu$ is its Mahalanobis distance using $S$.

Definition: Let $\rho:[0,\infty)\to\mathbb{R}$ be a differentiable radial loss, with $\rho(0)$ finite. A radial location M-estimator is a minimizer, or a stationary solution, of $$ Q(\mu)=\sum_{i=1}^n \rho\bigl(d_i(\mu;S)\bigr), \qquad d_i(\mu;S)=\sqrt{(x_i-\mu)^\mathsf{T}S^{-1}(x_i-\mu)}. $$ Writing $\psi(d)=\rho'(d)$, define $u(d)=\psi(d)/d$ for $d>0$. When $u$ has a finite limit at zero, use that limit as $u(0)$.

The function $\rho$ penalizes a residual distance, while $u$ is the weight that appears in the location equation. When $\rho(d)=d^2/2$, we have $u(d)=1$, so the equation reduces to the ordinary mean equation. Robust choices make $u(d)$ decrease for sufficiently large distances.

Deriving the Estimating Equation

Set $r_i=x_i-\mu$. At any residual with $d_i>0$, differentiation gives $$ \nabla_\mu d_i=-\frac{S^{-1}r_i}{d_i}. $$ The chain rule therefore yields the gradient of the objective. If a residual is zero, the same result follows by continuity when $u(0)$ is finite.

\nabla_\mu Q(\mu) =-\sum_{i=1}^n \frac{\psi(d_i)}{d_i}S^{-1}r_i =-S^{-1}\sum_{i=1}^n u(d_i)r_i.

Since $S^{-1}$ is nonsingular, a stationary location must solve the following estimating equation. Conversely, every solution of this equation is stationary for $Q$ under the stated differentiability conditions. A stationary solution need not be a global minimizer: a nonconvex radial loss can have multiple stationary solutions.

Theorem: For fixed positive-definite $S$, any differentiable stationary solution of $Q$ satisfies $$ \sum_{i=1}^n u\bigl(d_i(\mu;S)\bigr)(x_i-\mu)=0. $$ If the weights at a proposed solution are nonnegative and have positive sum, this equation is equivalent, at that solution, to $$ \mu= \frac{\sum_{i=1}^n u\bigl(d_i(\mu;S)\bigr)x_i} {\sum_{i=1}^n u\bigl(d_i(\mu;S)\bigr)}. $$

To verify the weighted-mean form, expand $x_i-\mu$ in the estimating equation: $$ \sum_i u(d_i)x_i-\mu\sum_i u(d_i)=0. $$ Division by the positive weight sum gives the displayed expression. The weights depend on $\mu$, so this is an implicit equation, not an ordinary formula with weights known in advance. Iteratively reweighted methods update the distances, recalculate weights, and repeat; an individual update is not automatically the final M-estimator.

Worked Example: One Reweighted Update in Two Dimensions

Take $S=I$, the identity matrix, and observations $(0,0),(2,0),(0,2),(8,8)$. Use the Huber radial loss with threshold $c=2$, whose score and weight are $$ \psi(d)=\min(d,2),\qquad u(d)= \begin{cases} 1,&d\leq 2,\\ 2/d,&d>2. \end{cases} $$ Start an iteration at $\mu^{(0)}=(0,0)$.

The four distances from this starting location are $0,2,2,\sqrt{128}=8\sqrt{2}$. Thus the weights are $1,1,1,2/(8\sqrt{2})=1/(4\sqrt{2})$. The weighted-mean update is $$ \mu^{(1)} = \frac{(2,0)+(0,2)+(1/(4\sqrt{2}))(8,8)} {3+1/(4\sqrt{2})} = \left( \frac{2+\sqrt{2}}{3+1/(4\sqrt{2})}, \frac{2+\sqrt{2}}{3+1/(4\sqrt{2})} \right) \approx(1.075,1.075). $$

The update is a solution to the weighted subproblem using weights computed at $\mu^{(0)}$. It is not yet a solution to the full implicit equation. For instance, its distance to $(8,8)$ is approximately $\sqrt{2}(8-1.075)=9.79$, so that point’s next weight is $2/9.79\approx0.204$, not the starting weight $1/(4\sqrt{2})\approx0.177$. A further iteration must recompute all distances and weights.

Residual Scores and Observation Influence

A point’s weight alone does not describe its contribution to the estimating equation. In coordinates whitened by $S$, define $z=S^{-1/2}(x-\mu)$, so that $d=\lVert z\rVert$. The radial score vector is $$ \Psi(z)=u(\lVert z\rVert)z. $$ Its length is $u(d)d=\psi(d)$ when the score is nonnegative. Thus a point can be far away and have a small weight while still exerting a bounded, nonzero score. For the Huber loss, $\lVert\Psi(z)\rVert=\min(d,c)$, so the score length never exceeds $c$.

Result: For the Huber radial loss, the score contribution of any residual is bounded in norm by $c$. Indeed, if $d\leq c$, then $\lVert\Psi(z)\rVert=\lVert z\rVert=d\leq c$. If $d>c$, then $u(d)=c/d$, so $\lVert\Psi(z)\rVert=(c/d)\lVert z\rVert=c$. Hence $\lVert\Psi(z)\rVert\leq c$ in both cases.

This bound is a statement about each residual’s score, not a guarantee that every fitted estimator is stable. The behavior of the whole estimate also depends on the slope of the estimating equation near its solution. The following local calculation makes that distinction explicit.

Result: Suppose a location functional $T(F)$ solves $E_F[\Psi(X-T(F))]=0$ in whitened coordinates, and consider a solution $\theta=T(F)$. If $\Psi$ is differentiable at the relevant residuals, differentiation under the expectation is valid in a neighborhood of $\theta$, and $$ A=E_F[D\Psi(X-\theta)] $$ is nonsingular, then its influence function at $x$ is $$ \operatorname{IF}(x;T,F)=A^{-1}\Psi(x-\theta). $$ Here $D\Psi$ is the matrix of first derivatives.

To derive the formula, contaminate $F$ by a point mass at $x$, writing $F_\varepsilon=(1-\varepsilon)F+\varepsilon\Delta_x$ and $\theta_\varepsilon=T(F_\varepsilon)$. The defining equation is $$ (1-\varepsilon)E_F[\Psi(X-\theta_\varepsilon)] +\varepsilon\Psi(x-\theta_\varepsilon)=0. $$ Differentiate at $\varepsilon=0$. The first expectation is zero at $\theta$; its derivative through $\theta_\varepsilon$ is $-A\theta'_0$. The derivative of the factor $1-\varepsilon$ multiplies that zero expectation, and the derivative of the second term is $\Psi(x-\theta)$. Consequently, $ -A\theta'_0+\Psi(x-\theta)=0$. Since $A$ is nonsingular, $\theta'_0=A^{-1}\Psi(x-\theta)$, as claimed. This derivation requires the stated differentiation and nonsingularity conditions; a bounded score alone does not supply them. $\square$

Worked Example: A Distant Point with Bounded Score

Take $S=I$, Huber threshold $c=1$, and the six observations $$ (0.5,0),\;(-0.5,0),\;(0,0.5),\;(0,-0.5),\;(10,0),\;(-10,0). $$ At $\mu=(0,0)$, the four inner residuals have distance $0.5$ and weight $1$. The two outer residuals have distance $10$ and weight $1/10$. By symmetry, the weighted residuals sum to zero, so the origin solves the estimating equation. The sum of weights is $4+0.2=4.2$, and the weighted-mean expression also gives $(0,0)$.

An inner point contributes a score of length $1(0.5)=0.5$. An outer point contributes a score of length $(1/10)(10)=1$, rather than the length $10$ that an unweighted residual would contribute. The outer observation is still influential, but its score has been capped at $c=1$.

We can also calculate the local sensitivity of the sample equation. For $d<c$, $\Psi(z)=z$, so $D\Psi=I$; the four inner points contribute $4I$. For $d>c$, $\Psi(z)=cz/d$, with derivative $$ D\Psi(z)=\frac{c}{d}\left(I-\frac{zz^\mathsf{T}}{d^2}\right). $$ Each outer point lies on the first coordinate axis, so its derivative is $\operatorname{diag}(0,0.1)$. The total derivative matrix is therefore $\operatorname{diag}(4,4.2)$. A score perturbation of one outer case in the positive first-coordinate direction is $(1,0)$; the corresponding first-order location shift is $$ \operatorname{diag}(4,4.2)^{-1}(1,0)^\mathsf{T}=(0.25,0). $$ This is a local linearized effect, not an exact replacement estimate. It shows both how the bounded score limits the pull of the distant point and how the equation’s derivative calibrates that pull.

When Radial Weights Can Mislead

The weight is not the same thing as influence. A large-distance residual can have a small weight but a score contribution $u(d)d$ that is still substantial. For Huber weighting that contribution reaches a constant bound; with a redescending score it may eventually decrease. The local influence calculation also involves $A^{-1}$, so a nearly singular derivative matrix can magnify score contributions.

Redescending losses require particular care. If $u(d)$ becomes zero for sufficiently large $d$, an observation can cease to affect the equation once its residual is large enough. This can be useful for limiting the effect of extreme points, but it can also create multiple stationary solutions. An algorithm’s starting value and convergence checks then matter; a root of the equation should not be assumed to be the unique global minimum.

Common pitfall: “The residual distance is ten times larger, so its influence must be ten times larger.” That reasoning is correct for an unweighted mean equation, but not for a radial M-estimator. The relevant score length is $u(d)d$, and local influence also depends on the inverse derivative matrix. Compute the weight, score, and local equation slope separately.

Check Your Understanding

Use these questions to check whether you can move from the loss to the estimating equation and interpret its weights.

For a fixed positive-definite scatter matrix $S$, derive the radial location equation from $Q(\mu)=\sum_i\rho(d_i(\mu;S))$.
For Huber radial loss with threshold $c$, what are the weight and score length for a residual distance $d>c$?
Why is one weighted-mean update generally not enough to establish that the implicit M-estimating equation has been solved?
In whitened coordinates, state the conditions in the influence-function result that ensure the local formula is defined.
Why can a redescending score make the choice of starting value important?

Radial M-Estimation Of Multivariate Location

What You'll Learn