Least Squares Adjustment for Control Networks in Python: Deterministic Gauss-Markov Implementation and Survey-Grade Validation

Least squares adjustment (LSA) is the sub-task within Algorithmic Math & Geodetic Workflows that reconciles redundant geodetic observations into a single, statistically defensible coordinate solution. When a control network carries more measurements than unknowns — the normal condition in cadastral and high-precision work — those surplus observations cannot all be satisfied exactly, so the adjustment distributes the inevitable misclosures according to each observation’s precision rather than discarding them. This guide implements that reconciliation in deterministic Python: a weighted Gauss-Markov estimator with explicit np.float64 arithmetic, Cholesky-to-SVD solver routing, a global chi-square test on the variance factor, and the standardized-residual screening and structured logging that survey-grade work demands under ISO 19111 coordinate reference system metadata expectations. Every residual it emits is traceable to an observation equation, and every numeric path is reproducible bit-for-bit across runs.

Figure — Cholesky-to-SVD solver routing with convergence gating.

The Gauss-Markov prerequisite: weighted normal equations #

The functional model states that each observation is a function of the unknown parameters; linearized about an approximate solution it becomes the Gauss-Markov observation model

where $\mathbf{v}$ is the residual vector, $\mathbf{A}$ is the design matrix of partial derivatives of each observation with respect to each parameter, $\mathbf{x}$ is the parameter-correction vector, and $\boldsymbol{\ell}$ is the misclosure (observed minus computed). Heterogeneous observations — slope distances, zenith angles, horizontal directions, levelled height differences, and GNSS baseline components — carry different uncertainty, so they are weighted by the inverse of their covariance, $\mathbf{W} = \mathbf{C}_\ell^{-1}$. Minimizing the weighted sum of squared residuals $\mathbf{v}^{\mathsf{T}}\mathbf{W}\mathbf{v}$ yields the normal equations and their solution

The dispersion of the estimated parameters follows from propagation of variance, scaled by the a-posteriori variance factor $\hat{\sigma}_0^2$ that absorbs any mismatch between the assumed and realized observation precision:

Here $n$ is the number of observation equations, $u$ the number of estimated parameters, and the redundancy $r = n-u$ is the degrees of freedom. A variance factor near $1.0$ confirms the stochastic model matches reality; values far from unity signal mis-weighted observations or an unmodeled datum shift. This same variance factor is the quantity that downstream error distribution modeling in Python propagates into the positional error ellipse of every derived point, so the adjustment treats $\hat{\sigma}_0^2$ and $\mathbf{N}^{-1}$ as non-negotiable outputs.

Step-by-step implementation #

Step 1 — Build the weight matrix from observation variances #

Precision enforcement begins at ingestion. Observation variances are converted to np.float64, validated against instrument manufacturer specifications and agency tolerance tables, and clamped away from zero so the weight matrix stays positive-definite during inversion.

import numpy as np
from numpy.typing import NDArray


def build_weight_matrix(
    variances: NDArray[np.float64],
    min_variance: float = 1e-12,
) -> NDArray[np.float64]:
    """Construct the diagonal weight matrix W = C_l^-1 from observation variances.

    Satisfies the Gauss-Markov stochastic model: ISO 19111 requires the
    uncertainty of every observation to be explicit, never assumed uniform.
    """
    var = np.asarray(variances, dtype=np.float64)   # enforce float64 at the boundary
    if np.any(var < 0):
        raise ValueError("Observation variances must be non-negative.")
    var = np.maximum(var, min_variance)             # clamp near-zero variances
    return np.diag(1.0 / var)

A full covariance matrix (correlated GNSS baseline components, for example) replaces the diagonal construction with a direct inverse of C_l; the rest of the pipeline is agnostic to how W was formed, provided it is symmetric positive-definite.

Step 2 — Solve the normal equations with deterministic fallback routing #

Cholesky factorization is the optimal route for well-conditioned, positive-definite normal equations. When rank deficiency, collinear control, or a poorly constrained network produces a near-singular $\mathbf{N}$, the solver routes to a singular-value-decomposition least-squares fallback. Routing is driven by an explicit condition-number threshold, not by catching an exception after the fact, so the decision is logged and reproducible. The right-hand side uses the current residual $(\boldsymbol{\ell} - \mathbf{A}\mathbf{x})$ so accumulated corrections converge instead of re-adding the full solution on every iteration — the form required when the functional model is genuinely non-linear.

