Multivariate Methods

Consider a person who is 6'8" tall and weighs 130 pounds. Their height is unusual but not extreme—professional basketball players routinely exceed this. Their weight is on the low side but not alarming—some people are naturally thin. Yet the combination of extreme height with low weight is highly anomalous—it suggests a medical condition, data entry error, or fraud.

This illustrates a fundamental limitation of univariate methods: anomalies often hide in the relationships between variables, not in individual variables. The Z-score, IQR, and Grubbs' tests we've studied can miss these multivariate outliers entirely.

Multivariate anomaly detection addresses this by considering the joint distribution of all variables simultaneously, accounting for correlations and dependencies that define what constitutes 'normal' in high-dimensional space.

Let's formalize why applying univariate tests to each variable independently can fail dramatically.

The Marginal-Joint Discrepancy

Consider a 2D dataset where variables $X$ and $Y$ are positively correlated (when $X$ is high, $Y$ tends to be high). The data forms an elliptical cloud in 2D space.

A point with moderately high $X$ and moderately low $Y$ might be:

• Univariately normal: Both $X$ and $Y$ values are within 2 standard deviations of their respective means
• Multivariately anomalous: The combination violates the correlation structure—the point lies far from the data cloud

Mathematically: If $X$ and $Y$ are marginally normal, applying Z-scores to each gives:

• $z_X = \frac{x - \mu_X}{\sigma_X}$
• $z_Y = \frac{y - \mu_Y}{\sigma_Y}$

Even if $|z_X| < 3$ and $|z_Y| < 3$, the point $(x, y)$ can be highly improbable under the joint distribution.

Geometric Intuition

Univariate thresholds define a rectangle (in 2D) or hypercube (in higher dimensions) aligned with the coordinate axes.

Joint probability contours for correlated normal data form an ellipse (in 2D) or ellipsoid (in higher dimensions) that can be rotated and stretched.

The rectangle and ellipse can differ dramatically:

• Points inside the rectangle but outside the ellipse are multivariate outliers missed by univariate tests
• Points outside the rectangle but inside the ellipse are false positives from univariate tests

As dimensionality increases, this discrepancy worsens exponentially.

The Curse of Dimensionality Preview

In high dimensions, most of the probability mass in a Gaussian distribution lies in a thin shell at a specific distance from the center—not at the center itself. Univariate reasoning fails to capture this phenomenon. Points that appear 'central' in each dimension can be extremely unlikely in combination.

The Mahalanobis distance, introduced by Prasanta Chandra Mahalanobis in 1936, is a distance measure that accounts for correlations between variables. It is the fundamental tool for multivariate anomaly detection.

Definition

For a point $\mathbf{x} \in \mathbb{R}^p$ and a distribution with mean vector $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$, the Mahalanobis distance is:

$$D_M(\mathbf{x}) = \sqrt{(\mathbf{x} - \boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})}$$

Often, we work with the squared Mahalanobis distance:

$$D_M^2(\mathbf{x}) = (\mathbf{x} - \boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})$$

Understanding the Formula

Let's unpack this expression:

1. $(\mathbf{x} - \boldsymbol{\mu})$: Centers the point relative to the mean (just like Z-score)

2. $\boldsymbol{\Sigma}^{-1}$: The inverse covariance matrix (precision matrix) standardizes and decorrelates

3. The quadratic form: Computes a weighted sum of squared deviations, where weights depend on variances and covariances

Special Cases

Case 1: Uncorrelated Variables

If variables are uncorrelated, $\boldsymbol{\Sigma}$ is diagonal: $$\boldsymbol{\Sigma} = \text{diag}(\sigma_1^2, \sigma_2^2, \ldots, \sigma_p^2)$$

The Mahalanobis distance simplifies to: $$D_M^2(\mathbf{x}) = \sum_{i=1}^{p} \left(\frac{x_i - \mu_i}{\sigma_i}\right)^2 = \sum_{i=1}^{p} z_i^2$$

This is the sum of squared Z-scores—the Euclidean distance in standardized space.

Case 2: Identity Covariance

If $\boldsymbol{\Sigma} = \mathbf{I}$ (standardized uncorrelated data), Mahalanobis distance equals the standard Euclidean distance from the mean: $$D_M(\mathbf{x}) = |\mathbf{x} - \boldsymbol{\mu}|_2$$

Geometric Interpretation

The covariance matrix $\boldsymbol{\Sigma}$ can be decomposed as $\boldsymbol{\Sigma} = \mathbf{V} \mathbf{D} \mathbf{V}^T$, where:

• $\mathbf{V}$ contains the eigenvectors (principal directions)
• $\mathbf{D}$ contains the eigenvalues (variances along principal directions)

Mahalanobis distance effectively:

1. Rotates the data so principal components align with coordinate axes (using $\mathbf{V}^T$)
2. Scales each rotated dimension by its standard deviation (using $\mathbf{D}^{-1/2}$)
3. Computes Euclidean distance in this transformed space

This transforms the elliptical probability contours into circular ones, making simple distance comparison meaningful.

Just as Z-scores have probabilistic interpretation under normality, Mahalanobis distances have a precise distributional characterization.

The Chi-Squared Connection

If $\mathbf{x} \sim \mathcal{N}_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})$ (multivariate normal), then the squared Mahalanobis distance follows a chi-squared distribution with $p$ degrees of freedom:

$$D_M^2(\mathbf{x}) \sim \chi^2_p$$

This result is fundamental—it allows us to set thresholds with known false positive rates.

Setting Thresholds

To achieve a false positive rate of $\alpha$, set the threshold at the $(1-\alpha)$ quantile of the chi-squared distribution:

