Adaptive & Robust
Kalman Filtering

Chapter 1: What Q and R Actually Encode

Before we can adapt Q and R, we need to feel what they DO. Not their definitions — kalman-filter.html gave you those — but their fingerprints: what exactly goes wrong, in which direction, when each one lies.

Reframe both matrices as confessions rather than parameters. Q says: "between two ticks, the world does things my model F did not predict — this much of them." It is the error bar on your physics. R says: "when the sensor speaks, this is the variance of its lie." It is the error bar on your evidence. The filter's entire behavior — every gain it ever computes — is the argument between these two confessions.

Let's make the argument quantitative with the simplest possible system: a random-walk state (F = 1) with Var(w) = q and Var(v) = r. The predicted variance obeys P_pred = P + q with the posterior P = P_pred r / (P_pred + r). At steady state these two lines close into a quadratic, and the quadratic has a closed-form solution:

Every claim in this chapter is this one formula, evaluated. Read its regimes: for q much smaller than r, P_pred ≈ √(qr), so K ≈ √(q/r). The gain depends on the RATIO of the two confessions — and only as a square root. Cut q by 100 and K falls by 10, not 100. Hold that thought: it has a sharp consequence for adaptation later, because steady-state behavior cannot identify q and r separately — only their ratio. Chapter 4 pays this bill.

Worked Numbers — the Gain Law, Evaluated

Take r = 1 and the two q values from chapter 0's disaster.

q = 0.01: P_pred = (0.01 + √(0.0001 + 0.04))/2 = (0.01 + 0.200250)/2 = 0.105125. Gain K = 0.105125/1.105125 = 0.095125. Response time 1/K = 10.5 steps.

q = 0.0001 (the chapter-0 mistuning): P_pred = (0.0001 + √(1e-8 + 4e-4))/2 = 0.010050, K = 0.010050/1.010050 = 0.009950. Response time 1/K = 100.5 steps.

The ratio of gains is 0.095125/0.009950 = 9.56 — cutting q by 100 cut K by roughly √100 = 10, confirming K ≈ √(q/r) in the small-q regime (√0.01 = 0.1 and √0.0001 = 0.01, both within 5% of the exact gains). Sanity-check the closed form by brute force: start P = 10 and iterate the two-line recursion 10,000 times with q = 0.0001; the iterated P_pred lands on 0.010050 — the same number, so the quadratic solution is the true fixed point.

The Four-Quadrant Symptom Table

Each mistuning direction has a distinct, diagnosable symptom — this table is the diagnostic heart of the chapter:

The diagnostic: Q errors and R errors both change smoothness, but their INNOVATION signatures differ — which is exactly what chapter 2 learns to read.

One more confession, about the confessions. Engineers commit a deliberate crime with Q: they use it as a tuning knob rather than a physical quantity. Real systems inflate Q to cover unmodeled dynamics — linearization error in the EKF, unmodeled accelerations in constant-velocity trackers — a lie told in a safe direction. The steady-state formula is the evidence for why the direction is safe: overestimating q costs you variance like a square root (gentle), underestimating costs you divergence (catastrophic). When forced to guess, guess fat.

Close with the identifiability foreshadow. Since K depends only on q/r in steady state, two engineers — one with the true (Q, R), one with both doubled — ship IDENTICAL estimate sequences. But NOT identical covariances: the computed P and the promised innovation variance S = P_pred + r both scale up with the common factor. The estimates cannot tell the difference; the innovation STATISTICS can. That asymmetry is the crack that chapters 2–3 wedge open.

Chapter 2: The Innovation — the Filter's Only Mirror

The filter cannot see its own error — that would require the truth. But there is one quantity it predicts BEFORE nature answers: the next measurement. The gap between prediction and answer is the innovation, and it is the only mirror the filter will ever get.

Define the object: ν_k = z_k − H x̂_pred is the part of the measurement the filter could NOT predict — the genuinely new information, hence the name. Everything predictable was already in x̂_pred; what remains should be pure surprise. You met innovations in sf-21-debugging-consistency.html as a way to TEST a filter on the bench. Here is the pivot this lesson needs: the innovation becomes the sensor that DRIVES adaptation, in flight.

A correctly-modeled filter makes a precise, three-part promise about its innovation stream:

Zero mean: no systematic surprise. Covariance exactly S: the filter knows HOW surprised it should be. And whiteness: innovations at different times are uncorrelated, because an optimal filter extracts every predictable pattern, leaving residue with no structure. The derivations live in sf-21; what we need here is what each violation MEANS.

Whiteness is the deep one, and the intuition is worth a PhD: if innovations are correlated in time, there is a pattern in the surprises — and a pattern in the surprises is information the filter failed to extract. A too-small gain leaves tracking error that persists across steps, so consecutive innovations share the same sign: positive autocorrelation. This is Mehra's classic 1970 observation, and it is directional evidence: fat-but-white innovations indict the noise levels; colored innovations indict the gain — the q/r ratio, or the model itself.

Worked Numbers — Auditing the Sleeping Filter

Quantify chapter 0's disaster. Truth: random walk with q_true = 0.25, measured with r_true = 1. The filter assumes q = 0.0025 — 100x too small.

Step 1 — the promise. The filter's internal steady state under its assumed q: P_pred = (0.0025 + √(0.0025² + 4×0.0025))/2 = 0.0513. So it PROMISES S = P_pred + r = 1.0513 and runs gain K = 0.0488.

Step 2 — the reality. For a fixed suboptimal gain the actual error dynamics are e_post = (1−K)e_pred + Kv with e_pred picking up a fresh w each step, which closes into a steady state P_actual = (K²r_true + q_true)/(1 − (1−K)²) = 2.6523.

