Computing Multivariate M-Estimators Reliably

Why Computation Is Part of the Estimator

A multivariate location-scatter M-estimator is defined by estimating equations, but in practice it is usually found by iteration. At each step, observations receive weights determined by their current Mahalanobis distances; those weights then define an updated location and scatter. This feedback makes computation consequential: a poor starting matrix can make distances unusable, and a converged iteration can settle at a root that does not represent the solution an analyst intended to find.

The goal here is not just to obtain a pair of matrices that stop changing numerically. A reliable computation checks that distances and weights are well-defined, the scatter remains positive definite, and the reported solution is stable under reasonable changes in initialization and numerical settings.

Fixed-Point Updates and Their Preconditions

Definition: Let \(x_1,\ldots,x_n\in\mathbb{R}^p\) be observations. For a location vector \(\mu\) and a positive-definite scatter matrix \(\Sigma\), define \(d_i^2(\mu,\Sigma)=(x_i-\mu)^\mathsf{T}\Sigma^{-1}(x_i-\mu)\). Given a nonnegative weight function \(u\) and a fixed multiplier \(b>0\), the location-scatter fixed-point equations are

\mu=\frac{\sum_{i=1}^n u(d_i^2)x_i}{\sum_{i=1}^n u(d_i^2)},\qquad \Sigma=\frac{b}{n}\sum_{i=1}^n u(d_i^2)(x_i-\mu)(x_i-\mu)^\mathsf{T}.

These equations require \(\sum_i u(d_i^2)>0\) for the location update and require the resulting scatter matrix to be positive definite if it is to be used in the next Mahalanobis-distance calculation. The multiplier \(b\) is held fixed during iteration; its calibration is a separate issue.

The standard iteratively reweighted scheme starts from \((\mu^{(0)},\Sigma^{(0)})\), computes all distances and weights from that pair, then updates the location and scatter:

w_i^{(t)}=u\!\left((x_i-\mu^{(t)})^\mathsf{T}(\Sigma^{(t)})^{-1}(x_i-\mu^{(t)})\right),\qquad \mu^{(t+1)}=\frac{\sum_i w_i^{(t)}x_i}{\sum_i w_i^{(t)}},

\Sigma^{(t+1)}=\frac{b}{n}\sum_i w_i^{(t)} (x_i-\mu^{(t+1)})(x_i-\mu^{(t+1)})^\mathsf{T}.

Thus the scatter update uses the old weights but centers residuals at the newly updated location. The next iteration recomputes distances using both updated quantities. An iteration is not defined if the current scatter matrix cannot be inverted, or if all weights vanish. Even when every step is defined, convergence is not automatic: redescending weight functions can make the estimating equations nonconvex and permit multiple roots.

What the Weighted Updates Guarantee

The location update is the weighted least-squares center for the current weights. For fixed nonnegative weights with positive sum, define \(Q(a)=\sum_i w_i\|x_i-a\|^2\). Expanding around the weighted mean \(\bar{x}_w=\sum_iw_ix_i/\sum_iw_i\) gives

Q(a)=\sum_iw_i\|x_i-\bar{x}_w\|^2+ \left(\sum_iw_i\right)\|a-\bar{x}_w\|^2.

The second term is nonnegative and equals zero only at \(a=\bar{x}_w\). This proves that the update uniquely minimizes \(Q(a)\) when the total weight is positive. It does not prove that the full location-scatter iteration decreases one common objective: the weights change with the parameters, so the fixed-weight minimization property alone is not a convergence theorem.

Result: For fixed weights \(w_i\geq 0\), the weighted scatter update \(S=(b/n)\sum_i w_i(x_i-a)(x_i-a)^\mathsf{T}\) is positive definite if and only if the residual vectors \(x_i-a\) with \(w_i>0\) span \(\mathbb{R}^p\).

To prove the result, take any vector \(v\in\mathbb{R}^p\). Then

v^\mathsf{T}Sv=\frac{b}{n}\sum_iw_i\{v^\mathsf{T}(x_i-a)\}^2.

Every term is nonnegative. If the positive-weight residuals span \(\mathbb{R}^p\), then for each nonzero \(v\) at least one positive-weight residual has a nonzero inner product with \(v\); hence \(v^\mathsf{T}Sv>0\), so \(S\) is positive definite. Conversely, if those residuals do not span \(\mathbb{R}^p\), there is a nonzero vector \(v\) orthogonal to all of them. Then \(v^\mathsf{T}Sv=0\), so \(S\) is not positive definite. This criterion makes the danger of too few positive weights precise: the effective observations must provide variation in every direction.

Worked Example: One Reweighted Update

