Incorporated relaxation into CovarianceMatching

examples/equipped.rs

@@ -183,7 +183,8 @@ fn main() -> Result<(), FilterError> {
         // Reject a residual too large to be this filter's own noise. 0.99 throws
         // away one good measurement in a hundred to do it.
         .with_gate(ResidualGate::new(0.99)?)
-        // Let the innovations set the size of Q. One free parameter, the window.
+        // Let the innovations set the size of Q. Two tuning parameters: the window
+        // length here, and a relaxation factor left at its default.
         .with_adaptation(CovarianceMatching::inflating(20)?)
         // And discount old measurements slightly, so the covariance keeps a floor
         // under it and the filter never stops listening.
@@ -275,14 +276,20 @@ fn main() -> Result<(), FilterError> {
     // hundred, and 120 steps bought about one. Tightening the gate to 0.999 would
     // buy that measurement back and start letting smaller outliers through.
     //
-    // **The Q scale does not glide up at the maneuver; it thrashes.** It is at the
-    // clamp of 1.0 most of the time and jumps into the hundreds for a few steps at
-    // a stretch. That is this mechanism's real behavior, and it is a feedback loop
-    // rather than noise: inflating the scale makes the next prediction's covariance
-    // larger, which makes the innovations look explicable again, which collapses
-    // the scale back to the clamp. A longer window damps the amplitude without
-    // removing the oscillation. It still works — see the RMS — but anyone plotting
-    // the scale expecting a tracking signal will be surprised.
+    // **The Q scale is not a smooth tracking signal**, even damped. It spends
+    // stretches at its floor of 1.0 and jumps into the tens or hundreds. Two
+    // separate reasons, and neither is a defect: with a finite window the solved
+    // covariance comes out negative routinely, so the floor is reached often
+    // rather than rarely; and inflating the scale enlarges the next prediction's
+    // covariance, which makes the innovations look explicable again, which pulls
+    // the scale back down.
+    //
+    // The relaxation factor damps that, and the default is damped for a reason —
+    // undamped, this mechanism leaves the filter overconfident by nearly a factor
+    // of two. What it does not do is give the scale a meaningful absolute level:
+    // the baseline depends on how visible the nominal Q is in the measurement and
+    // is not 1. Do not threshold it to detect faults; that is what NIS is for,
+    // and `diagnostics.rs` shows how.
     //
     // **The plain filter's final 1 sigma is the number to be suspicious of.** It is
     // *smaller* than the equipped filter's while its RMS is nearly three times

readme.md

@@ -79,9 +79,11 @@ let sensor = Position { deviation: 1.0 };
 
 // A gated filter with adaptive process noise, started from a loose prior.
 let prior = Gaussian::from_deviations(Vector2::new(0.0, 0.0), Vector2::new(10.0, 5.0));
+// The window and the relaxation are both tuning parameters; `inflating` picks a
+// default relaxation that is a starting point, not an answer.
 let mut filter = Ukf::new(prior)
     .with_gate(ResidualGate::new(0.99)?)
-    .with_adaptation(CovarianceMatching::inflating(20)?);
+    .with_adaptation(CovarianceMatching::inflating(20)?.with_relaxation(0.1)?);
 
 for step in 0..100 {
     let time = step as f64 * 0.5;
@@ -174,7 +176,7 @@ let delay = fired.map(|index| index.saturating_sub(onset));
   `[1.9533, 2.0473]`, and a 95% coverage rate of 0.9536.
 - **Each variant against its fault.** Gating against injected 30-sigma outliers:
   position RMS 4.035 to 1.602. Adaptive noise against an unmodeled 0.3 m/s^2:
-  mean NEES 30.7 to 9.5. Consider covariance against an unestimated scale factor:
+  mean NEES 30.7 to 4.2. Consider covariance against an unestimated scale factor:
   mean NEES 650 to 1.67, with the accuracy unchanged, which is the point. Bias
   augmentation against a 5-unit sensor offset: RMS 2.567 to 0.503 and the offset
   recovered as 4.985. An IMM across a mode switch: whole-arc RMS 0.770 against

src/adaptive.rs

@@ -10,9 +10,10 @@
 //!
 //! - [`CovarianceMatching`] scales a nominal `Q` by one number, chosen so that the
 //!   innovation covariance the filter predicts matches the one it is measuring.
-//!   It is derived rather than tuned, has a single free parameter (the window
-//!   length), and keeps whatever structure the nominal `Q` has. This is the one to
-//!   reach for first.
+//!   The scale is solved for rather than guessed, and it keeps whatever structure
+//!   the nominal `Q` has. This is the one to reach for first. **It has two tuning
+//!   parameters**, the window length and the relaxation factor, and both change
+//!   its behavior enough to need choosing deliberately.
 //! - [`InnovationNoise`] estimates `Q` outright from the state corrections the
 //!   filter has been making, after Mehra. It assumes no structure, which is both
 //!   its appeal and its weakness: the estimate has rank at most the measurement
@@ -21,14 +22,10 @@
 //! - [`NisInflation`] watches the normalized innovation squared against a
 //!   chi-squared threshold and inflates on exceedance, decaying back when things
 //!   are quiet. It reacts to a single step rather than to a window, which is the
-//!   reason to want it, at the cost of three parameters that interact. Note that
-//!   reacting per step is not the same as responding sooner: it multiplies by its
-//!   inflation factor per exceedance, so reaching a tenfold process noise at a
-//!   factor of two takes more than three consecutive exceedances, while the
-//!   covariance match solves for the scale and can get there in one update. On a
-//!   step change in the dynamics, measured over 120 runs, the match took a mean of
-//!   2.98 steps to reach tenfold and the inflater 7.34. See
-//!   `tests/variants.rs`.
+//!   reason to want it, at the cost of three parameters that interact. Its scale
+//!   has a meaningful baseline of exactly 1, which the covariance match's does
+//!   not — see the conditioning note on [`CovarianceMatching`] — so a threshold on
+//!   it means what it looks like.
 //!
 //! # What a mechanism gets to see
 //!
@@ -45,6 +42,22 @@ use serde::{Deserialize, Serialize};
 use crate::chi2;
 use crate::error::{FilterError, FilterResult};
 
+/// The default relaxation factor for [`CovarianceMatching`].
+///
+/// Not 1.0, which would be the plain solve. Solved from scratch the scale
+/// oscillates hard enough to leave the filter overconfident: measured over 80 runs
+/// of 120 steps with a step change in the dynamics, mean NEES came out at 3.69
+/// against an expectation of 2, where damping brought it to 1.77. The value where
+/// NEES crossed 2 on that problem was about 0.15.
+///
+/// **This is a tuning parameter with a default, not a constant of nature.** 0.1 is
+/// a reasonable starting point on the one problem it was measured on; the right
+/// value depends on how fast the process noise actually changes relative to the
+/// measurement rate, and it trades responsiveness against steadiness. Choose it
+/// deliberately for anything that matters, and
+/// [`CovarianceMatching::with_relaxation`] says how.
+pub const DEFAULT_RELAXATION: f64 = 0.1;
+
 /// What one update observed, reduced to what an adapter can use.
 ///
 /// The traces are of the measurement-space quantities projected through the
@@ -159,16 +172,57 @@ impl<const N: usize> AdaptiveProcessNoise<N> for FixedProcessNoise {
 ///
 /// # What it costs and what it assumes
 ///
-/// One free parameter beyond the clamps: the window length, which sets how fast
-/// the estimate responds against how noisy it is. The clamps are guards, not
-/// tuning: `minimum` keeps a quiet stretch from driving `Q` to zero and
-/// `maximum` keeps one outlier from opening the filter up.
+/// **Two tuning parameters**, and neither has an obviously right value:
+///
+/// - the **window length**, which trades how fast the estimate responds against
+///   how noisy it is
+/// - the **relaxation factor**, which trades the same thing again by a different
+///   route, and which defaults to [`DEFAULT_RELAXATION`] rather than to the plain
+///   solve for the reasons given there
+///
+/// They are not independent — both control the same trade — so tuning one while
+/// holding the other at a default is usually enough, and tuning both against a
+/// known fault will overfit to it.
+///
+/// The clamps are not tuning. `minimum` keeps a quiet stretch from driving `Q` to
+/// zero and `maximum` keeps one outlier from opening the filter up. Note that
+/// `minimum` is reached often rather than rarely: with a finite window the solved
+/// covariance comes out negative routinely, which is expected and is why clamping
+/// it is the conventional response.
 ///
 /// It assumes the nominal `Q` has the right *shape* and only the wrong size. Where
 /// the mismodeling is in a direction the nominal `Q` does not cover, the scale
 /// cannot find it; [`InnovationNoise`] can, at the cost of assuming nothing and
 /// therefore estimating much more from the same data.
 ///
+/// # It is ill-conditioned when the nominal `Q` is barely visible
+///
+/// `tr(H Q0 H')` is the denominator of the solve, so the scale's sensitivity to
+/// noise in the numerator goes as `tr(R) / tr(H Q0 H')`. A nominal `Q` that the
+/// measurement can hardly see makes the solved scale enormous and nearly
+/// meaningless — not wrong, but dominated by the sampling error of the numerator.
+///
+/// Measured on a constant-velocity filter against a one-meter position sensor,
+/// with nothing wrong at all, reporting the scale's mean level over 120 runs:
+///
+/// | acceleration noise intensity | mean scale with no fault |
+/// |---|---|
+/// | 0.01 | 228.3 |
+/// | 0.1 | 27.8 |
+/// | 1.0 | 4.84 |
+/// | 10.0 | 1.82 |
+///
+/// At the top of that table the mechanism reports a process noise two hundred
+/// times nominal on a filter whose model is exactly right. The scale is not
+/// reporting a fault; it is reporting that one unit of nominal `Q` is almost
+/// invisible against a measurement variance of 1, so the ratio is huge and noisy.
+///
+/// Two consequences. A threshold on the scale — "inflated means something is
+/// wrong" — is only meaningful once the baseline level is known for that `Q` and
+/// that sensor, and the baseline is not 1. And a very small nominal `Q` is better
+/// made larger than left to be rescued by adaptation, because the adaptation
+/// cannot resolve it.
+///
 /// # Feeding it more than one sensor
 ///
 /// Measurements of different dimensions can share one window. Each sample carries
@@ -188,6 +242,7 @@ pub struct CovarianceMatching {
     window: usize,
     minimum: f64,
     maximum: f64,
+    relaxation: f64,
     scale: f64,
     history: VecDeque<[f64; 4]>,
 }
@@ -226,6 +281,7 @@ impl CovarianceMatching {
             window,
             minimum,
             maximum,
+            relaxation: DEFAULT_RELAXATION,
             scale: 1.0,
             history: VecDeque::with_capacity(window),
         })
@@ -246,11 +302,94 @@ impl CovarianceMatching {
         Self::new(window, 1.0, 1e6)
     }
 
+    /// Sets the relaxation factor: how far toward the newly solved scale each
+    /// update moves.
+    ///
+    /// **A tuning parameter.** It defaults to [`DEFAULT_RELAXATION`], which is a
+    /// starting point measured on one problem and not a value to rely on
+    /// unexamined. Pass `1.0` to get the undamped solve.
+    ///
+    /// With `relaxation = r`, each update sets
+    /// `scale <- (1 - r) * scale + r * solved`, so `r = 1` replaces the scale
+    /// outright — the default, and the behavior of this type without this call —
+    /// and smaller values move it a fraction of the way.
+    ///
+    /// # Why this exists
+    ///
+    /// Solved from scratch, the scale is not a smooth tracking signal. Measured
+    /// over 60 runs of 120 steps with a step change in the dynamics: its standard
+    /// deviation is two to six times its mean at every window length, and it sits
+    /// at `minimum` between 62% and 74% of the time. Two things cause that, and
+    /// only one of them is noise.
+    ///
+    /// The first is expected and documented in the literature: with a finite
+    /// window the solved covariance frequently comes out negative, and clamping it
+    /// is the recommended response rather than a defect. The clamp is therefore
+    /// load-bearing most of the time, not a rarely-fired guard.
+    ///
+    /// The second is a feedback loop. Inflating the scale enlarges the next
+    /// prediction's covariance, which enlarges `tr(H P_sigma H')`, which shrinks
+    /// the numerator, which collapses the solved scale back to the clamp — and the
+    /// innovations then grow again. The mechanism is correcting against a quantity
+    /// it moved itself, one step earlier.
+    ///
+    /// # How to choose it
+    ///
+    /// The published trade is between smoothness and speed: a factor near zero
+    /// gives steady estimates and lags a real change in the dynamics, and one near
+    /// one responds immediately and keeps the oscillation. Measured on a step
+    /// change in the dynamics, 80 runs of 120 steps at a window of 20:
+    ///
+    /// | relaxation | position RMS | mean NEES |
+    /// |---|---|---|
+    /// | 1.00 | 0.804 | 3.69 |
+    /// | 0.25 | 0.777 | 2.28 |
+    /// | 0.10 | 0.754 | 1.77 |
+    /// | 0.05 | 0.745 | 1.62 |
+    ///
+    /// Accuracy barely moves; consistency moves a great deal, and NEES passing
+    /// through its expectation of 2 near `0.15` is the useful signal. That
+    /// suggests a way to pick it without tuning against the fault: run a nominal
+    /// ensemble and take the value at which the ensemble NEES matches the state
+    /// dimension, which [`crate::diagnostics::ensemble_test`] reports directly.
+    /// Choosing it by best RMS on the faulted case instead will overfit to that
+    /// fault.
+    ///
+    /// # Relation to the Sage-Husa filter
+    ///
+    /// This is the Sage-Husa adaptive filter's recursion,
+    /// `Q_k = (1 - d_k) Q_{k-1} + d_k Qhat_k`, applied to a scalar rather than to
+    /// the whole matrix. Sage-Husa ramps its weight as
+    /// `d_k = (1 - b) / (1 - b^(k+1))`, which converges to `1 - b`; this takes the
+    /// limit directly and holds it fixed, so early updates are weighted the same
+    /// as later ones. The ramp matters when the estimate starts from nothing, and
+    /// here it starts from a nominal `Q` a caller chose deliberately.
+    ///
+    /// # Errors
+    ///
+    /// [`FilterError::Configuration`] unless `relaxation` lies in `(0, 1]`.
+    pub fn with_relaxation(mut self, relaxation: f64) -> FilterResult<Self> {
+        if relaxation.is_nan() || relaxation <= 0.0 || relaxation > 1.0 {
+            return Err(FilterError::configuration(
+                "relaxation",
+                relaxation,
+                "inside (0, 1]",
+            ));
+        }
+        self.relaxation = relaxation;
+        Ok(self)
+    }
+
     /// The scale currently applied to the nominal process noise.
     pub fn scale(&self) -> f64 {
         self.scale
     }
 
+    /// The relaxation factor. `1.0` means the scale is replaced outright.
+    pub fn relaxation(&self) -> f64 {
+        self.relaxation
+    }
+
     /// How many updates the window holds.
     pub fn window(&self) -> usize {
         self.window
@@ -303,10 +442,15 @@ impl<const N: usize> AdaptiveProcessNoise<N> for CovarianceMatching {
         // All four are sums over the same window, so the ratio of sums is the
         // ratio of means and the window length cancels.
         let solved = (measured - noise - propagated) / template;
-        self.scale = solved.clamp(self.minimum, self.maximum);
-        if self.scale.is_nan() {
+        if solved.is_nan() {
             self.scale = self.minimum;
+            return;
         }
+        // Clamp before relaxing, so that a wild solved value cannot drag the
+        // carried scale outside the bounds on its way past.
+        let target = solved.clamp(self.minimum, self.maximum);
+        self.scale = ((1.0 - self.relaxation) * self.scale + self.relaxation * target)
+            .clamp(self.minimum, self.maximum);
     }
 
     fn process_noise(&self, nominal: &SMatrix<f64, N, N>) -> SMatrix<f64, N, N> {
@@ -427,18 +571,22 @@ impl<const N: usize> AdaptiveProcessNoise<N> for InnovationNoise<N> {
 ///
 /// # What it costs
 ///
-/// It reacts to a single bad step rather than waiting for a window. That is the
-/// reason to want it, but it does not follow that it responds sooner: the scale
-/// moves by one factor of `inflation` per exceedance, so reaching tenfold at
+/// It reacts to a single bad step rather than waiting for a window, which is the
+/// reason to want it. It does not follow that it responds sooner: the scale moves
+/// by one factor of `inflation` per exceedance, so reaching tenfold at
 /// `inflation = 2` takes better than three consecutive exceedances, and `decay`
-/// gives ground back on every step that does not exceed. Measured on a step change
-/// in the dynamics over 120 runs, [`CovarianceMatching`] reached a tenfold scale in
-/// a mean of 2.98 steps against this mechanism's 7.34.
+/// gives ground back on every step that does not exceed.
 ///
-/// It also has three parameters — threshold, inflation and decay — that interact,
-/// and all three set how large the response to a given fault is. Tuning them
-/// against a known fault and then reporting the filter as adaptive overstates what
-/// the mechanism does on a fault it was not tuned for.
+/// Three parameters — threshold, inflation and decay — that interact, and all
+/// three set how large the response to a given fault is. Tuning them against a
+/// known fault and then reporting the filter as adaptive overstates what the
+/// mechanism does on a fault it was not tuned for.
+///
+/// What it has over [`CovarianceMatching`] is a scale whose baseline is exactly
+/// one when nothing is wrong, whatever the nominal `Q`. The covariance match's
+/// baseline depends on how visible that `Q` is in the measurement and can sit in
+/// the hundreds with a perfectly correct model, so the two scales are not
+/// comparable quantities and only this one supports a fixed threshold.
 ///
 /// Note the interaction with a residual gate: a gate rejects exactly the
 /// measurements this mechanism would inflate on, and rejected measurements are not
@@ -549,10 +697,17 @@ mod tests {
     }
 
     /// Fed innovations consistent with a process noise four times the nominal, the
-    /// matcher returns a scale of four.
+    /// solve returns a scale of four.
+    ///
+    /// Undamped on purpose: this is a test of the moment-matching solve, so the
+    /// relaxation is turned off to look at it directly. What the default does to
+    /// the approach is `the_default_damps_the_approach` below.
     #[test]
     fn covariance_matching_recovers_a_known_scale() {
-        let mut matcher = CovarianceMatching::new(20, 1e-6, 1e6).unwrap();
+        let mut matcher = CovarianceMatching::new(20, 1e-6, 1e6)
+            .unwrap()
+            .with_relaxation(1.0)
+            .unwrap();
         // tr(H P H') = 2, tr(R) = 1, tr(H Q H') = 0.5. With Q scaled by 4 the
         // innovations carry tr(S) = 2 + 4*0.5 + 1 = 5.
         for _ in 0..20 {
@@ -590,7 +745,10 @@ mod tests {
         let a = sample(1, 5.0, 1.0, 2.0, 0.5);
         let b = sample(2, 71.0, 17.0, 30.0, 6.0);
 
-        let mut pooled = CovarianceMatching::new(40, 1e-9, 1e9).unwrap();
+        let mut pooled = CovarianceMatching::new(40, 1e-9, 1e9)
+            .unwrap()
+            .with_relaxation(1.0)
+            .unwrap();
         for _ in 0..20 {
             pooled.observe(&a);
             pooled.observe(&b);
@@ -601,7 +759,10 @@ mod tests {
         // template-trace-weighted mean of them, not an average of the two scales.
         // Here A implies 4 and B implies 1, with weights 0.5 and 6.0.
         let disagreeing = sample(2, 53.0, 17.0, 30.0, 6.0);
-        let mut mixed = CovarianceMatching::new(40, 1e-9, 1e9).unwrap();
+        let mut mixed = CovarianceMatching::new(40, 1e-9, 1e9)
+            .unwrap()
+            .with_relaxation(1.0)
+            .unwrap();
         for _ in 0..20 {
             mixed.observe(&a);
             mixed.observe(&disagreeing);
@@ -624,6 +785,119 @@ mod tests {
         assert_eq!(matcher.scale(), 1.0);
     }
 
+    /// The default is [`DEFAULT_RELAXATION`], and passing it explicitly is the same
+    /// thing.
+    ///
+    /// Asserted so that changing the constant is a deliberate act with a failing
+    /// test attached, rather than something that quietly moves every measured
+    /// number in this crate — which is what it did when it moved from 1.0 to 0.1.
+    #[test]
+    fn the_default_relaxation_is_the_documented_constant() {
+        let default = CovarianceMatching::new(10, 1e-6, 1e6).unwrap();
+        assert_eq!(default.relaxation(), DEFAULT_RELAXATION);
+        assert_eq!(DEFAULT_RELAXATION, 0.1);
+        assert_eq!(
+            CovarianceMatching::inflating(10).unwrap().relaxation(),
+            DEFAULT_RELAXATION
+        );
+
+        let mut implicit = CovarianceMatching::new(10, 1e-6, 1e6).unwrap();
+        let mut explicit = CovarianceMatching::new(10, 1e-6, 1e6)
+            .unwrap()
+            .with_relaxation(DEFAULT_RELAXATION)
+            .unwrap();
+        for step in 0..30 {
+            // A varying sequence, so a difference in handling would show.
+            let squared = 5.0 + 3.0 * (step as f64 * 0.7).sin();
+            let sample = record(squared, 2.0, 0.5);
+            implicit.observe(&sample);
+            explicit.observe(&sample);
+            assert_eq!(implicit.scale(), explicit.scale(), "parted at step {step}");
+        }
+    }
+
+    /// The default reaches the same answer as the undamped solve, more slowly.
+    #[test]
+    fn the_default_damps_the_approach() {
+        // Innovations consistent with a scale of 4.
+        let sample = record(5.0, 2.0, 0.5);
+
+        let mut undamped = CovarianceMatching::new(1, 1e-6, 1e6)
+            .unwrap()
+            .with_relaxation(1.0)
+            .unwrap();
+        let mut damped = CovarianceMatching::new(1, 1e-6, 1e6).unwrap();
+
+        undamped.observe(&sample);
+        damped.observe(&sample);
+        // One update takes the undamped solve all the way and the default a tenth
+        // of the way: 1 + 0.1 * (4 - 1) = 1.3.
+        assert_relative_eq!(undamped.scale(), 4.0, max_relative = 1e-12);
+        assert_relative_eq!(damped.scale(), 1.3, max_relative = 1e-12);
+
+        // Given enough updates it gets there too. The gap closes by 0.9 each time,
+        // so after n it is 3 * 0.9^n; 300 puts that below 1e-13.
+        for _ in 0..300 {
+            damped.observe(&sample);
+        }
+        assert!(3.0 * 0.9f64.powi(300) < 1e-13);
+        assert_relative_eq!(damped.scale(), 4.0, max_relative = 1e-12);
+    }
+
+    /// Relaxation approaches the same answer, geometrically rather than at once.
+    #[test]
+    fn relaxation_damps_the_approach_without_moving_the_answer() {
+        // Innovations consistent with a scale of 4, from a starting scale of 1.
+        let sample = record(5.0, 2.0, 0.5);
+
+        let mut sharp = CovarianceMatching::new(1, 1e-6, 1e6)
+            .unwrap()
+            .with_relaxation(1.0)
+            .unwrap();
+        sharp.observe(&sample);
+        assert_relative_eq!(sharp.scale(), 4.0, max_relative = 1e-12);
+
+        let mut relaxed = CovarianceMatching::new(1, 1e-6, 1e6)
+            .unwrap()
+            .with_relaxation(0.25)
+            .unwrap();
+
+        // One step moves a quarter of the way: 1 + 0.25 * (4 - 1) = 1.75.
+        relaxed.observe(&sample);
+        assert_relative_eq!(relaxed.scale(), 1.75, max_relative = 1e-12);
+
+        // And it keeps closing the remaining gap by a quarter each time, so after
+        // `n` further updates the gap is `3 * 0.75^n`. That is the tolerance: at
+        // 40 updates it is 2e-5, which is why this runs to 120 and asserts 1e-12
+        // rather than asserting a round number at an arbitrary step count.
+        const FURTHER: i32 = 120;
+        let mut previous = relaxed.scale();
+        for _ in 0..FURTHER {
+            relaxed.observe(&sample);
+            assert!(relaxed.scale() > previous, "should still be climbing");
+            previous = relaxed.scale();
+        }
+        let remaining_gap = 3.0 * 0.75f64.powi(FURTHER) / 4.0;
+        assert!(remaining_gap < 1e-12);
+        assert_relative_eq!(relaxed.scale(), 4.0, max_relative = 1e-12);
+    }
+
+    /// Relaxing toward a target that is itself outside the clamps cannot carry the
+    /// scale outside them.
+    #[test]
+    fn relaxation_respects_the_clamps() {
+        let mut relaxed = CovarianceMatching::inflating(1)
+            .unwrap()
+            .with_relaxation(0.5)
+            .unwrap();
+        // Innovations smaller than predicted solve for a negative scale, which
+        // clamps to the floor of 1 before the relaxation sees it.
+        for _ in 0..20 {
+            relaxed.observe(&record(0.1, 2.0, 0.5));
+        }
+        assert_eq!(relaxed.scale(), 1.0);
+    }
+
     #[test]
     fn covariance_matching_ignores_an_invisible_nominal() {
         let mut matcher = CovarianceMatching::new(5, 1e-6, 1e6).unwrap();
@@ -681,6 +955,13 @@ mod tests {
         assert!(CovarianceMatching::new(0, 1.0, 2.0).is_err());
         assert!(CovarianceMatching::new(5, -1.0, 2.0).is_err());
         assert!(CovarianceMatching::new(5, 3.0, 2.0).is_err());
+        let matcher = || CovarianceMatching::new(5, 1.0, 2.0).unwrap();
+        assert!(matcher().with_relaxation(0.0).is_err());
+        assert!(matcher().with_relaxation(1.5).is_err());
+        assert!(matcher().with_relaxation(-0.5).is_err());
+        assert!(matcher().with_relaxation(f64::NAN).is_err());
+        assert!(matcher().with_relaxation(1.0).is_ok());
+        assert!(matcher().with_relaxation(0.1).is_ok());
         assert!(InnovationNoise::<2>::new(0, 0.5).is_err());
         assert!(InnovationNoise::<2>::new(5, 1.5).is_err());
         assert!(NisInflation::new(1.0, 2.0, 0.5, 10.0).is_err());

tests/variants.rs

@@ -179,10 +179,14 @@ fn adaptive_noise_against_an_unmodeled_acceleration() {
     let adaptive_mean_nees = adaptive_nees.iter().sum::<f64>() / adaptive_nees.len() as f64;
     let mean_scale = final_scales.iter().sum::<f64>() / final_scales.len() as f64;
 
-    // Measured: position RMS 1.791 plain against 0.938 adaptive, mean NEES 30.7
-    // against 9.5, and the scale settled at a mean of 752. The NEES is the number
+    // Measured: position RMS 1.791 plain against 0.845 adaptive, mean NEES 30.7
+    // against 4.2, and the scale settled at a mean of 356. The NEES is the number
     // that matters: both filters are inconsistent on a fault neither models, but
-    // the plain one is wrong about its uncertainty by a further factor of three.
+    // the plain one is wrong about its uncertainty by a further factor of seven.
+    //
+    // Those figures are at the default relaxation. Undamped they were 0.938 and
+    // 9.5 — the accuracy barely moved and the consistency more than halved, which
+    // is why the default is what it is.
     assert!(
         adaptive_rms * 1.5 < plain_rms,
         "adaptive RMS {adaptive_rms:.3} against plain {plain_rms:.3}"
@@ -670,25 +674,52 @@ fn fading_memory_against_a_drift() {
     );
 }
 
-/// The covariance match reaches a tenfold process noise sooner than the NIS
-/// inflater does.
+/// Both adaptive mechanisms respond to a step change in the dynamics, and the
+/// covariance match's scale has no fixed baseline to measure that against.
 ///
-/// This test was written expecting the opposite, on the reasoning that a
-/// threshold-triggered rule reacts to one bad step while a windowed one has to
-/// wait for its window. Measured over 120 runs of a step change in the dynamics:
-/// the covariance match takes a mean of 2.98 steps, the NIS inflater 7.34.
+/// This test previously claimed the covariance match "reaches a tenfold scale
+/// sooner", on a measured mean of 2.98 steps against the inflater's 7.34. That
+/// was an artifact and the claim has been withdrawn. The nominal acceleration
+/// intensity was 0.01 against a one-meter position sensor, which makes
+/// `tr(H Q0 H')` about 4e-4 against `tr(R)` of 1 — so the solved scale is
+/// `tr(R) / tr(H Q0 H')` times as sensitive to numerator noise as it is to the
+/// signal, and it sat above tenfold in 82% of runs *before the fault began*. The
+/// "delay" being short was mostly the scale already being there.
 ///
-/// The reasoning was wrong about what each rule does per step. The match *solves*
-/// for the scale, so a window holding a few large innovations puts it where it
-/// needs to be in one update, window or not. The inflater multiplies by its
-/// inflation factor per exceedance, so at a factor of 2 it needs `log2(10)`,
-/// better than three consecutive exceedances to reach tenfold, and its decay gives
-/// ground back on every step that does not exceed. Which is quicker is therefore a
-/// property of the inflation factor and the fault size, not of windowed against
-/// thresholded, and the direction can be reversed by changing either.
+/// Measured across nominal intensities, 120 runs each, scale level sampled at the
+/// onset step with nothing yet wrong:
+///
+/// | intensity | mean scale, no fault | above 10 pre-fault | match delay | inflater delay |
+/// |---|---|---|---|---|
+/// | 0.01 | 228.3 | 82% | 0.70 | 7.34 |
+/// | 0.1 | 27.8 | 48% | 2.93 | 8.93 |
+/// | 1.0 | 4.84 | 12% | 11.88 | never |
+/// | 10.0 | 1.82 | 0% | 17.00 | never |
+///
+/// So which mechanism looks faster is set by the nominal `Q`'s visibility and by
+/// the threshold chosen, not by a property of the mechanisms. At a
+/// well-conditioned intensity the inflater never reaches tenfold at all — not
+/// because it fails, but because a larger `Q` absorbs the fault and the
+/// normalized innovation rarely exceeds its chi-squared threshold.
+///
+/// What is left to assert is the part that holds: at a well-conditioned nominal
+/// `Q`, both mechanisms move the process noise up after the fault. Measured as the
+/// ratio of the mean scale after the fault to the mean before it, so that the
+/// baseline cancels: the covariance match rises 5.79x and the NIS inflater 1.51x.
+///
+/// The inflater's modest rise is not a failure of the inflater, and it is the
+/// interaction worth taking away from this test. A nominal `Q` large enough to
+/// make the covariance match well-conditioned is also large enough to absorb much
+/// of the fault, so the normalized innovation stays inside the inflater's
+/// chi-squared threshold and it has little to react to. The two mechanisms want
+/// opposite things from the nominal `Q`, and no single scenario shows both at
+/// their best.
 #[test]
-fn the_covariance_match_reaches_a_tenfold_scale_sooner() {
+fn both_adaptive_mechanisms_respond_to_a_step_change() {
+    const ONSET: usize = 60;
+
     struct Stepping {
+        intensity: f64,
         onset: f64,
         acceleration: f64,
     }
@@ -704,26 +735,29 @@ fn the_covariance_match_reaches_a_tenfold_scale_sooner() {
                 state[1] + acceleration * step,
             )
         }
-        fn noise(
-            &self,
-            _time: f64,
-            step: f64,
-            state: &nalgebra::SVector<f64, 2>,
-        ) -> SMatrix<f64, 2, 2> {
-            ConstantVelocity { intensity: 0.01 }.noise(0.0, step, state)
+        fn noise(&self, time: f64, step: f64, state: &SVector<f64, 2>) -> SMatrix<f64, 2, 2> {
+            ConstantVelocity {
+                intensity: self.intensity,
+            }
+            .noise(time, step, state)
         }
     }
 
-    let onset_step = 60usize;
+    // Well conditioned: at this intensity the match's no-fault baseline is 4.8
+    // rather than 228, so a rise above it means something.
+    const INTENSITY: f64 = 1.0;
     let truth_process = Stepping {
-        onset: onset_step as f64 * STEP,
+        intensity: INTENSITY,
+        onset: ONSET as f64 * STEP,
         acceleration: 2.0,
     };
-    let filter_process = ConstantVelocity { intensity: 0.01 };
+    let filter_process = ConstantVelocity {
+        intensity: INTENSITY,
+    };
     let measurement = Position { deviation: 1.0 };
 
-    let mut matching_delays = Vec::new();
-    let mut inflation_delays = Vec::new();
+    let mut matching_rise = Vec::new();
+    let mut inflation_rise = Vec::new();
 
     for seed in 0..SEEDS {
         let mut rng = Rng::new(0x2b2b_0000 + seed as u64);
@@ -743,8 +777,10 @@ fn the_covariance_match_reaches_a_tenfold_scale_sooner() {
         let mut inflating =
             Ukf::new(walk_prior()).with_adaptation(NisInflation::new(0.99, 2.0, 0.9, 1e6).unwrap());
 
-        let mut matching_delay = None;
-        let mut inflation_delay = None;
+        // The mean scale over the quiet stretch, against the mean over the faulted
+        // one. A ratio rather than a threshold, so the baseline cancels.
+        let mut quiet = (0.0, 0.0, 0usize);
+        let mut faulted = (0.0, 0.0, 0usize);
         for step in 1..run.len() {
             let time = run.times[step - 1];
             matching.predict(&filter_process, time, STEP).unwrap();
@@ -756,32 +792,37 @@ fn the_covariance_match_reaches_a_tenfold_scale_sooner() {
                 .update(&measurement, run.times[step], &run.measurements[step])
                 .unwrap();
 
-            if step >= onset_step {
-                if matching_delay.is_none() && matching.adaptation().scale() > 10.0 {
-                    matching_delay = Some(step - onset_step);
-                }
-                if inflation_delay.is_none() && inflating.adaptation().scale() > 10.0 {
-                    inflation_delay = Some(step - onset_step);
-                }
+            // Skip the 20 steps either side of the onset: the window straddles it.
+            let bucket = if (20..ONSET - 10).contains(&step) {
+                Some(&mut quiet)
+            } else if step > ONSET + 20 {
+                Some(&mut faulted)
+            } else {
+                None
+            };
+            if let Some(bucket) = bucket {
+                bucket.0 += matching.adaptation().scale();
+                bucket.1 += inflating.adaptation().scale();
+                bucket.2 += 1;
             }
         }
-        if let Some(delay) = matching_delay {
-            matching_delays.push(delay);
-        }
-        if let Some(delay) = inflation_delay {
-            inflation_delays.push(delay);
-        }
+        matching_rise.push((faulted.0 / faulted.2 as f64) / (quiet.0 / quiet.2 as f64));
+        inflation_rise.push((faulted.1 / faulted.2 as f64) / (quiet.1 / quiet.2 as f64));
     }
 
-    let mean = |values: &[usize]| values.iter().sum::<usize>() as f64 / values.len() as f64;
-    // Measured: covariance match 2.98 steps, NIS inflater 7.34, both responding
-    // in all 120 runs.
-    assert_eq!(matching_delays.len(), SEEDS);
-    assert_eq!(inflation_delays.len(), SEEDS);
+    let mean = |values: &[f64]| values.iter().sum::<f64>() / values.len() as f64;
+    // Measured: 5.79x for the match, 1.51x for the inflater. Different bars,
+    // because a nominal Q chosen to condition the match well leaves the inflater
+    // with little to trigger on — see the note above. The ordering between the two
+    // is not a stable property of the mechanisms.
     assert!(
-        mean(&matching_delays) < mean(&inflation_delays),
-        "matching {:.2} steps, inflation {:.2} steps",
-        mean(&matching_delays),
-        mean(&inflation_delays)
+        mean(&matching_rise) > 3.0,
+        "the covariance match should raise its scale after the fault: {:.2}x",
+        mean(&matching_rise)
+    );
+    assert!(
+        mean(&inflation_rise) > 1.2,
+        "the NIS inflater should raise its scale after the fault: {:.2}x",
+        mean(&inflation_rise)
     );
 }