Step 3 — the broken promise. The ACTUAL innovation variance is C = P_actual + r_true = 3.6523. The filter promised 1.0513 and delivers 3.6523 — a broken promise by the factor C/S = 3.47. That ratio is exactly what a running NIS average would read.

Step 4 — Monte Carlo confirms. A 200,000-step seeded run measures innovation variance 3.634, and lag-1 autocorrelation 0.674 — the surprises are massively patterned, because the sleepy gain leaves persistent tracking error that leaks into consecutive innovations.

Step 5 — the control. Same data, filter assuming the TRUE q = 0.25: promised S = 1.6404, measured 1.648, lag-1 autocorrelation −0.002. Right scale, white residue. The mirror works.

The promise hands you a free gift on the way out: the innovation gate. Since ν should be N(0, S), a measurement with ν²/S > 9 (a 3-σ event) is a once-in-370 occurrence; reject it as an outlier before it touches the state. This is the cheapest robustness in all of estimation — and it is only as good as S is honest. An overconfident filter (small S) starts gating TRUE measurements, which is how the confident-lie failure cascades. Depth on outlier theory — influence functions, Huber losses, when hard gates fail — belongs to robust-estimation.html; the gate itself belongs here, because S does.

Chapter 3: Innovation-Based Adaptation — Closing the Loop

Chapter 2 built the diagnosis: measured surprise versus promised surprise. The adaptive move is almost embarrassingly direct — if the promise does not match the measurement, change the promise until it does.

Set up the estimator. Over a sliding window of the last W innovations, form the sample covariance Ĉ = (1/W) Σ νν^T. This is the measured side. The promised side is S = H P_pred H^T + R. The innovation identity says that IF the model were right, these would agree in expectation. So read the identity backwards, as an equation FOR R:

Everything on the right is observable. Nothing needs truth. This is innovation-based adaptive estimation (IAE), the correlation-family member of Mehra's classic taxonomy of adaptive filters — Bayesian, maximum-likelihood, correlation-matching, covariance-matching (chapter 4 takes the last).

Now the subtlety that separates this from a homework exercise: the P_pred in the formula is the FILTER'S internal covariance, computed under its current wrong R. So the first solve gives an IMPROVED R, not the true one. The equation is implicit: changing R changes the gain, which changes P_pred, which changes what R̂ solves to. The honest algorithm is a fixed-point iteration — and the fixed point is exactly the truth, because when assumed R equals true R the filter is consistent and C = H P_pred H^T + R holds with equality.

Worked Numbers — Walking the Fixed Point

Scalar random walk with q = 0.09 (known and correct) and true measurement variance r = 4.0. The filter starts believing r = 1.0. At each iteration: compute the filter's internal steady-state P_pred and gain under its assumed r, compute the ACTUAL innovation variance C that gain produces (closed form, exactly as in chapter 2), then update via R̂ = C − P_pred.

Iteration 1: assumed r = 1.0 gives internal P_pred = 0.3484, gain K = 0.2584. Actual innovation variance: C = P_actual + r_true where P_actual = (K²r_true + q)/(1−(1−K)²), giving C = 4.7934. Update: R̂ = 4.7934 − 0.3484 = 4.4450. An overshoot — the overconfident filter underestimates its own P_pred, so too much of C gets attributed to R.

Iteration 2: assumed r = 4.4450 gives P_pred = 0.6791 and C = 4.6475, so R̂ = 3.9684. Now slightly under — the loop is contracting.

Iterations 3–5: R̂ = 4.0024, then 3.9998, then 4.0000. Geometric convergence onto the truth.

Fixed-point check: at assumed r = 4.0 exactly, P_pred = 0.6467, the filter is consistent so C = P_pred + 4.0 = 4.6467, and R̂ = C − P_pred = 4.0 — the update maps the truth to itself.

Window size W is the knob. W = 30–100 is typical. Smaller tracks sensor changes faster but injects estimation noise into the gain: the variance of a sample variance of Gaussian data is 2σ⁴/W, so R̂'s relative jitter is roughly √(2/W) — about 20% at W = 50, and that jitter feeds straight into K. Larger W is calmer but slower to notice a degrading sensor. You are spending one statistical budget two ways.

What about Q? Same door, harder passage. The innovation identity constrains H P_pred H^T + R as a SUM — with R trusted, surplus innovation variance can be attributed to Q: raise q until S matches Ĉ, using P_pred's dependence on q from chapter 1's Riccati closed form. But H projects the state: innovations only see Q through the measured subspace, so process noise on unobserved states is weakly identified, slowly, through the dynamics. Practical systems therefore adapt a scalar inflation on a fixed Q shape rather than a full matrix — estimate the level, trust the structure.

And the boundary every adaptive filter must fear: adapt R or Q from the same innovation stream and you are fine; adapt both freely and you have one equation (the S level) with two unknowns — chapter 1's ratio identifiability, now with teeth. The standard discipline: fix whichever side you can characterize offline (R from sensor bench data, usually), adapt the other online, and let whiteness — not just scale — arbitrate when both are suspect. Colored innovations indict the Q side; white-but-fat indicts the R side. Chapter 4 shows what happens when this discipline is skipped.

From Scratch: the Complete IAE Loop

The whole chapter in 40 lines. A mistuned filter (assumed R four times too small) versus the same filter with sliding-window adaptation, on the same seeded data. The frozen filter's mean NIS reads ~3.65 — chapter 2's broken-promise ratio; the adaptive one heals to ~1.03 with R̂ ≈ 3.85 against a true 4.0. No ground truth used anywhere.

python
import numpy as np

# Scalar random-walk tracker: q known and correct, R badly wrong.
rng = np.random.default_rng(11)
q_true, r_true = 0.09, 4.0
N = 4000
xs = np.cumsum(rng.normal(0, np.sqrt(q_true), N))   # truth
zs = xs + rng.normal(0, np.sqrt(r_true), N)          # measurements