Take four observations in two dimensions: \(x_1=(0,0)^\mathsf{T}\), \(x_2=(1,0)^\mathsf{T}\), \(x_3=(0,1)^\mathsf{T}\), and \(x_4=(3,0)^\mathsf{T}\). Start at \(\mu^{(0)}=(0,0)^\mathsf{T}\), \(\Sigma^{(0)}=I\), and use \(u(t)=\min(1,4/t)\), with \(u(0)=1\). Let \(b=1\).

The initial squared distances are \(0,1,1,9\), so the weights are \(1,1,1,4/9\), respectively. Their sum is \(31/9\). The updated location is therefore

\mu^{(1)}= \left(\frac{1+3(4/9)}{31/9},\frac{1}{31/9}\right)^\mathsf{T} =\left(\frac{21}{31},\frac{9}{31}\right)^\mathsf{T}.

For the scatter, the weighted sum of the uncentered outer products is \(\operatorname{diag}(5,1)\). Centering at the weighted mean subtracts \((31/9)\mu^{(1)}(\mu^{(1)})^\mathsf{T}\), giving

\sum_iw_i(x_i-\mu^{(1)})(x_i-\mu^{(1)})^\mathsf{T} = \begin{pmatrix}106/31&-21/31\\-21/31&22/31\end{pmatrix}.

Dividing by \(n=4\) gives \(\Sigma^{(1)}=\begin{pmatrix}53/62&-21/124\\-21/124&11/62\end{pmatrix}\). Its determinant is \(53(11)/62^2-(21/124)^2=0.1444\) to four decimal places, which is positive; its leading diagonal entry is positive as well, so this two-by-two matrix is positive definite. This is one valid update, not evidence that the iteration has converged. The next weights must be calculated using \(\mu^{(1)}\) and \(\Sigma^{(1)}\).

Singular Starts and Unstable Weights

Mahalanobis distances require an invertible scatter matrix. A singular start is not repaired by attempting the inverse anyway: it makes the distance undefined in at least one direction. Near singularity can be nearly as damaging. In an eigenvector direction with a very small eigenvalue, the inverse scatter magnifies even modest residual components, often turning many observations into extreme outliers according to the weights.

Worked Example: A Nearly Singular Starting Matrix

Consider the starting matrix \(\Sigma_\varepsilon=\begin{pmatrix}1+\varepsilon&1\\1&1+\varepsilon\end{pmatrix}\), with \(\varepsilon>0\), and the observation \(x=(1,0)^\mathsf{T}\) relative to starting location zero. The eigenvalues of this matrix in directions \((1,1)^\mathsf{T}\) and \((1,-1)^\mathsf{T}\) are \(2+\varepsilon\) and \(\varepsilon\). Resolving \(x\) into those two unit directions gives squared distance

d^2=\frac{1}{2}\left(\frac{1}{2+\varepsilon}+\frac{1}{\varepsilon}\right).

At \(\varepsilon=0.01\), this is approximately \(50.249\). With the weight \(u(t)=\min(1,4/t)\), the observation receives weight approximately \(4/50.249=0.0796\). If \(\varepsilon\) decreases further, this weight approaches zero. At \(\varepsilon=0\), the matrix is singular and the distance is undefined. A tiny eigenvalue can therefore change the effective data sharply, even though the entries of the starting matrix appear innocuous.

A practical response is to obtain a full-rank initial scatter estimate, check its eigenvalues and condition number, and reject or repair a numerically singular start before calculating distances. Adding a small ridge term \(\lambda I\) can make an inverse computable, but it changes the starting geometry and is not affine-equivariant in general. Report and assess that choice rather than treating it as a neutral numerical fix.

Weights also need direct inspection. A robust weight function may legitimately assign a small weight to an outlier, but an abrupt pattern of nearly zero weights can leave too few observations to span the space. A useful iteration log includes the smallest eigenvalue of \(\Sigma^{(t)}\), its condition number, the number of positive or practically nonzero weights, and the range of the weights. These checks can reveal a failing update before a matrix inversion or a final estimate becomes misleading.

Convergence Checks, Damping, and Multiple Solutions

A stopping rule should check both location and scatter, because stability in one does not imply stability in the other. One possible relative criterion is

\frac{\|\mu^{(t+1)}-\mu^{(t)}\|}{1+\|\mu^{(t)}\|}<\tau_\mu,\qquad \frac{\|\Sigma^{(t+1)}-\Sigma^{(t)}\|_F}{1+\|\Sigma^{(t)}\|_F}<\tau_\Sigma,

