Detecting Outliers with Normalized Residuals

Once a least-squares control-network adjustment has converged, the next obligation is to prove that no single observation is a blunder corrupting the whole solution — the disciplined way to do that is Baarda data snooping on normalized residuals, a technique that sits directly under least-squares adjustment for control networks within the broader Algorithmic Math & Geodetic Workflows reference. A raw residual of five millimetres means nothing on its own; it only becomes actionable once it is scaled by how well the network geometry could actually detect an error at that observation.

Figure — data snooping rejects at most one observation per pass, then re-adjusts.

What Normalizing the Residual Actually Does #

Least-squares distributes an error across many observations, so a gross blunder in one measurement leaks into the residuals of its well-connected neighbours. The residual vector $\mathbf{v}$ therefore does not point cleanly at the culprit, and its entries do not share a common scale — a residual at a high-redundancy observation is far more diagnostic than one at a barely-checked observation. Baarda’s insight was to divide each residual by its own standard deviation so that, under the null hypothesis of no blunder, every scaled residual follows the same standard normal distribution. The normalized (studentized) residual is

where $v_i$ is the $i$-th residual, $\sigma_0$ is the a-priori reference standard deviation (unit weight), and $q_{v_i v_i}$ is the $i$-th diagonal entry of the cofactor matrix of the residuals, $\mathbf{Q}_{vv} = \mathbf{Q}_{\ell\ell} - \mathbf{A}\,\mathbf{N}^{-1}\mathbf{A}^{\mathsf T}$. That diagonal is the local redundancy number carried in cofactor form: it shrinks toward zero where the network cannot check an observation, which is exactly why raw residuals mislead. Each $w_i$ is then compared against a critical value $c$ — commonly $3.29$ for a per-observation significance of $\alpha = 0.001$ — and the observation with the largest $|w_i|$ above $c$ is flagged as the most likely single blunder. This is the same variance-scaling discipline that underpins error distribution modeling in Python, where residuals are propagated through the full covariance rather than judged in raw metres.

Complete Runnable Implementation #

The routine below takes the residual vector, the residual cofactor matrix $\mathbf{Q}_{vv}$, and the reference standard deviation, then returns each normalized residual and a single-blunder flag following the one-at-a-time data-snooping rule. It never mutates its inputs and clamps near-zero cofactor diagonals so an unchecked observation cannot masquerade as a huge normalized residual.

from __future__ import annotations

import logging
from dataclasses import dataclass

import numpy as np

logger = logging.getLogger("data_snooping")

# Baarda critical values for common per-observation significance levels.
BAARDA_CRITICAL: dict[float, float] = {0.05: 1.96, 0.01: 2.58, 0.001: 3.29}


@dataclass(frozen=True)
class SnoopingResult:
    normalized_residuals: np.ndarray  # w_i, one per observation
    flagged_index: int                # index of the worst blunder, -1 if none
    flagged_w: float                  # |w| at that index (0.0 if none)
    critical_value: float             # c used for the test


def detect_outlier(
    residuals: np.ndarray,
    q_vv: np.ndarray,
    sigma0: float = 1.0,
    alpha: float = 0.001,
    zero_floor: float = 1e-12,
) -> SnoopingResult:
    """Flag the single most likely blunder by Baarda data snooping.

    Computes the normalized residual w_i = v_i / (sigma0 * sqrt(q_vv,i))
    for every observation and returns the index of the largest |w_i| that
    exceeds the critical value, or -1 when the network is clean.

    Args:
        residuals: length-n residual vector v from the adjustment (metres).
        q_vv: n-by-n residual cofactor matrix Q_vv, or its length-n diagonal.
        sigma0: a-priori reference standard deviation of unit weight.
        alpha: per-observation significance level (0.05, 0.01, or 0.001).
        zero_floor: cofactor-diagonal values at or below this are treated as
            uncheckable (redundancy ~ 0) and excluded from the test.

    Returns:
        A SnoopingResult with the full w vector and the worst-blunder flag.
    """
    v = np.asarray(residuals, dtype=np.float64).ravel()
    q = np.asarray(q_vv, dtype=np.float64)
    q_diag = np.diag(q) if q.ndim == 2 else q.ravel()
    if q_diag.shape != v.shape:
        raise ValueError(
            f"residuals ({v.shape}) and q_vv diagonal ({q_diag.shape}) disagree"
        )
    if sigma0 <= 0.0:
        raise ValueError("sigma0 must be positive")

    checkable = q_diag > zero_floor
    w = np.zeros_like(v)
    denom = sigma0 * np.sqrt(np.where(checkable, q_diag, 1.0))
    w[checkable] = v[checkable] / denom[checkable]

    critical = BAARDA_CRITICAL.get(alpha)
    if critical is None:
        raise ValueError(f"no tabulated critical value for alpha={alpha}")

    abs_w = np.abs(w)
    worst = int(np.argmax(abs_w)) if abs_w.size else -1
    if worst < 0 or abs_w[worst] <= critical:
        return SnoopingResult(w, -1, 0.0, critical)
    logger.info("blunder at obs %d: |w|=%.3f > c=%.2f", worst, abs_w[worst], critical)
    return SnoopingResult(w, worst, float(abs_w[worst]), critical)

Parameter Reference #

Worked Example #

A five-observation levelling loop with one deliberately corrupted measurement makes the scaling visible: the blunder does not have the largest raw residual, yet it has the largest normalized residual.

import numpy as np

v = np.array([0.004, -0.003, 0.021, 0.005, -0.002])          # metres
q_vv = np.array([0.62, 0.55, 0.71, 0.10, 0.48])              # cofactor diagonal
res = detect_outlier(v, q_vv, sigma0=0.005, alpha=0.001)
print(np.round(res.normalized_residuals, 2))
print(res.flagged_index, round(res.flagged_w, 2))
# [ 1.02 -0.81  4.98  3.16 -0.58]
# 2 4.98

Validation Check #

Gate the flag before it is used to reject an observation, and re-adjust after every single removal rather than dropping several at once.

assert res.flagged_index == 2, "expected the corrupted observation to be flagged"
assert res.flagged_w > res.critical_value, "flagged w must exceed the critical value"
assert np.isfinite(res.normalized_residuals).all(), "non-finite normalized residual"

Common Mistakes #

Related #

• Least-squares adjustment for control networks — the parent adjustment that produces the residuals and cofactor tested here.
• Error distribution modeling in Python — how the surviving residuals are propagated through the full covariance.
• Python least-squares adjustment for control points — the affine fit whose residuals are screened the same way.
• Algorithmic Math & Geodetic Workflows — the pipeline that adjustment, snooping, and covariance propagation feed.