Propagation with the Covariance Matrix

Statement

Given a random input vector x ∈ ℝⁿ with mean μ and covariance matrix Σ_x, and a differentiable vector-valued map y = f(x) with f: ℝⁿ → ℝ^m, the covariance of the output to first order is Σ_y = J Σ_x J^T, where J is the m×n Jacobian J_ij = ∂f_i/∂x_j evaluated at μ.

Why it matters

Every measured quantity carries uncertainty, and real experiments almost never measure the final quantity of interest directly — they measure several correlated inputs and combine them. The covariance matrix is the complete second-order description of that joint uncertainty, and the Jacobian law tells you exactly how it transforms when you push the inputs through any smooth function. It is the workhorse of least-squares fitting, Kalman filtering, GPS trilateration, particle-physics track reconstruction, and any calibration chain.

Crucially, it handles correlations. Naive scalar error propagation, applied component-by-component, silently assumes the inputs are independent and the outputs are uncorrelated. Both assumptions are usually false, and both can move an error bar by a large factor. The matrix form is the honest bookkeeping that keeps off-diagonal terms where they belong.

Assumptions

Derivation

Result

Limiting cases

• Scalar output (m = 1): J is a row vector g^T = (∂f/∂x₁, …) and the formula collapses to σ_y² = g^TΣ_xg = ∑_jk (∂f/∂x_j)(∂f/∂x_k)(Σ_x)_jk, the standard scalar propagation law with covariance terms.
• Uncorrelated inputs: Σ_x = diag(σ₁², …, σ_n²) and the scalar case reduces further to σ_y² = ∑_j (∂f/∂x_j)²σ_j², the familiar “add in quadrature” rule.
• Linear map f(x) = A x + b: the Jacobian is exactly A everywhere, so Σ_y = A Σ_x A^T holds without any approximation and for any input distribution.
• Change of coordinates / rotation: if J = R is orthogonal (R^TR = I), the transformation is a similarity that rotates the uncertainty ellipsoid; its eigenvalues (variances along principal axes) are preserved.
• Single input, single output: J = f′(μ) and σ_y = |f′(μ)| σ_x, the elementary one-variable rule.

Breaks when

Failure modes

• Dropping the off-diagonal terms of Σ_x. Treating correlated inputs as independent and adding in quadrature; for strongly anti-correlated inputs this can overestimate the true error by a large factor (a difference of two positively correlated quantities is less uncertain than quadrature suggests).
• Forgetting the induced output correlations. Reporting only the diagonal of Σ_y and then feeding those “independent” variances into a further calculation, discarding the off-diagonal information that the shared inputs created.
• Evaluating J at the wrong point. Using the measured sample value instead of the best estimate of the mean, or (worse) leaving J symbolic and never numerically evaluating it at μ.
• Transpose slip. Writing J^TΣ_xJ (which has the wrong dimensions n×n when m ≠ n) instead of JΣ_xJ^T.
• Mixing standard deviations and variances. Placing σ rather than σ² on the diagonal of Σ_x, or reporting √(diag Σ_y) as a variance.
• Unit inconsistency in J. Building the Jacobian from a formula written in one unit system while Σ_x is expressed in another (degrees vs radians is the perennial offender in angle derivatives).

Discussion

The formula is a statement about how quadratic forms transform under linear maps. The input uncertainty is geometrically an ellipsoid δx^TΣ_x⁻¹δx = const in input space; the Jacobian maps that ellipsoid, to first order, onto a new ellipsoid in output space, and JΣ_xJ^T is precisely the covariance whose inverse defines the image. That is why the construction always returns a symmetric positive-semidefinite matrix: for any vector v, v^T(JΣ_xJ^T)v = (J^Tv)^TΣ_x(J^Tv) ≥ 0 because Σ_x is itself positive-semidefinite. Positive-semidefiniteness (not strict positivity) is important: if J has a nontrivial null space — fewer effective outputs than inputs — some directions in output space carry zero variance, and Σ_y is genuinely singular.

The result generalises two prior results into one object. The scalar first-order error-propagation formula and the variance-of-linear-combinations rule are both special cases: the former is the single-output projection, the latter is the exact statement for the linear map f = A x. Seeing them as diagonal-and-scalar shadows of the matrix law is the conceptual payoff — you stop memorising three formulas and carry one.

Connections run deep in applied statistics. In linear least squares, the fitted parameters are θ̂ = (A^TW A)⁻¹A^TW y, a linear map of the data, so their covariance Σ_θ = (A^TW A)⁻¹ (with W = Σ_y⁻¹) is exactly this law applied to the design matrix. In the Kalman filter, the predict step is P⁻ = F P F^T + Q — the same congruence by the state-transition Jacobian, with process noise added. The extended Kalman filter is nothing more than this formula with F the Jacobian of a nonlinear dynamics evaluated at the current estimate, which is also where its well-known divergence under strong nonlinearity comes from.