where \(\tau_\mu,\tau_\Sigma>0\) are declared tolerances and \(\|\cdot\|_F\) is the Frobenius norm. These inequalities only say that consecutive iterates are close. Before reporting convergence, also verify that the estimating equations have small residuals, the scatter is positive definite and adequately conditioned, and the weights and distances are finite. If the algorithm stops at a maximum iteration count without these checks passing, report nonconvergence rather than labeling the last iterate a solution.

Damping can reduce abrupt changes. Given an undamped proposal \((\mu^\star,\Sigma^\star)\), one may set \(\mu_{\text{new}}=(1-\alpha)\mu^{(t)}+\alpha\mu^\star\) and \(\Sigma_{\text{new}}=(1-\alpha)\Sigma^{(t)}+\alpha\Sigma^\star\), where \(0<\alpha\leq1\). If both scatter matrices are positive definite, their convex combination is positive definite: for every nonzero \(v\), \(v^\mathsf{T}\Sigma_{\text{new}}v=(1-\alpha)v^\mathsf{T}\Sigma^{(t)}v+\alpha v^\mathsf{T}\Sigma^\star v>0\). Damping therefore preserves positive definiteness under that condition. It does not guarantee convergence to a desired root, and the weights must still be recomputed consistently from the damped parameters.

Worked Example: Two Distinct Fixed Points

This construction shows why a converged root need not be the intended one. Use the redescending weight \(u(t)=(1-t)^2\) for \(0\leq t<1\), and \(u(t)=0\) for \(t\geq1\). Let \(b=27\), \(n=8\), and take the observations \((c+1,0),(c-1,0),(c,1),(c,-1)\) for each of the two centers \(c=10\) and \(c=-10\). Consider first \(\mu=(10,0)^\mathsf{T}\), \(\Sigma=3I\).

The four observations centered at \(10\) each have squared distance \(1/3\), hence weight \((1-1/3)^2=4/9\). The four observations centered at \(-10\) have first-coordinate residuals with absolute values \(19\) or \(21\), so their squared distances exceed \(19^2/3>1\) and their weights are zero. The positive-weight observations have weighted mean \((10,0)^\mathsf{T}\), so the location update leaves \(\mu\) unchanged. Their residual outer products sum to \(2I\). The scatter update is

\frac{27}{8}\frac{4}{9}(2I)=3I.

Thus \(((10,0)^\mathsf{T},3I)\) is exactly a fixed point. By symmetry, \(((-10,0)^\mathsf{T},3I)\) is another. An algorithm initialized at either fixed point will remain there, even though the data contain two separated clusters and neither root alone summarizes both. This is not a numerical roundoff error: the redescending weights exclude the other cluster exactly. In practice, use multiple plausible starts, compare the resulting roots and their weight patterns, and use the study question to decide which solution is relevant. Multiple starts expose sensitivity; they do not by themselves prove which root is scientifically preferable.

Practical Reliability Checklist

Validate the start.
Check that the initial scatter is symmetric positive definite and numerically invertible. Use a defensible full-rank start rather than relying on an unstable inverse.

Check each reweighting step.
Confirm that distances and weights are finite, the total weight is positive, and positive-weight residuals span the full observation space.

Monitor both estimates.
Track relative changes in location and scatter, scatter eigenvalues, conditioning, and the distribution of weights. Apply damping only as an explicit stabilization choice.

Assess the returned root.
Check equation residuals and repeat from distinct starts. Report materially different converged roots or failed starts rather than concealing that sensitivity.

The central computational distinction is between a defined update, a numerically stable iteration, and a scientifically appropriate solution. Each is a separate claim and should be checked separately.

Pitfall: treating small changes as proof of a valid estimate. An iteration can appear stationary because weights have collapsed or the algorithm is trapped near an unintended root. Check equation residuals, positive definiteness, and the active weights as well as successive differences.
Pitfall: assuming a ridge term is harmless. Replacing a singular matrix by \(\Sigma+\lambda I\) makes inversion possible but changes distances and can change which observations are downweighted. Record the intervention and test sensitivity to it.
Pitfall: expecting every M-estimation iteration to decrease an objective. Weighted least squares minimizes a criterion with weights held fixed; it does not establish monotonicity after the weights change. Do not infer convergence from that minimization fact alone.

Check Your Understanding

Use the equations and diagnostics above to answer these questions.

Why must the scatter matrix be positive definite before the next Mahalanobis distances can be computed?
State the condition on positive-weight residuals that makes the weighted scatter update positive definite, and explain why it is necessary.
How can a nearly singular starting matrix alter robust weights even when its entries look reasonable?
What does damping guarantee when both scatter matrices being combined are positive definite, and what does it not guarantee?
Why should an analyst use multiple starts when the weight function is redescending?