from scipy import linalg
from typing import Dict
import logging

logger = logging.getLogger(__name__)


def solve_control_network_lsa(
    A: NDArray[np.float64],
    l: NDArray[np.float64],
    W: NDArray[np.float64],
    tol_x: float = 1e-9,
    tol_v: float = 1e-6,
    max_iter: int = 15,
    cond_threshold: float = 1e12,
) -> Dict[str, object]:
    """Deterministic weighted least squares adjustment with explicit float64
    precision, condition monitoring, and Cholesky-to-SVD fallback routing.

    Implements the Gauss-Markov model v = A x - l; aligns with ISO 19111 CRS
    metadata and EPSG parameter propagation for the derived covariance.
    """
    A = np.asarray(A, dtype=np.float64)
    l = np.asarray(l, dtype=np.float64)
    W = np.asarray(W, dtype=np.float64)

    n_obs, n_params = A.shape
    if n_obs <= n_params:
        raise ValueError("Network is underdetermined; redundancy required for validation.")

    x_corr = np.zeros(n_params, dtype=np.float64)
    v_prev = np.full(n_obs, np.nan, dtype=np.float64)
    iteration = 0
    converged = False
    cond_num = np.nan

    while iteration < max_iter:
        # Normal equations from the *current* residual so corrections accumulate.
        N = A.T @ W @ A
        U = A.T @ W @ (l - A @ x_corr)

        cond_num = float(np.linalg.cond(N))
        if cond_num > cond_threshold:
            logger.warning("cond(N)=%.2e exceeds threshold; routing to SVD fallback.", cond_num)
            x_step, *_ = linalg.lstsq(N, U, lapack_driver="gelsd")  # divide-and-conquer SVD
        else:
            try:
                cho = linalg.cho_factor(N, lower=True)              # primary Cholesky route
                x_step = linalg.cho_solve(cho, U)
            except linalg.LinAlgError as exc:
                logger.warning("Cholesky failed (%s); routing to SVD fallback.", exc)
                x_step, *_ = linalg.lstsq(N, U, lapack_driver="gelsd")

        x_corr += x_step
        v = A @ x_corr - l

        dx_max = float(np.max(np.abs(x_step)))
        dv_norm = float(np.linalg.norm(v))
        dv_change = abs(dv_norm - np.linalg.norm(v_prev)) if not np.isnan(v_prev[0]) else np.inf

        if dx_max < tol_x and dv_change < tol_v:
            converged = True
            break

        v_prev = v.copy()
        iteration += 1

    if not converged:
        logger.error("Max iterations reached without convergence; audit network geometry.")

    dof = n_obs - n_params
    sigma0_sq = float((v.T @ W @ v) / dof)
    try:
        cov_x = sigma0_sq * np.linalg.inv(N)
    except np.linalg.LinAlgError:
        cov_x = sigma0_sq * np.linalg.pinv(N, hermitian=True)

    return {
        "x_corrections": x_corr,
        "residuals": v,
        "sigma0_sq": sigma0_sq,
        "covariance_params": cov_x,
        "iterations": iteration,
        "converged": converged,
        "condition_number": cond_num,
        "dof": dof,
    }

Step 3 — Apply the global chi-square test on the variance factor #

Convergence alone does not certify a network. The global test compares the realized $\hat{\sigma}_0^2$ against its a-priori expectation of $1.0$. Because $\mathbf{v}^{\mathsf{T}}\mathbf{W}\mathbf{v} = r\,\hat{\sigma}_0^2$ is chi-square distributed with $r$ degrees of freedom, the statistic must fall inside the two-sided acceptance band at the chosen confidence level. Rejection above the upper bound points to over-optimistic weights or undetected blunders; rejection below the lower bound points to pessimistic weights.

from scipy import stats


def global_chi_square_test(sigma0_sq: float, dof: int, alpha: float = 0.05) -> Dict[str, object]:
    """Two-sided chi-square test of the a-posteriori variance factor against 1.0.

    Agency QA: a network that fails this test must not be committed to the
    cadastral store regardless of how small the individual residuals look.
    """
    test_statistic = dof * sigma0_sq                      # v^T W v under H0: sigma0=1
    lower = float(stats.chi2.ppf(alpha / 2.0, dof))
    upper = float(stats.chi2.ppf(1.0 - alpha / 2.0, dof))
    return {
        "test_statistic": test_statistic,
        "lower_bound": lower,
        "upper_bound": upper,
        "passed": bool(lower <= test_statistic <= upper),
    }