A subtle and important point is that JΣ_xJ^T propagates only the second moment, and for a nonlinear f the output mean is not simply f(μ). The next term in the expansion contributes a deterministic bias E[y_i] ≈ f_i(μ) + ½ tr(H_i Σ_x), where H_i is the Hessian of the i-th component. Practitioners who need the mean to matching order, or who face appreciable curvature, abandon the linearisation for the unscented transform (propagate a deterministic set of sigma points and reconstruct mean and covariance) or Monte Carlo. These methods reproduce JΣ_xJ^T in the small-fluctuation limit but remain correct when it fails, which is exactly why they exist.

Common misconceptions. The formula does not require the inputs to be Gaussian — it is a statement about second moments alone and holds for any distribution with a finite covariance. What Gaussianity buys you is that, for a linear map, the output is again exactly Gaussian and so Σ_y is the complete description; for a nonlinear map even Gaussian inputs give non-Gaussian outputs. A second misconception is that uncorrelated inputs give uncorrelated outputs: sharing input variables through the rows of J generically produces off-diagonal terms in Σ_y.

Worked examples

Example 1 — Polar to Cartesian with correlated-looking output. A radar returns range r and bearing θ. Convert to Cartesian (x, y) = (r cosθ, r sinθ) and propagate the covariance. Take r = 100 m with σ_r = 2 m, θ = 30° with σ_θ = 0.02 rad, uncorrelated inputs.

Example 2 — Sum and difference of two correlated measurements. Two voltmeters read V₁ = 5.00 V and V₂ = 3.00 V, each with σ = 0.10 V, but they share a common reference so their readings are positively correlated with correlation coefficient ρ = 0.6. Compute the covariance of the sum s = V₁ + V₂ and difference d = V₁ − V₂.

Problems

1. (Warm-up) For y = x₁x₂ with x₁ = 4.0 ± 0.2, x₂ = 5.0 ± 0.3 (uncorrelated), find σ_y using the Jacobian law.
   Solution

   J = (∂y/∂x₁, ∂y/∂x₂) = (x₂, x₁) = (5.0, 4.0). With Σ_x = diag(0.04, 0.09): σ_y² = (5.0)²(0.04) + (4.0)²(0.09) = 1.00 + 1.44 = 2.44, so σ_y = 1.56. Since y = 20.0, the relative error 1.56/20.0 = 7.8% equals √((0.2/4)² + (0.3/5)²) = √(0.0025+0.0036) as expected for a product.
2. (Correlation matters) Repeat problem 1 but now with correlation ρ = −0.5 between x₁ and x₂. Does σ_y rise or fall?
   Solution

   Cov = ρσ₁σ₂ = (−0.5)(0.2)(0.3) = −0.030. The scalar law gains a cross term 2J₁J₂Cov = 2(5.0)(4.0)(−0.030) = −1.20. So σ_y² = 2.44 − 1.20 = 1.24, σ_y = 1.11. Negative correlation reduces the product's uncertainty from 1.56 to 1.11.
3. (Two outputs) Let (u, v) = (x + y, x²) with x = 2.0, y = 1.0, σ_x = 0.1, σ_y = 0.2, uncorrelated. Find the full 2×2 output covariance, including Cov(u,v).
   Solution

   J = ( 1   1 ; 2x   0 ) = ( 1   1 ; 4   0 ). Σ_x = diag(0.01, 0.04). Var(u) = 1·0.01 + 1·0.04 = 0.05. Var(v) = 4²·0.01 + 0 = 0.16. Cov(u,v) = (1)(4)(0.01) + (1)(0)(0.04) = 0.04. So Σ_y = ( 0.05   0.04 ; 0.04   0.16 ). Note the outputs are correlated (through the shared variable x) despite independent inputs; correlation coefficient = 0.04/√(0.05·0.16) = 0.45.
4. (Rotation preserves area) An uncertainty ellipse has Σ_x = diag(9, 1) (units cm²). Apply a 45° rotation R = (1/√2)( 1   −1 ; 1   1 ). Find Σ_y and verify its determinant (and hence the ellipse area) is unchanged.
   Solution

   Σ_y = RΣ_xR^T. Compute RΣ_x = (1/√2)( 9   −1 ; 9   1 ); then times R^T = (1/√2)( 1   1 ; −1   1 ): Σ_y = ½( 9+1   9−1 ; 9−1   9+1 ) = ( 5   4 ; 4   5 ). det Σ_y = 25 − 16 = 9 = det Σ_x (9·1). Rotation preserves the determinant because det R = 1, so det(RΣ_xR^T) = det Σ_x; the ellipse area π√det is invariant.
5. (Where it breaks) For y = x² with x = 0.0 ± 0.5, the Jacobian law predicts σ_y = 0. Explain why, and compute the true standard deviation of y assuming x is Gaussian.
   Solution

   J = 2x = 0 at x = 0, so the first-order estimate gives σ_y = |2x|σ_x = 0 — a singular Jacobian (stationary point), a documented failure regime. The true spread lives entirely in the discarded quadratic term. For x ~ N(0, σ²) with σ = 0.5: E[y] = E[x²] = σ² = 0.25; E[y²] = E[x⁴] = 3σ⁴ = 3(0.0625) = 0.1875; Var(y) = 0.1875 − 0.0625 = 0.125, so σ_y = 0.354 (= √2 σ²), definitely not zero. The lesson: at a stationary point of f, propagate with the second-order term or by direct computation.