$$\tau = \chi^2_{p, 1-\alpha}$$

Points with $D_M^2(\mathbf{x}) > \tau$ are flagged as outliers.

Interpretation

Notice how the threshold increases with dimensionality. In 2D, $D_M^2 > 5.99$ is an outlier at α=0.05. In 100D, you need $D_M^2 > 124$ for the same significance.

Why? In high dimensions, even normal points are expected to be 'far' from the center in Mahalanobis distance. The chi-squared distribution with more degrees of freedom has a larger mean and wider spread.

Mean of $\chi^2_p$: $\mathbb{E}[\chi^2_p] = p$

Variance of $\chi^2_p$: $\text{Var}(\chi^2_p) = 2p$

So in 100 dimensions, the average squared Mahalanobis distance is 100, not near zero.

Algorithm: Mahalanobis Distance-Based Outlier Detection

Input: Dataset $\mathbf{X} \in \mathbb{R}^{n \times p}$ (n observations, p variables), significance level $\alpha$

Output: Boolean array indicating outliers

1. Estimate mean vector: μ̂ = (1/n) Σᵢ xᵢ
2. Estimate covariance matrix: Σ̂ = (1/(n-1)) Σᵢ (xᵢ - μ̂)(xᵢ - μ̂)ᵀ
3. Compute precision matrix: Σ̂⁻¹ (with regularization if needed)
4. For each observation i:
   D²ᵢ = (xᵢ - μ̂)ᵀ Σ̂⁻¹ (xᵢ - μ̂)
5. Compute threshold: τ = χ²_p,1-α
6. Flag outliers: D²ᵢ > τ

Computational Complexity:

• Covariance estimation: $O(np^2)$
• Matrix inversion: $O(p^3)$
• Distance computation per point: $O(p^2)$
• Total: $O(np^2 + p^3)$

The masking problem we encountered with univariate methods becomes even more severe in the multivariate setting.

How Outliers Corrupt Covariance Estimates

When computing the sample covariance matrix, outliers can:

1. Inflate variances: Extreme values stretch the ellipsoid along certain axes
2. Distort correlations: Outliers can create spurious correlations or mask true ones
3. Shift the mean: Multiple outliers in one direction bias the center

The combined effect: The estimated ellipsoid expands and distorts to 'accommodate' the outliers, making their Mahalanobis distances appear smaller than they should be.

Illustration

Consider a clean dataset forming a tight ellipse. Add a single extreme outlier:

• Without the outlier: Mahalanobis distance of outlier position would be huge
• With the outlier in the data: The covariance estimate expands toward the outlier, so its computed Mahalanobis distance is much smaller

This is why directly applying Mahalanobis distance with sample estimates can miss exactly the outliers we're trying to detect.

Breakdown Point Implications

The sample mean has breakdown point 1/n. The sample covariance matrix has breakdown point approximately 1/n as well for certain types of contamination. In p dimensions, a single outlier can arbitrarily distort the covariance estimate in p(p+1)/2 ways (the number of unique covariance parameters).

To address the masking problem, we need robust estimators of location and scatter that are not influenced by outliers.

The Minimum Covariance Determinant (MCD) Estimator

The MCD estimator, introduced by Rousseeuw (1984), finds the subset of observations that gives the smallest determinant of the sample covariance matrix.

Definition:

Given $n$ observations in $\mathbb{R}^p$, find the subset $H$ of size $h$ (where $h > (n+p+1)/2$) that minimizes:

$$\det(\hat{\boldsymbol{\Sigma}}_H)$$

Where $\hat{\boldsymbol{\Sigma}}_H$ is the sample covariance of the subset $H$.

Intuition: The subset with smallest covariance determinant is the 'tightest' cluster—the core of non-outlying observations. Outliers, by definition, lie far from this core and would inflate the determinant if included.

Properties of MCD

High breakdown point: If $h \approx n/2$, the MCD can resist up to ~50% contamination.

Affine equivariance: The estimator behaves correctly under linear transformations of the data.

Statistical efficiency: For normal data, MCD is less efficient than classical estimates, but this is the price of robustness.

When the number of features $p$ approaches or exceeds the number of observations $n$, classical covariance estimation breaks down entirely.

The Singularity Problem

The sample covariance matrix requires $n > p$ to be non-singular (invertible). When $n \leq p$, the matrix is rank-deficient and cannot be inverted directly.

Even when $n > p$, if $n$ is not substantially larger than $p$, the estimates are highly unstable. A common rule of thumb is $n \geq 5p$ to $10p$ for reliable estimation.

Solutions for High Dimensions

1. Regularization

Add a regularization term to the covariance matrix to ensure invertibility: $$\hat{\boldsymbol{\Sigma}}_{\text{reg}} = \hat{\boldsymbol{\Sigma}} + \lambda \mathbf{I}$$

Where $\lambda > 0$ is a regularization parameter. This shrinks eigenvalues toward a common value.

2. Shrinkage Estimators

Ledoit-Wolf shrinkage interpolates between the sample covariance and a structured target (e.g., diagonal or scaled identity): $$\hat{\boldsymbol{\Sigma}}_{\text{shrunk}} = (1 - \alpha) \hat{\boldsymbol{\Sigma}} + \alpha \mathbf{T}$$

Where $\alpha$ is optimally chosen to minimize estimation error.

3. Dimensionality Reduction First

Apply PCA to reduce dimensions to $k << n$, then perform Mahalanobis-based detection in the reduced space. The first $k$ principal components often capture the main structure; anomalies may appear in minor components too.