Step 4 — Screen local reliability with standardized residuals #

A network can pass the global test while hiding a single blunder. Local screening normalizes each residual by its own standard deviation, taken from the residual cofactor matrix $\mathbf{Q}_{vv} = \mathbf{W}^{-1} - \mathbf{A}\mathbf{N}^{-1}\mathbf{A}^{\mathsf{T}}$. Standardized residuals with $|w_i| > 3.0$ are flagged for manual review or automated down-weighting before re-adjustment.

def standardized_residuals(
    A: NDArray[np.float64],
    W: NDArray[np.float64],
    v: NDArray[np.float64],
    sigma0: float,
) -> NDArray[np.float64]:
    """Standardized (data-snooping) residuals w_i = v_i / (sigma0 * sqrt(q_vv_ii)).

    q_vv is the diagonal of the residual cofactor matrix; |w_i| > 3 marks a
    candidate outlier under Baarda data snooping at the 99.7% level.
    """
    A = np.asarray(A, dtype=np.float64)
    W = np.asarray(W, dtype=np.float64)
    cofactor_ll = np.linalg.inv(W)
    N = A.T @ W @ A
    cofactor_vv = cofactor_ll - A @ np.linalg.inv(N) @ A.T
    q_vv = np.clip(np.diag(cofactor_vv), 1e-15, None)     # guard tiny negative drift
    return v / (sigma0 * np.sqrt(q_vv))

Parameter and return-value reference #

The table lists every key input and output of the solver, its type, units, valid range, and cadastral significance. It scrolls horizontally on narrow viewports.

Worked example: a NAVD88 levelling network #

Consider a four-benchmark vertical control network referenced to NAVD88 (EPSG:5703). Benchmark A is held fixed at $H_A = 100.000$ m; the heights of B, C, and D are unknown. Five levelled height differences — including a closing loop and one redundant diagonal — give a redundancy of two. Each row of A carries $+1$ at the upgradient station and $-1$ at the downgradient station; the fixed benchmark moves to the misclosure vector. Datum context for the EPSG codes and the fixed-station treatment is covered in setting up high-precision coordinate reference systems.

H_A = 100.000  # fixed benchmark, NAVD88 orthometric height (EPSG:5703), metres

# Parameters x = [H_B, H_C, H_D]; one row per levelled run A->B, B->C, C->D, D->A, A->C.
A = np.array([
    [1.0,  0.0, 0.0],
    [-1.0, 1.0, 0.0],
    [0.0, -1.0, 1.0],
    [0.0,  0.0, -1.0],
    [0.0,  1.0, 0.0],
], dtype=np.float64)

# Misclosure: observed dh with the fixed benchmark folded in where it appears.
l = np.array([5.102 + H_A, 3.621, -2.548, -6.205 - H_A, 8.730 + H_A], dtype=np.float64)

# A-priori standard deviations (metres); the longer A->C run is weighted weaker.
sigmas = np.array([0.008, 0.008, 0.008, 0.008, 0.012], dtype=np.float64)
W = build_weight_matrix(sigmas ** 2)

result = solve_control_network_lsa(A, l, W)
print("adjusted heights:", np.round(result["x_corrections"], 4))
print("sigma0^2:", round(result["sigma0_sq"], 4), "| dof:", result["dof"])

Running this prints the adjusted heights and the variance factor:

adjusted heights: [105.1083 108.7355 106.1963]
sigma0^2: 1.9117 | dof: 2

The adjusted height of B is $105.1083$ m against a single-run expectation of $H_A + 5.102 = 105.102$ m — the 6 mm shift is the adjustment distributing the loop misclosure across the redundant observations. Comparing the adjusted heights to an independent check benchmark is the acceptance step formalized in validating datum alignment with control points, and it is what separates a converged number from a defensible one.

Verification and residual analysis #

The verification stage runs the global test, screens local reliability, and serializes a structured audit record. The chi-square statistic of $3.82$ falls inside the 95% band $[0.051,\,7.378]$, so the network passes globally; the largest standardized residual is $1.38$, comfortably under the $3.0$ data-snooping limit, so no single observation dominates.