def run(adapt, W=50, r_assumed=1.0):
    """One filter pass. adapt=True closes the IAE loop:
    R_hat = mean(nu^2 over window) - P_pred  (innovation identity,
    read backwards), floored to stay positive."""
    xh, Pf, Rh = 0.0, 1.0, r_assumed
    window, nis = [], []
    for k in range(N):
        Pp = Pf + q_true                 # predict (F = 1)
        S = Pp + Rh                      # the promise
        nu = zs[k] - xh                  # the surprise
        nis.append(nu*nu/S)              # running self-audit (sf-21)
        window.append(nu*nu)
        if len(window) > W:
            window.pop(0)
        if adapt and len(window) == W:
            Chat = np.mean(window)       # measured surprise level
            Rh = max(1e-3, Chat - Pp)    # solve the promise for R
        K = Pp/(Pp + Rh)                 # gain under current belief
        xh = xh + K*nu
        Pf = (1-K)*Pp                    # (short form is fine here:
    return float(np.mean(nis[500:])), Rh #  scalar, exact-arithmetic K)

nis_frozen, _ = run(adapt=False)
nis_adapt, R_final = run(adapt=True)
print('frozen  (R assumed 1, true 4): mean NIS =', round(nis_frozen, 2))
print('adaptive (window W=50):        mean NIS =', round(nis_adapt, 2))
print('final R_hat =', round(R_final, 2), ' (true 4.0)')
# frozen mean NIS ~3.65 (the broken promise, cf. chapter 2's ratio);
# adaptive ~1.03 with R_hat ~3.85: consistency restored from innovations
# alone -- no ground truth was used anywhere above.

Chapter 4: Covariance Matching & Sage-Husa — Recursive Noise Estimation

A sliding window stores W innovations and recomputes a covariance every tick. An engineer with a microcontroller asks the obvious question: can I get the same adaptation in one line, no buffer? Sage and Husa answered in 1969 — with a formula and, hidden inside it, a trap.

Derive the recursive form from the window. A sample mean over a growing window obeys m_k = m_k−1 + (1/k)(new − m_k−1). Replace 1/k with a forgetting schedule d_k = (1−b)/(1−b^k+1) built from a forgetting factor b (typically 0.95–0.999): early on, d_k is near 1/k (use all data); asymptotically it settles at 1−b (an exponential window of effective length 1/(1−b)). Apply this averaging to chapter 3's IAE quantity and you get the Sage-Husa measurement-noise estimator:

Worked Numbers — One Update by Hand, Then a Seeded Run

The hand step. Forgetting factor b = 0.95, so the asymptotic weight is d = 0.05. Previous R̂ = 2.5, innovation ν = 3.1, H P_pred H^T = 0.6. Instantaneous evidence: ν² − 0.6 = 9.61 − 0.6 = 9.01 — a one-sample estimate of R, huge because this particular innovation was large. Blend with the memory: R̂ = 0.95 × 2.5 + 0.05 × 9.01 = 2.375 + 0.4505 = 2.8255. The estimate moves 5% of the way toward the instantaneous evidence — that 5% IS the forgetting factor's asymptotic weight d = 1−b. One multiply-accumulate per tick; the buffer is gone.

The seeded run. The chapter-3 system (q = 0.09 known, true r = 4.0, filter starting at R̂ = 1.0), b = 0.995, 5000 steps, seed 7, R̂ floored at 10⁻⁶. The early schedule d_k ≈ 1/k learns fast from nothing: R̂ climbs to 3.54 by step 500 and 3.99 by step 2000 — converged within noise. But it does not STAY there: over steps 3000–5000 the estimate wanders with mean 3.78, because an exponential window of effective length 1/(1−b) = 200 has irreducible relative jitter of roughly √(2/200) ≈ 10%. The wander is not a bug to fix but a budget to spend: shrink it only by growing the memory, which slows real-change tracking by the same factor.

Be honest about the estimator's texture: each instantaneous term ν² − H P_pred H^T is a ONE-SAMPLE variance estimate — brutally noisy, and frequently negative (whenever a small innovation arrives, ν² < H P_pred H^T). The exponential average smooths it, but two engineering guards are mandatory: floor R̂ at a small positive value (a negative or zero R breaks the gain arithmetic and the gate), and treat persistent flooring as an alarm, not a state.

The forgetting factor b controls the same tradeoff as chapter 3's W, because it is the same object: b → 1 gives long memory — low estimator variance, slow response to a genuinely degrading sensor; small b gives fast tracking that jitters. The wander scales like √(2(1−b)) relative — at b = 0.995, about 10% — and that jitter feeds straight into the gain. You are always choosing where to spend the same statistical budget: confidence in the noise estimate versus freshness of it.

Now the trap the title promised. Sage and Husa also published the dual Q estimator — Q̂ from the state-correction statistics (Kνν^TK^T plus the covariance-decrease terms). Run BOTH estimators simultaneously and the scheme can diverge, and quals examiners love asking why: chapter 1 proved the steady-state innovation statistics constrain roughly the SUM H P_pred(q) H^T + R — one observable level, two unknowns. The joint estimator has a ridge of near-fixed-points: raise Q̂, lower R̂ compatibly, and the innovations look statistically identical. Estimation noise performs a random walk along the ridge until the gain is nonsense. The whiteness statistic breaks the tie in principle — Q-mistuning colors innovations, pure R-scale mistuning does not — but the classic recursions never use it.