import json

sigma0 = result["sigma0_sq"] ** 0.5
chi = global_chi_square_test(result["sigma0_sq"], result["dof"])
w = standardized_residuals(A, W, result["residuals"], sigma0)

audit = {
    "converged": result["converged"],
    "iterations": result["iterations"],
    "condition_number": round(result["condition_number"], 3),
    "sigma0_sq": round(result["sigma0_sq"], 4),
    "dof": result["dof"],
    "chi_square": {k: (round(v, 3) if isinstance(v, float) else v) for k, v in chi.items()},
    "max_std_residual": round(float(np.max(np.abs(w))), 3),
    "worst_observation": int(np.argmax(np.abs(w))),
    "param_std_m": np.round(np.sqrt(np.diag(result["covariance_params"])), 4).tolist(),
}
logger.info("lsa audit: %s", json.dumps(audit))

# The gate: a network that fails the global test or hides a blunder is not committed.
assert chi["passed"], "Global chi-square test failed; review weights and blunders."
assert np.max(np.abs(w)) < 3.0, "Standardized residual exceeds 3-sigma; flag observation."

The two assert statements are the commit gate. A failure here stops the pipeline rather than writing a non-compliant coordinate to the cadastral store. Every adjustment run should emit, at minimum, the global chi-square statistic and p-value, the maximum standardized residual and its observation ID, the condition number of the normal matrix, the final $\hat{\sigma}_0^2$ with its degrees of freedom, and the propagated coordinate covariance. The numeric thresholds themselves — confidence level, condition limit, residual band — belong in one place; centralizing them is the subject of tuning transformation thresholds for survey-grade work.

Troubleshooting and gotchas #

Singular or near-singular normal matrix. Collinear or tightly clustered control inflates cond(N) and corrupts $\mathbf{N}^{-1}$. The cond_threshold gate routes to the SVD fallback, but a fallback solution on a rank-deficient network has an arbitrary component in its null space — add geometrically independent control rather than lowering the threshold to mask the warning.

Single- vs double-precision promotion. Passing float32 arrays (common from binary I/O paths) silently degrades the variance factor and can flip small eigenvalues of $\mathbf{N}$ negative during factorization, breaking Cholesky. Every entry point casts with np.asarray(..., dtype=np.float64); never strip those casts for speed.

Mis-scaled observation covariance. If variances arrive in the wrong units — millimetres² instead of metres² — the solver still converges but $\hat{\sigma}_0^2$ lands orders of magnitude from $1.0$ and the global test rejects. Treat any variance factor outside roughly $0.1$–$10$ as a stochastic-model error, not a measurement outlier, and check units before re-weighting.

Re-adding the full solution each iteration. Forming the right-hand side as $\mathbf{A}^{\mathsf{T}}\mathbf{W}\boldsymbol{\ell}$ instead of $\mathbf{A}^{\mathsf{T}}\mathbf{W}(\boldsymbol{\ell} - \mathbf{A}\mathbf{x})$ makes a non-linear network oscillate or diverge because each step re-applies corrections already absorbed. The implementation above uses the current residual; keep it that way for trilateration and resection models that need real iteration.

Zero redundancy. With exactly the minimum observations ($n = u$) the redundancy is zero, $\hat{\sigma}_0^2$ is undefined, and the constructor raises. A network with no surplus cannot be validated statistically — add redundant ties before treating the covariance as defensible.

Untreated systematic distortion. When standardized residuals show a coherent spatial pattern rather than random scatter, the cause is usually unmodeled regional deformation, not blunders. Removing it with polynomial shift algorithms for regional adjustments before re-adjusting restores a variance factor near unity; for purely local distortion, implementing affine transformations for local grids is the lighter-weight pre-step.

Related guides #

• Algorithmic Math & Geodetic Workflows — parent guide to the deterministic transformation pipeline
• Error distribution modeling in Python — propagating this adjustment’s variance factor into positional error ellipses
• Polynomial shift algorithms for regional adjustments — removing spatially correlated distortion before re-adjustment
• Implementing affine transformations for local grids — the linear pre-step that minimizes misclosure on local networks
• Validating datum alignment with control points — comparing adjusted coordinates against independent control
• Tuning transformation thresholds for survey-grade work — centralizing the confidence, condition, and residual limits used above