Know the map beyond the classics. Offline, autocovariance least squares (ALS, Odelson-Rawlings 2006) solves for Q and R jointly from logged innovation autocorrelations at MULTIPLE lags — the extra lags restore the identifiability that the single-lag level lacks. Online, variational Bayes adaptive filters put an inverse-Wishart posterior on R and update it jointly with the state (Sarkka-Nummenmaa 2009) — principled uncertainty on the noise itself. Sage-Husa survives in practice because it costs nothing; know what it approximates, and where its blind spot is.

Chapter 5: Fading Memory — Forgetting as a Feature

Every divergence story in this lesson ends the same way: the filter becomes too sure of a past that no longer describes the present. The bluntest cure is also the most elegant: make the filter constitutionally incapable of remembering too much.

The mechanism is one line. In the predict step, replace P_pred = F P F^T + Q with

with a fading factor α (Fagin 1964, Sorenson-Sacks 1971). Inflating covariance is the filter's language for humility: each tick, the filter is forced to concede a fixed fraction of its accumulated certainty.

"Inflate P" sounds like a hack until you see what it means in data terms. Unroll the recursion: a measurement j steps old enters the current estimate with weight proportional to α^−2j. The fading filter is EXACTLY exponentially-weighted least squares. The effective data window is N_eff = α²/(α² − 1): α = 1.01 remembers about 51 samples, α = 1.05 about 11. One knob, direct control over memory length.

Worked Numbers — the Gain-Floor Theorem

Take the WORST case for responsiveness: q = 0. The model claims the state never changes, so the standard filter's gain decays to zero and it falls asleep permanently. With fading (measurement variance r):

Steady state: P_pred = α²(1 − K) P_pred with K = P_pred/(P_pred + r). Substituting K and solving: P_pred + r = α² r, so P_pred = (α² − 1) r — the predicted variance refuses to fall below a floor set purely by α.

The gain floor follows: K_∞ = P_pred/(P_pred + r) = (α² − 1)/α². For α = 1.05: α² = 1.1025, K_∞ = 0.1025/1.1025 = 0.092971. For α = 1.01: K_∞ = 0.019704.

Effective memory: N_eff = α²/(α² − 1): α = 1.05 gives 10.76 samples; α = 1.01 gives 50.75. The two views agree: K_∞ ≈ 1/N_eff — a gain floor IS a memory ceiling.

Verification by iteration: starting P = 5 with α = 1.05, q = 0, the faded Riccati recursion settles at K = 0.092971 and P_pred = 0.1025 = (α² − 1) r exactly — the index-card theorem survives contact with the recursion.

The insurance premium: the un-faded filter on a truly constant state would drive variance to 0; the faded one pays a permanent steady-state P_pred of 0.1025. That is the fixed cost of never being fully asleep.

Place the tool against chapters 3–4. Adaptive estimation MEASURES the mismatch and corrects the model — accurate, but with a lag, an estimator variance, and failure modes (ridges, gated windows). Fading memory measures nothing. It is open-loop humility: always on, zero parameters estimated, working against ANY slow mismatch — drifting biases, aging sensors, unmodeled dynamics — including ones your adaptive estimator was not designed to see. The price is permanent: steady-state variance is inflated even when the model happens to be right. An insurance premium, paid every tick.

You will meet this idea wearing other costumes. Fading with factor α² is the same medicine as adding pseudo-noise Q_extra proportional to F P F^T (state-proportional inflation), and the same idea as the forgetting factor in recursive least squares (λ = α⁻² — EE269's adaptive filters are this lesson's cousins). In the EKF it earns its keep twice over: linearization error grows with P, and fading keeps the filter responsive enough to re-linearize before the error compounds — many production EKFs run α = 1.01–1.05 as a matter of policy (see ekf.html).

The honest tuning doctrine: choose α from the timescale of trust — if you believe your model for roughly N steps, set α² = N/(N−1). Watch the failure directions: α too large discards information you paid for (the estimate jitters at the gain floor), and fading is NOT a substitute for gross mistuning — it bounds the sleep depth but cannot fix a wrong H or a biased sensor. It composes cleanly with chapter 3: adapt R with IAE, fade against everything else — belt and suspenders, and the showcase runs exactly that stack.

Chapter 6: A Bank of Filters — MMAE and the Road to IMM

Chapters 3–5 assumed the model family was right and only the noise levels lied. But sometimes the disagreement is structural: is the target cruising or maneuvering? Is the sensor healthy or degraded? You cannot slide continuously between hypotheses. So stop sliding — run them all.

Reframe adaptation as hypothesis testing. Postulate a finite set of models M₁..M_N — here, the same F and H with different Q hypotheses (q, 10q, 100q), though the same machinery covers different R, F, or H. Each model gets its own complete Kalman filter running on the same measurements. This is the multiple-model adaptive estimator (MMAE, Magill 1965) — the Bayesian answer to "which confession is true?"

The referee takes two lines to derive. The filter matched to the true model produces innovations that are N(0, S_i) — that is chapter 2's promise. So the likelihood of the observed innovation under model i is:

Note carefully what is being scored: not smallness of the innovation, but AGREEMENT WITH THE PROMISE. A model that predicted big surprises and got them outscores a model that promised small surprises and got big ones. The exponent penalizes surprise relative to the promise; the 1/√(2πS) normalization taxes over-wide promises. The referee scores calibration — the same virtue NIS tests.

Worked Numbers — the Referee's Arithmetic

A two-filter bank with equal priors p = (0.5, 0.5). The quiet filter promises innovation variance S₁ = 1.05; the maneuver filter promises S₂ = 4.0. A maneuver begins, and innovations ν = 2.5, −2.2, 2.8 arrive.

First innovation ν = 2.5: L₁ = exp(−6.25/2.1)/√(2π × 1.05) = 0.01985; L₂ = exp(−6.25/8)/√(2π × 4) = 0.09132. The quiet model called a 2.4-σ event; the maneuver model called a comfortable 1.25-σ event. Bayes: unnormalized (0.5 × 0.01985, 0.5 × 0.09132), normalized: p = (0.1786, 0.8214). One measurement moved the bank from agnostic to 82% maneuver.

Second innovation ν = −2.2: the likelihood ratio again favors the wide promise; p = (0.0719, 0.9281).

Third innovation ν = 2.8: p = (0.0095, 0.9905). Three consistent surprises produce 99% confidence — multiplicative Bayes compounds evidence exponentially fast. No threshold was tuned, no detector designed; the promise-versus-evidence ratio did all the work.

Read the design lesson in both directions: the same compounding that detects the maneuver in 3 ticks would, after 1000 quiet ticks, leave p(maneuver) around 10⁻⁴⁰ if unfloored. Detection speed and saturation depth are the same mechanism — which is exactly what IMM's Markov mixing repairs.

The bank has two outputs, and you should know when to use which. The MMSE estimate is the probability-weighted blend Σ p_i x̂_i — smooth, hedges during ambiguity. The MAP choice takes the highest-probability filter's estimate — commits, better when models are genuinely exclusive. And the blended covariance must include the spread-of-means term, Σ p_i(P_i + (x̂_i − x̂)(x̂_i − x̂)^T), because disagreement between filters IS uncertainty — dropping that term reproduces this lesson's original sin: overconfidence.

Now the structural flaw that makes this chapter a bridge rather than a destination. STATIC MMAE assumes one model is true for all time. The probability update is multiplicative, so log-odds accumulate; after a long cruise, p(maneuver) decays exponentially toward zero — and when the target finally does turn, the bank must climb out of a log-probability pit that deepened for the entire cruise. The referee saturates. Band-aids exist — floor the probabilities, exponentially discount old likelihoods (chapter 5's forgetting trick applied to log-odds) — but they treat the symptom.

Name the disease precisely: "which model?" is not a static unknown. It is a STATE that switches — cruise, then turn, then cruise — with its own dynamics. The principled fix gives model probability its own transition model: a Markov chain over modes, mixing the bank's estimates before each cycle so no hypothesis ever starves. That is the Interacting Multiple Model (IMM) filter, and it is the entire subject of the sibling lesson imm-estimation.html: this chapter hands you the bank and the referee; IMM adds the mode dynamics, the mixing step, and the full worked cycle.

Chapter 7: Numerics I — Roundoff, the Joseph Form, and Sequential Updates

An adaptive filter changes its gain constantly, from noisy estimates, in finite-precision arithmetic. Time to confess something about the covariance update you have been using since kalman-filter.html: it silently assumes the gain is EXACTLY right — and yours never is.

Derive the general truth first. If you update with ANY gain K, the posterior error is e_post = (I − KH)e_pred − Kv, with e_pred and v independent. Its true covariance is therefore:

No optimality invoked anywhere: this formula is just bookkeeping for whatever estimator you actually built. Substitute the OPTIMAL K and the cross terms collapse algebraically to the short form P = (I − KH)P_pred. The short form is a special case that happens to be cheaper — and everyone shipped it.

Now expose the fragility with a derivative. Perturb K to K + δ. The Joseph form, being the true variance, changes only at second order — dP/dK vanishes at the optimum, because the optimal gain minimizes posterior variance; that is what optimal means. The short form is LINEAR in K: its error is first-order in δ, and worse, in the wrong direction — increase K past optimal and (I − KH)P_pred DECREASES, claiming improvement while the actual variance grows. The short form does not report your estimator's variance; it reports what the variance would have been had your K been perfect.

Worked Numbers — the Impossible Improvement and the Negative Variance

Scalar update, exact arithmetic: P_pred = 1, R = 1, so S = 2 and the optimal gain is K = 0.5 with true posterior variance 0.5.

At K = 0.5 both forms agree: short (1−0.5)×1 = 0.5; Joseph (1−0.5)²×1 + 0.5²×1 = 0.25 + 0.25 = 0.5. At the exact optimum the forms are algebraically identical.

At K = 0.6 (a 20% gain error — one bad adaptive step): the short form claims (1−0.6)×1 = 0.40 — an improvement over the optimal 0.5, which is impossible: no gain can beat the optimum. The short form is not measuring the estimator; it is extrapolating a formula outside its validity. Joseph: (0.4)²×1 + (0.6)²×1 = 0.16 + 0.36 = 0.52 — which IS the true variance of the K = 0.6 estimator, since Var[(1−K)e_pred + Kv] = (1−K)²P_pred + K²R. A 2-million-sample Monte Carlo agrees: 0.5196. Note 0.52 = 0.5 + 2×(0.1)² — second-order in the gain error, the optimality of K = 0.5 made visible.

At K = 1.2 (an adaptive overshoot): the short form gives (1−1.2)×1 = −0.20 — a negative variance, which downstream becomes a negative S, an imaginary gate threshold, and a sign-flipped gain. The filter is not degraded; it is insane. Joseph: (−0.2)² + (1.2)² = 0.04 + 1.44 = 1.48 — a terrible estimator, honestly reported.

Floating point makes "wrong gain" not exotic but guaranteed: even the optimal K is stored rounded, so the short form is ALWAYS evaluated off-optimum. Its first-order error accumulates over thousands of updates, and because (I − KH)P is a product of non-symmetric factors, P drifts asymmetric (P₀₁ ≠ P₁₀) and eventually indefinite. In float32 — still the currency of embedded autopilots and GPU-resident filter banks — a strong measurement (small R) can push the naive update's smallest eigenvalue negative within a handful of steps. Joseph's quadratic structure plus an explicit symmetrization P ← (P + P^T)/2 degrades gracefully under the same conditions. Cost accounting: Joseph is roughly 2x the flops of the short form — the cheapest robustness you will buy this lesson.

While we are at the update step, take the second numerics upgrade: sequential processing. When R is diagonal (independent sensor channels), a vector measurement need not be absorbed in one matrix-inverse gulp. Process each scalar component through its own scalar update — S_i = H_i P H_i^T + r_i, one division each, no matrix inversion at all. The posterior is exactly the same: conditioning on independent evidence commutes. Three wins: no inverse to go ill-conditioned; per-channel innovation gates (chapter 2's gate now isolates WHICH channel is lying — one bad GPS pseudorange gets rejected without discarding the epoch); and per-channel adaptive-R hooks for chapter 3's machinery.

From Scratch: One Update, Three Ways

A vector KF update written as short form, Joseph, and sequential scalar channels — agreeing at the optimum, then deliberately fed an overshooting adaptive gain. Watch the assertions: the short form's covariance goes indefinite (a negative eigenvalue); Joseph reports the honest, larger covariance of the gain actually used.

python
import numpy as np

rng = np.random.default_rng(3)
n = 2
P = np.array([[1.0, 0.3], [0.3, 0.5]])   # predicted covariance
H = np.array([[1.0, 0.0], [0.0, 1.0]])   # two independent channels
R = np.diag([1.0, 2.25])                  # diagonal -> sequential OK

def joseph(P, K, H, R):
    IKH = np.eye(n) - K @ H
    Pn = IKH @ P @ IKH.T + K @ R @ K.T
    return 0.5*(Pn + Pn.T)                # explicit symmetrization

# --- optimal gain: all forms agree ---
S = H @ P @ H.T + R
K = P @ H.T @ np.linalg.inv(S)
P_short = (np.eye(n) - K @ H) @ P
P_jos = joseph(P, K, H, R)
assert np.allclose(P_short, P_jos), 'forms agree at the optimum'

# --- sequential scalar processing: same posterior, no inversion ---
Ps = P.copy()
for i in range(2):                        # one channel at a time
    h = H[i:i+1, :]
    Si = float(h @ Ps @ h.T) + R[i, i]    # scalar S: one division
    Ki = (Ps @ h.T) / Si
    IKHi = np.eye(n) - Ki @ h
    Ps = IKHi @ Ps @ IKHi.T + (Ki * R[i, i]) @ Ki.T   # Joseph per channel
    Ps = 0.5*(Ps + Ps.T)
assert np.allclose(Ps, P_jos, atol=1e-12), 'sequential == batch'

# --- perturbed gain (one bad adaptive transient): forms diverge ---
Kbad = K * 2.2                            # gain overshoot, 2.2x optimal
P_short_bad = (np.eye(n) - Kbad @ H) @ P
P_jos_bad = joseph(P, Kbad, H, R)
print('short-form eigenvalues :', np.round(np.linalg.eigvalsh(P_short_bad), 4))
print('joseph eigenvalues     :', np.round(np.linalg.eigvalsh(P_jos_bad), 4))
assert np.linalg.eigvalsh(P_short_bad).min() < 0   # short form: indefinite!
assert np.linalg.eigvalsh(P_jos_bad).min() > 0     # joseph: honest and PD
# prints short-form eigenvalues [-0.138, 0.2524]: the short form under an
# overshooting gain reports a covariance with a NEGATIVE eigenvalue -- a
# filter that believes an impossible thing. Joseph reports the honest
# (larger) covariance of the gain actually used, eigenvalues
# [0.3841, 1.393]: that honesty is what an adaptive filter, whose gain is
# always slightly wrong, is buying for its ~2x flops.

Set the boundary honestly, to hand off to chapter 8: Joseph plus symmetrization plus sequential updates fixes the UPDATE STEP for realistic conditioning, and it is the right default for every filter in this lesson. But when the condition number of P itself approaches the reciprocal of machine precision — strong measurements, near-perfect correlations from an adaptive transient, float32 budgets — no arrangement of P-arithmetic survives, because P can no longer even REPRESENT its own small eigenvalues. Then you stop computing with P, and start computing with its square root.

Chapter 8: Numerics II — Square-Root and UD Filtering

Chapter 7 patched the update formula. But there is a deeper problem no formula fixes: at extreme conditioning, the matrix P cannot faithfully HOLD its own small eigenvalues in finite precision. The 1969 Apollo-era answer remains the state of the art: never store P at all.

Diagnose the representation problem, not the algorithm. P's eigenvalues span its condition number κ(P); float32 carries about 7.2 decimal digits; and the entries of P mix the large and small eigendirections, so the small eigenvalue lives in the trailing digits of every entry. Quantify with this chapter's running example: P = [[1, 0.999],[0.999, 1]] has eigenvalues 1.999 and 0.001 — κ = 1999, so log₁₀(κ) = 3.3 of your 7.2 digits are spent before any algorithm runs. Perturb a single entry by 0.1% (0.999 → 1.000) and the small eigenvalue is annihilated: the matrix goes exactly singular. The weak direction was never safely stored.

The cure is a factorization with a theorem attached. Store S with P = S S^T — any square root; Cholesky's lower-triangular factor is the canonical choice. Since P's eigenvalues are the squares of S's singular values:

The log-condition halves. For the running example: 44.71 versus 1999. In digit terms: 1.65 digits spent instead of 3.3. The rule of thumb this buys is famous and worth stating exactly: square-root filtering achieves the same effective precision as covariance filtering with TWICE the word length — float32 square-root roughly matches float64 naive.

Worked Numbers — Conditioning, then a Float32 Stress Test

(a) Conditioning. Eigenvalues of P: 1.999 and 0.001, so κ(P) = 1999.0. Cholesky: S = [[1, 0],[0.999, 0.04471]] — note the weak direction gets its OWN entry, 0.04471 = √(1 − 0.999²). κ(S) = 44.71 = √1999: the doubling theorem, numerically. Now perturb P's 0.999 to 1.000 — a 0.1% change to one entry: the eigenvalues become (0, 2). The small eigenvalue is 100% destroyed; P is singular. The same-sized relative perturbation applied to S's off-diagonal entry moves the reconstructed small eigenvalue by 0.1% — three orders of magnitude less damage from the same-sized wound. Digit accounting: at κ(P) = 10⁷, P in float32 has NO digits left for the weak direction while S retains ~3.7 — which predicts exactly what (b) shows.

(b) Precision — the Grewal-Andrews stress test. From P₀ = I, absorb two near-collinear strong measurements (H = [1,1] then [1, 1+ε], R = ε², ε = 2.5×10⁻⁴). Naive float32 covariance update: the smallest eigenvalue of P comes out exactly 0.000e+00 — float32 could not represent the surviving sliver of uncertainty, so the filter now claims PERFECT knowledge along one direction, and the next measurement in that direction divides by a vanishing S. The exact float64 answer is 1.562e-08. A Potter square-root filter in the SAME float32: smallest eigenvalue 2.668e-08 — not exact (the problem sits at float32's edge) but positive and the right order of magnitude. Same arithmetic budget, qualitatively different integrity.

Quals expect the names and the deltas of the classical family. Potter (1963, flew on Apollo): measurement update directly on S for scalar measurements, one square root per update. Carlson: triangular-factor variant. Bierman UD (1977): factor P = U D U^T with U unit-upper-triangular and D diagonal; the measurement update is a rank-one modification needing NO square roots at all — pure multiply-add-divide, the reason it owned 1970s–80s flight computers and still owns radiation-hardened and fixed-point targets. Thornton: the companion time update, propagating U and D through F and Q via modified weighted Gram-Schmidt. UD is algebraically a square-root method (S = U √D) with the numerical benefits and integer-friendly arithmetic.

The modern framing connects to least squares. Today's preferred formulation is the array (QR) form: stack the propagated factor and noise factors into a tall block matrix, orthogonally triangularize, read off the updated factor. Numerically this is solving the underlying least-squares problem the way numerical analysts insist all least squares be solved — by orthogonal transforms on the square root of the information, never by forming the normal equations. P-arithmetic IS the normal equations, squared condition number and all. This is the same principle that makes factor-graph smoothers do QR on square-root information matrices — the thread runs from Apollo to modern SLAM unchanged.

Production Framing: Bierman's UD Update, As Flight Software Ships It

python
import numpy as np

def ud_factorize(P):
    """P = U D U' with U unit upper-triangular, D diagonal (one-time)."""
    n = P.shape[0]
    U, D, P = np.eye(n), np.zeros(n), P.copy()
    for j in range(n-1, -1, -1):
        D[j] = P[j, j]
        for i in range(j):
            U[i, j] = P[i, j] / D[j]
            for k in range(i+1):
                P[k, i] -= U[k, j] * D[j] * U[i, j]
    return U, D

def bierman_update(x, U, D, h, R, z):
    """Scalar-measurement UD update (Bierman 1977). No square roots:
    multiply/add/divide only -- fixed-point and certification friendly.
    D stays positive BY CONSTRUCTION: indefiniteness is unrepresentable."""
    f = U.T @ h                # f = U'h
    v = D * f                  # v_i = D_i f_i
    alpha = R                  # innovation variance accumulator
    b = np.zeros(len(x))
    for j in range(len(x)):
        alpha_prev = alpha
        alpha += f[j]*v[j]
        D[j] *= alpha_prev/alpha
        p = -f[j]/alpha_prev
        bj = v[j]
        for i in range(j):
            Uij = U[i, j]
            U[i, j] = Uij + b[i]*p
            b[i] += Uij*bj
        b[j] = bj
    K = b/alpha                # gain, free byproduct
    return x + K*(z - h @ x), U, D

# verify against the textbook update
rng = np.random.default_rng(3)
A = rng.normal(size=(3, 3)); P = A @ A.T + 3*np.eye(3)
x = rng.normal(size=3); h = rng.normal(size=3); R, z = 0.7, 1.9
U, D = ud_factorize(P)
assert np.allclose(U @ np.diag(D) @ U.T, P)
x2, U2, D2 = bierman_update(x.copy(), U.copy(), D.copy(), h, R, z)
S = h @ P @ h + R; K = P @ h / S
x_ref = x + K*(z - h @ x)
P_ref = P - np.outer(K, h @ P)
assert np.allclose(x2, x_ref)
assert np.allclose(U2 @ np.diag(D2) @ U2.T, P_ref)
assert D2.min() > 0
print('UD update == reference KF: PASS   min(D) =', round(D2.min(), 4))
# Production notes: pair with Thornton's MWGS time update to propagate
# U, D through F and Q; process vector measurements as sequential scalar
# channels (decorrelate R first if needed). Modern alternative: the
# array/QR square-root form -- orthogonal triangularization of stacked
# factors -- the same numerics that power square-root information
# smoothers in factor-graph SLAM.

Close with the engineering decision table and the adaptive tie-in. Modern float64 plus Joseph plus symmetrization is adequate for most well-conditioned filters — square-root costs roughly 1.5–2x and is the wrong complexity to buy reflexively. Buy it when: float32 or fixed-point hardware (embedded, GPU-resident filter banks); condition numbers pushed by strong-sensor/weak-sensor mixes or by adaptive transients (an adaptive filter slamming R̂ down 4x in one window is a conditioning event BY DESIGN); long-mission covariance integrity requirements; or certification demands provable PSD. And note the adaptive synergy runs both ways: chapter 3's machinery keeps S honest statistically; this chapter keeps it honest arithmetically — the showcase runs the safe update forms underneath the adaptive loop.

Chapter 9: Showcase — the Self-Tuning Tracker

You have the instruments (chapter 2), the estimators (chapters 3–4), the insurance (chapter 5), the referee (chapter 6), and trustworthy arithmetic (chapters 7–8). Time to fly the whole stack — and to break it on purpose, from the cockpit.

The cockpit: a 2D arena where a target flies constant-velocity legs punctuated by coordinated turns — the classic maneuvering profile. A position sensor reports every tick with true measurement variance r_true. The filter is a 4-state CV Kalman tracker — exactly kalman-filter.html's chapter-8 machine — but its ASSUMED q and r are set by two mischief sliders (0.01x to 100x truth). The adaptive stack sits behind one master toggle: innovation-window R estimation (chapter 3), a fading factor on the predict (chapter 5), Joseph-form sequential updates underneath (chapter 7).

One honest detail before you fly. The 1.0x baseline q is not the raw random-acceleration variance of the truth — it is chapter 1's "deliberate crime" performed correctly: q inflated about 8x so the coordinated turns live inside the process-noise budget. That is the tuning a deployed CV tracker actually ships with, and it is why 1.0x/1.0x reads consistent: YES on a maneuvering target — with only brief instrument blips as each turn passes through the windows. Everything the mischief sliders do is measured against that covering tuning.

Experiment 1 — recreate chapter 0, honestly. Drag Q to 0.01x, adaptation OFF. During straight legs the picture is seductive: glassy estimate, tiny ellipse. Watch the instruments instead: mean NIS creeps above the 2-dof 95% band [0.05, 7.38] (sf-21's chart, running live), and the lag-1 autocorrelation bar pushes outside its whiteness band. Then the first turn: the estimate cuts the corner wide, NIS spikes past 20, and — the detail worth staring at — the ellipse does not grow, because the filter's Q says turns do not happen. Confidence and error have fully decoupled.

Experiment 2 — the opposite quadrant. Q back to truth, R to 0.01x. The estimate shivers, gluing itself to every measurement. Now read the instrument correctly: with R assumed tiny, the promised S is small, so NIS = ν²/S runs HIGH — the filter is more surprised than its promise allows. Contrast with R assumed huge (100x): S is fat, NIS sinks BELOW the band, and the filter free-runs. The panel teaches the diagnosis direction: NIS high means promised S too small (overconfident somewhere); NIS low means S too fat (underconfident); the autocorrelation bars say whether the GAIN — the q/r balance — is also wrong.

Experiment 3 — the payoff. Leave the sliders wrong. Flip ADAPT on. The R̂ readout walks from its wrong start toward r_true within roughly two window-lengths; the fading factor keeps the gain floor up so maneuvers stop being catastrophic; mean NIS bends back into the band; the whiteness bars sink into their thresholds. Then fight it: yank the sliders mid-flight, hit Storm (r_true doubles), trigger extra turns. The stack re-converges each time, with the lag the window slider predicts.

Experiment 4 — find the limit. Set BOTH Q and R wrong by the same factor (both 10x). The estimates are fine — chapter 1's ratio invariance; the gain was never wrong. NIS is broken by the common scale, and adaptation fixes NIS — but look at the R̂ readout against truth: the stack lands on A consistent scale, not necessarily THE true split. The identifiability ridge from chapter 4, now experienced: consistency is restorable from innovations alone; the true decomposition is not, without more structure. The display says so explicitly — it reports "consistent" and "correct" separately — rather than pretending.

Chapter 10: Quals Arsenal

The machinery is built. Now sit the exam: seven questions of the kind that open an estimation quals, each expecting derivation plus limits plus judgment — and the numbers you verified in this lesson are your answer key.

Protocol: attempt each question on paper or the Teach whiteboard BEFORE expanding the worked answer. The graded skill is not recall but the chain derivation → limit → decision — every answer here ends with a judgment call, because every examiner's follow-up is "so what would you ship?"

A map of the seven, so a miss tells you which chapter to revisit: gain-and-identifiability theory (Q1, Q3), the adaptive estimators and their fixed points (Q2, Q4), robustness mechanisms (Q5 fading, Q6 MMAE), and numerics (Q7 Joseph/square-root).

Debates — Oral-Exam Finishers

Defensible positions exist on both sides of each. But they must be argued with the detection-lag / estimator-variance / identifiability vocabulary this lesson built — "it felt more robust" is a failing answer.

Chapter 11: Connections — the Filter That Doubts Itself

Strip the algebra and this was a lesson about one signal: the gap between how surprised the filter said it would be and how surprised it was. The innovation promise S = H P_pred H^T + R is the filter's testable claim about itself. IAE and Sage-Husa CLOSE THE LOOP on it; fading memory BOUNDS the damage when the loop is too slow; MMAE runs the loop over DISCRETE hypotheses; and Joseph/square-root arithmetic guarantees the claim is at least computed honestly. One signal, five disciplines.

Method Selection

Habitats

GNSS/INS integration adapts R against multipath and jamming while fading covers IMU aging — sf-14-ins-gnss-coupling.html is the natural habitat. EKF-based VIO inflates Q for linearization error and runs Joseph religiously (ekf.html). Target trackers hand the maneuver problem to full IMM — the sibling lesson imm-estimation.html — with this lesson's MMAE as the vestibule. Flight software still ships UD filtering, and GPU-resident filter banks are re-discovering float32 square-root discipline (chapter 8).

What This Lesson Deliberately Did Not Cover

The research frontier argues in this lesson's vocabulary: learned noise models trained through differentiable filters (sf-20); adaptive filtering under adversarial spoofing, where the innovation stream itself is attacked; and interacting multiple-model banks over learned dynamics. Promises, identifiability, gain floors — the words do not change.