Merge branch 'feat/joint-handle'

Factorise the joint once with History::joint, so a batch of queries pays
the O(n^3) Cholesky once rather than once per question.

Closes #51

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ

Cargo.toml

@@ -79,3 +79,7 @@ debug = true
 
 [profile.dev]
 debug = true
+
+[[bench]]
+name = "joint"
+harness = false

benches/joint.rs

@@ -0,0 +1,71 @@
+//! Cost of the joint posterior: factorising versus querying.
+//!
+//! The split is the whole point of `History::joint`. Factorising is `O(n^3)` in
+//! the history's appearances and depends only on the fit; a query is `O(n^2)`
+//! and depends only on the question. `posterior_of_one_shot` pays both every
+//! time, `joint_query` pays only the second.
+
+use criterion::{Criterion, criterion_group, criterion_main};
+use smallvec::smallvec;
+use trueskill_tt::{ConstantDrift, ConvergenceOptions, Event, History, Member, Outcome, Team};
+
+/// 30 slices of 8 duels: 480 appearances over 100 competitors.
+fn fitted() -> History<i64, ConstantDrift, trueskill_tt::NullObserver, String> {
+    let mut h: History<i64, ConstantDrift, _, String> = History::builder_with_key()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.05))
+        .convergence(ConvergenceOptions {
+            max_iter: 30,
+            epsilon: 1e-10,
+            alpha: 1.0,
+        })
+        .build();
+
+    let mut events: Vec<Event<i64, String>> = Vec::new();
+    let mut k = 0usize;
+    for t in 0..30i64 {
+        for _ in 0..8 {
+            k += 1;
+            events.push(Event {
+                time: t,
+                teams: smallvec![
+                    Team::with_members([Member::new(format!("p{}", k % 100))]),
+                    Team::with_members([Member::new(format!("p{}", (k + 37) % 100))]),
+                ],
+                outcome: Outcome::scores([
+                    (k as f64 * 0.3).sin().abs() * 20.0,
+                    (k as f64 * 0.3).cos().abs() * 20.0,
+                ]),
+            });
+        }
+    }
+    h.add_events(events).unwrap();
+    let _ = h.converge().unwrap();
+    h
+}
+
+fn bench_joint(c: &mut Criterion) {
+    let h = fitted();
+    let a = "p0".to_string();
+    let b = "p1".to_string();
+    let terms = [(&a, 1.0), (&b, -1.0)];
+
+    c.bench_function("joint_factorise_480_appearances", |bencher| {
+        bencher.iter(|| std::hint::black_box(h.joint().unwrap().variables()));
+    });
+
+    c.bench_function("posterior_of_one_shot_480_appearances", |bencher| {
+        bencher.iter(|| std::hint::black_box(h.posterior_of(&terms).unwrap()));
+    });
+
+    let joint = h.joint().unwrap();
+    c.bench_function("joint_query_480_appearances", |bencher| {
+        bencher.iter(|| std::hint::black_box(joint.posterior_of(&terms).unwrap()));
+    });
+}
+
+criterion_group!(benches, bench_joint);
+criterion_main!(benches);

src/history.rs

@@ -950,6 +950,14 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
     /// per-day or per-event slices, competitors are rarely all present in any
     /// one of them. Use [`History::posterior_of_at`] to pin a time instead.
     ///
+    /// # Asking more than one question
+    ///
+    /// This factorises the joint, uses it once, and throws it away. The
+    /// factorisation is the expensive part and it depends only on the fit, so
+    /// asking `n` questions this way pays for it `n` times. Take a
+    /// [`Joint`] with [`History::joint`] instead — the answers are identical,
+    /// and only the first one pays.
+    ///
     /// # Limitations
     ///
     /// Exact only for a history whose events are all scored, because a scored
@@ -959,10 +967,6 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
     /// so a history containing ranked events returns `JointUnavailable` rather
     /// than a plausible wrong number.
     ///
-    /// Cost is a dense solve over the history's *appearances*, not its
-    /// competitors: a competitor contributes one variable per slice it appears
-    /// in, minus any consecutive pair with no drift between them.
-    ///
     /// # Errors
     ///
     /// `UnknownKey` for a competitor the history has never seen, and
@@ -972,42 +976,7 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
     where
         K: std::fmt::Debug,
     {
-        if self.time_slices.is_empty() {
+        self.joint()?.posterior_of(terms)
-            return Err(InferenceError::JointUnavailable {
-                reason: "the history has no events",
-            });
-        }
-        if !self.time_slices.iter().all(TimeSlice::all_scored) {
-            return Err(InferenceError::JointUnavailable {
-                reason: "the history contains ranked events, whose EP factors are \
-                         not retained after convergence",
-            });
-        }
-
-        let TimeExpanded {
-            lambda,
-            latest,
-            width,
-            ..
-        } = self.time_expanded_joint();
-
-        let ResolvedTerms {
-            contrast,
-            unseen,
-            mean,
-        } = self.resolve_terms(terms, width, |index| latest.get(&index).copied())?;
-
-        let z =
-            crate::joint::solve_spd(lambda, &contrast).ok_or(InferenceError::JointUnavailable {
-                reason: "the precision matrix is not positive-definite, which means \
-                         a competitor has neither a proper prior nor any evidence",
-            })?;
-
-        let prior_var = self.sigma * self.sigma;
-        let variance: f64 = contrast.iter().zip(&z).map(|(c, z)| c * z).sum::<f64>()
-            + unseen.values().map(|c| c * c * prior_var).sum::<f64>();
-
-        Ok(Gaussian::from_mv(mean, variance))
     }
 
     /// Posterior of a linear combination, read as of `time`.
@@ -1018,6 +987,9 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
     /// stand at the end of last season" — rather than to wherever each
     /// competitor was last seen.
     ///
+    /// As with [`History::posterior_of`], this factorises the joint for one
+    /// question; [`History::joint`] amortises that across many.
+    ///
     /// # Errors
     ///
     /// As [`History::posterior_of`], plus `UnknownKey` for a competitor with no
@@ -1026,54 +998,7 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
     where
         K: std::fmt::Debug,
     {
-        if self.time_slices.is_empty() {
+        self.joint()?.posterior_of_at(time, terms)
-            return Err(InferenceError::JointUnavailable {
-                reason: "the history has no events",
-            });
-        }
-        if !self.time_slices.iter().all(TimeSlice::all_scored) {
-            return Err(InferenceError::JointUnavailable {
-                reason: "the history contains ranked events, whose EP factors are \
-                         not retained after convergence",
-            });
-        }
-
-        let TimeExpanded {
-            lambda,
-            at_slice,
-            width,
-            ..
-        } = self.time_expanded_joint();
-
-        // Latest appearance at or before `time`, per competitor.
-        let mut as_of: HashMap<Index, (usize, usize)> = HashMap::new();
-        for (slice_idx, slice) in self.time_slices.iter().enumerate() {
-            if slice.time > time {
-                break;
-            }
-            for (agent, _) in slice.appearances() {
-                if let Some(row) = at_slice.get(&(agent, slice_idx)) {
-                    as_of.insert(agent, (*row, slice_idx));
-                }
-            }
-        }
-
-        let ResolvedTerms {
-            contrast,
-            unseen,
-            mean,
-        } = self.resolve_terms(terms, width, |index| as_of.get(&index).copied())?;
-
-        let z =
-            crate::joint::solve_spd(lambda, &contrast).ok_or(InferenceError::JointUnavailable {
-                reason: "the precision matrix is not positive-definite",
-            })?;
-
-        let prior_var = self.sigma * self.sigma;
-        let variance: f64 = contrast.iter().zip(&z).map(|(c, z)| c * z).sum::<f64>()
-            + unseen.values().map(|c| c * c * prior_var).sum::<f64>();
-
-        Ok(Gaussian::from_mv(mean, variance))
     }
 
     /// How much observing this matchup would shrink the variance of `target`.
@@ -1089,6 +1014,11 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
     /// score. It is also far cheaper: one linear solve rather than a full
     /// inference pass per possible outcome.
     ///
+    /// Scoring a field of candidates is the whole point of this call, and each
+    /// candidate is one question against an unchanged fit — so use
+    /// [`Joint::expected_variance_reduction`] for anything past a single
+    /// candidate, or pay for the factorisation once per candidate.
+    ///
     /// # There is no expectation to take
     ///
     /// Observing a scored event is a rank-one update to the precision matrix,
@@ -1118,14 +1048,56 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
     where
         K: std::fmt::Debug,
     {
-        if teams.len() != 2 {
+        self.joint()?.expected_variance_reduction(teams, target)
-            return Err(InferenceError::MismatchedShape {
+    }
-                kind: "expected_variance_reduction takes exactly 2 teams",
-                expected: 2,
-                got: teams.len(),
-            });
-        }
 
+    /// Factorise the joint posterior once, to answer many questions against it.
+    ///
+    /// [`History::posterior_of`] and its neighbours each build and factorise
+    /// the joint, use it once, and drop it. The factorisation is `O(n^3)` in
+    /// the history's *appearances* and depends only on the fit, so a caller
+    /// asking about every pair in a standings table, every cell in a grid, or
+    /// every candidate in an active-learning sweep pays for the same
+    /// factorisation once per question.
+    ///
+    /// A `Joint` pays it once. Each subsequent query is `O(n^2)` — one forward
+    /// substitution — and returns exactly what the one-shot call would.
+    ///
+    /// ```
+    /// # use smallvec::smallvec;
+    /// # use trueskill_tt::{Event, History, Member, Outcome, Team};
+    /// # let mut h = History::builder().score_sigma(1.0).build();
+    /// # let round = |x, y, sx, sy, t| Event {
+    /// #     time: t,
+    /// #     teams: smallvec![
+    /// #         Team::with_members([Member::new(x)]),
+    /// #         Team::with_members([Member::new(y)]),
+    /// #     ],
+    /// #     outcome: Outcome::scores([sx, sy]),
+    /// # };
+    /// # h.add_events(vec![
+    /// #     round("a", "b", 3.0, 1.0, 1),
+    /// #     round("b", "c", 2.0, 1.0, 2),
+    /// # ]).unwrap();
+    /// # h.converge().unwrap();
+    /// let joint = h.joint()?;
+    /// for (a, b) in [("a", "b"), ("a", "c"), ("b", "c")] {
+    ///     let gap = joint.posterior_of(&[(&a, 1.0), (&b, -1.0)])?;
+    ///     println!("{a} - {b}: {:.3} +/- {:.3}", gap.mu(), gap.sigma());
+    /// }
+    /// # Ok::<(), trueskill_tt::InferenceError>(())
+    /// ```
+    ///
+    /// The handle borrows the history, so the borrow checker enforces what a
+    /// cache would otherwise have to invalidate: no events can be added and no
+    /// refit can run while it is alive. Drop it to release the factorisation,
+    /// which is `n^2` floats and is the largest thing this crate allocates.
+    ///
+    /// # Errors
+    ///
+    /// `JointUnavailable` if the history is empty, contains ranked events, or
+    /// yields a precision matrix that is not positive-definite.
+    pub fn joint(&self) -> Result<Joint<'_, T, D, O, K>, InferenceError> {
         if self.time_slices.is_empty() {
             return Err(InferenceError::JointUnavailable {
                 reason: "the history has no events",
@@ -1138,70 +1110,28 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
             });
         }
 
-        // The candidate matchup, expressed as the same kind of linear
-        // functional as the target.
-        let mut matchup: Vec<(&K, f64)> = Vec::new();
-        let mut noise = self.score_sigma * self.score_sigma;
-        for (team_idx, team) in teams.iter().enumerate() {
-            if team.is_empty() {
-                return Err(InferenceError::EmptyTeam { team: team_idx });
-            }
-            let sign = if team_idx == 0 { 1.0 } else { -1.0 };
-            for key in team.iter() {
-                matchup.push((*key, sign));
-                let beta = self
-                    .keys
-                    .get(*key)
-                    .map_or(self.beta, |index| self.agents[index].rating.beta);
-                noise += beta * beta;
-            }
-        }
-
         let TimeExpanded {
             lambda,
             latest,
+            at_slice,
             width,
-            ..
         } = self.time_expanded_joint();
 
-        let target = self.resolve_terms(target, width, |i| latest.get(&i).copied())?;
+        let cholesky = crate::joint::Cholesky::factor(lambda, width).ok_or(
-        let matchup = self.resolve_terms(&matchup, width, |i| latest.get(&i).copied())?;
-        let (target_contrast, target_unseen) = (target.contrast, target.unseen);
-        let (matchup_contrast, matchup_unseen) = (matchup.contrast, matchup.unseen);
-
-        // One solve: z = L^-1 a serves both inner products, since
-        // c^T L^-1 a = c^T z and a^T L^-1 a = a^T z.
-        let z = crate::joint::solve_spd(lambda, &matchup_contrast).ok_or(
             InferenceError::JointUnavailable {
-                reason: "the precision matrix is not positive-definite",
+                reason: "the precision matrix is not positive-definite, which means \
+                         a competitor has neither a proper prior nor any evidence",
             },
         )?;
 
-        let prior_var = self.sigma * self.sigma;
+        Ok(Joint {
-        // Competitors outside the history are independent, so they contribute
+            history: self,
-        // only where the same key appears in both functionals.
+            cholesky,
-        let cross_unseen: f64 = target_unseen
+            latest,
-            .iter()
+            at_slice,
-            .map(|(k, tc)| tc * matchup_unseen.get(k).copied().unwrap_or(0.0) * prior_var)
+            width,
-            .sum();
+        })
-        let self_unseen: f64 = matchup_unseen.values().map(|c| c * c * prior_var).sum();
-
-        let cross: f64 = target_contrast
-            .iter()
-            .zip(&z)
-            .map(|(c, z)| c * z)
-            .sum::<f64>()
-            + cross_unseen;
-        let matchup_var: f64 = matchup_contrast
-            .iter()
-            .zip(&z)
-            .map(|(a, z)| a * z)
-            .sum::<f64>()
-            + self_unseen;
-
-        Ok(cross * cross / (noise + matchup_var))
     }
-
     /// Predictive distribution of the score margin between two teams.
     ///
     /// Answers "what will the gap be, and how wide is that interval" for a
@@ -1980,6 +1910,195 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
     }
 }
 
+/// A factorised joint posterior, reusable across many queries.
+///
+/// Built by [`History::joint`]. Every question the joint answers — the width of
+/// a contrast, the covariance of two, how much a candidate matchup would
+/// sharpen either — is a bilinear form in the inverse precision matrix, and all
+/// of them share one factorisation. That factorisation is the whole cost:
+/// `O(n^3)` in the history's appearances to build, `O(n^2)` per question after.
+///
+/// The handle borrows the history, so no refit can run and no events can be
+/// added while it is alive. That is what makes it correct without any
+/// invalidation logic: there is no window in which the factorisation could
+/// describe a fit that no longer exists.
+pub struct Joint<'h, T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> {
+    history: &'h History<T, D, O, K>,
+    cholesky: crate::joint::Cholesky,
+    /// `(row, slice)` of each competitor's latest appearance.
+    latest: HashMap<Index, (usize, usize)>,
+    /// Row of each `(competitor, slice)` appearance.
+    at_slice: HashMap<(Index, usize), usize>,
+    /// Side length of the precision matrix.
+    width: usize,
+}
+
+/// Deliberately does not print the factorisation, which is `n^2` floats and
+/// would make a `{:?}` of a large joint unreadable and slow.
+impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> std::fmt::Debug
+    for Joint<'_, T, D, O, K>
+{
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        f.debug_struct("Joint")
+            .field("variables", &self.width)
+            .finish_non_exhaustive()
+    }
+}
+
+impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> Joint<'_, T, D, O, K> {
+    /// Number of variables in the joint: the history's appearances, after
+    /// collapsing consecutive pairs a competitor does not drift between.
+    ///
+    /// This is what the cost scales in, and it is not the competitor count — a
+    /// competitor contributes one variable per slice it appears in. Worth
+    /// checking before asking for a joint over a long history.
+    #[must_use]
+    pub fn variables(&self) -> usize {
+        self.width
+    }
+
+    /// Turn a resolved functional into its posterior.
+    ///
+    /// The variance is `|L^-1 c|^2` over the competitors the history knows,
+    /// plus an independent prior variance for each competitor it does not —
+    /// unseen competitors are uncorrelated with everything by construction.
+    fn distribution(&self, resolved: &ResolvedTerms) -> Gaussian {
+        let y = self.cholesky.whiten(&resolved.contrast);
+        let prior_var = self.history.sigma * self.history.sigma;
+        let variance = crate::joint::bilinear(&y, &y)
+            + resolved
+                .unseen
+                .values()
+                .map(|c| c * c * prior_var)
+                .sum::<f64>();
+        Gaussian::from_mv(resolved.mean, variance)
+    }
+
+    /// Posterior of a linear combination of competitors' skills.
+    ///
+    /// Identical to [`History::posterior_of`], including which appearance each
+    /// competitor is read at, without re-paying the factorisation.
+    ///
+    /// # Errors
+    ///
+    /// `UnknownKey` for a competitor the history has never seen.
+    pub fn posterior_of(&self, terms: &[(&K, f64)]) -> Result<Gaussian, InferenceError>
+    where
+        K: std::fmt::Debug,
+    {
+        let resolved = self
+            .history
+            .resolve_terms(terms, self.width, |index| self.latest.get(&index).copied())?;
+        Ok(self.distribution(&resolved))
+    }
+
+    /// Posterior of a linear combination, read as of `time`.
+    ///
+    /// Identical to [`History::posterior_of_at`] without re-paying the
+    /// factorisation.
+    ///
+    /// # Errors
+    ///
+    /// `UnknownKey` for a competitor with no appearance at or before `time`.
+    pub fn posterior_of_at(&self, time: T, terms: &[(&K, f64)]) -> Result<Gaussian, InferenceError>
+    where
+        K: std::fmt::Debug,
+    {
+        let as_of = self.rows_as_of(time);
+        let resolved = self
+            .history
+            .resolve_terms(terms, self.width, |index| as_of.get(&index).copied())?;
+        Ok(self.distribution(&resolved))
+    }
+
+    /// Latest appearance at or before `time`, per competitor.
+    fn rows_as_of(&self, time: T) -> HashMap<Index, (usize, usize)> {
+        let mut as_of: HashMap<Index, (usize, usize)> = HashMap::new();
+        for (slice_idx, slice) in self.history.time_slices.iter().enumerate() {
+            if slice.time > time {
+                break;
+            }
+            for (agent, _) in slice.appearances() {
+                if let Some(row) = self.at_slice.get(&(agent, slice_idx)) {
+                    as_of.insert(agent, (*row, slice_idx));
+                }
+            }
+        }
+        as_of
+    }
+
+    /// How much observing this matchup would shrink the variance of `target`.
+    ///
+    /// Identical to [`History::expected_variance_reduction`] without re-paying
+    /// the factorisation, which is the shape this call is normally used in:
+    /// one target, a field of candidate matchups, one unchanged fit.
+    ///
+    /// # Errors
+    ///
+    /// `MismatchedShape` unless exactly two teams are supplied, `EmptyTeam` for
+    /// an empty one, and `UnknownKey` for an unseen competitor.
+    pub fn expected_variance_reduction(
+        &self,
+        teams: &[&[&K]],
+        target: &[(&K, f64)],
+    ) -> Result<f64, InferenceError>
+    where
+        K: std::fmt::Debug,
+    {
+        if teams.len() != 2 {
+            return Err(InferenceError::MismatchedShape {
+                kind: "expected_variance_reduction takes exactly 2 teams",
+                expected: 2,
+                got: teams.len(),
+            });
+        }
+
+        // The candidate matchup, expressed as the same kind of linear
+        // functional as the target.
+        let mut matchup: Vec<(&K, f64)> = Vec::new();
+        let mut noise = self.history.score_sigma * self.history.score_sigma;
+        for (team_idx, team) in teams.iter().enumerate() {
+            if team.is_empty() {
+                return Err(InferenceError::EmptyTeam { team: team_idx });
+            }
+            let sign = if team_idx == 0 { 1.0 } else { -1.0 };
+            for key in team.iter() {
+                matchup.push((*key, sign));
+                let beta = self
+                    .history
+                    .keys
+                    .get(*key)
+                    .map_or(self.history.beta, |index| {
+                        self.history.agents[index].rating.beta
+                    });
+                noise += beta * beta;
+            }
+        }
+
+        let row_for = |index: Index| self.latest.get(&index).copied();
+        let target = self.history.resolve_terms(target, self.width, row_for)?;
+        let matchup = self.history.resolve_terms(&matchup, self.width, row_for)?;
+
+        let y_target = self.cholesky.whiten(&target.contrast);
+        let y_matchup = self.cholesky.whiten(&matchup.contrast);
+
+        let prior_var = self.history.sigma * self.history.sigma;
+        // Competitors outside the history are independent, so they contribute
+        // only where the same key appears in both functionals.
+        let cross_unseen: f64 = target
+            .unseen
+            .iter()
+            .map(|(k, tc)| tc * matchup.unseen.get(k).copied().unwrap_or(0.0) * prior_var)
+            .sum();
+        let self_unseen: f64 = matchup.unseen.values().map(|c| c * c * prior_var).sum();
+
+        let cross = crate::joint::bilinear(&y_target, &y_matchup) + cross_unseen;
+        let matchup_var = crate::joint::bilinear(&y_matchup, &y_matchup) + self_unseen;
+
+        Ok(cross * cross / (noise + matchup_var))
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use approx::assert_ulps_eq;

src/joint.rs

@@ -1,99 +1,152 @@
-//! Posterior of a linear combination of competitors.
+//! Cholesky factorisation of a joint precision matrix.
 //!
-//! Every accessor on `History` returns a per-competitor marginal, and almost
+//! Every question the joint answers is a *bilinear form* in the precision
-//! nothing a consumer publishes is one competitor: "can we tell these two
+//! matrix's inverse — the variance of a contrast is `c^T L^-1 c`, and the
-//! apart" is a difference, "what was this round worth" is a sum. Combining
+//! covariance of two contrasts is `c^T L^-1 a`. None of them wants `L^-1 c`
-//! marginals means assuming the competitors are independent, and they are
+//! itself, which is what makes the shape here worth stating explicitly.
-//! correlated through every event they share — which is the mechanism the model
-//! exists to exploit.
 //!
-//! Measured on a five-competitor round robin, the exact correlation is +0.857,
+//! Writing the precision as `A = L L^T`,
-//! so `sqrt(sa^2 + sb^2)` overstates the width of a difference by 2.6x.
+//!
+//! ```text
+//! c^T A^-1 a = c^T L^-T L^-1 a = (L^-1 c) . (L^-1 a)
+//! ```
+//!
+//! so a single forward substitution per contrast answers everything, and the
+//! back substitution a general solve would do is wasted work. That halves the
+//! cost of a query, and it removes a failure mode: a variance computed as
+//! `c . (A^-1 c)` is a difference of products that can round to a small
+//! negative number, where the same quantity as `|L^-1 c|^2` is a sum of
+//! squares and cannot.
+//!
+//! Factorising is `O(n^3)` and whitening is `O(n^2)`, so the split also
+//! matters structurally: the expensive half depends only on the fit, and is
+//! shared across every query a [`Joint`](crate::Joint) answers.
 
-/// Solve `A z = b` for a symmetric positive-definite `A`, by Cholesky.
+/// A factorised symmetric positive-definite matrix, reusable across queries.
-///
+pub(crate) struct Cholesky {
-/// `a` is row-major and is consumed as scratch.
+    /// Lower triangle of `L`, row-major `n * n`. The upper triangle is
-///
+    /// leftover scratch from the factorisation and is never read.
-/// Returns `None` if the matrix is not positive-definite, which for a precision
+    l: Vec<f64>,
-/// matrix means the model is improper — a competitor with no prior and no
+    n: usize,
-/// evidence.
+}
-pub(crate) fn solve_spd(mut a: Vec<f64>, b: &[f64]) -> Option<Vec<f64>> {
-    let n = b.len();
-    debug_assert_eq!(a.len(), n * n);
 
-    // In-place Cholesky: A = L L^T, lower triangle.
+impl Cholesky {
-    for j in 0..n {
+    /// Factorise `a` (row-major, `n * n`, symmetric) into `L L^T`.
-        let mut d = a[j * n + j];
+    ///
-        for k in 0..j {
+    /// `a` is consumed as scratch.
-            d -= a[j * n + k] * a[j * n + k];
+    ///
-        }
+    /// Returns `None` if the matrix is not positive-definite, which for a
-        // Explicit rather than `!(d > 0.0)`: a NaN pivot must fail here too,
+    /// precision matrix means the model is improper — a competitor with
-        // and a negated comparison would let it through as "not positive".
+    /// neither a proper prior nor any evidence.
-        if d.is_nan() || d <= 0.0 {
+    pub(crate) fn factor(mut a: Vec<f64>, n: usize) -> Option<Self> {
-            return None;
+        debug_assert_eq!(a.len(), n * n);
-        }
-        let d = d.sqrt();
-        a[j * n + j] = d;
 
-        for i in j + 1..n {
+        for j in 0..n {
-            let mut s = a[i * n + j];
+            let mut d = a[j * n + j];
             for k in 0..j {
-                s -= a[i * n + k] * a[j * n + k];
+                d -= a[j * n + k] * a[j * n + k];
+            }
+            // Explicit rather than `!(d > 0.0)`: a NaN pivot must fail here
+            // too, and a negated comparison would let it through as "not
+            // positive".
+            if d.is_nan() || d <= 0.0 {
+                return None;
+            }
+            let d = d.sqrt();
+            a[j * n + j] = d;
+
+            for i in j + 1..n {
+                let mut s = a[i * n + j];
+                for k in 0..j {
+                    s -= a[i * n + k] * a[j * n + k];
+                }
+                a[i * n + j] = s / d;
             }
-            a[i * n + j] = s / d;
         }
+
+        Some(Self { l: a, n })
     }
 
-    // Forward substitution, then back substitution.
+    /// Whiten a contrast: `y = L^-1 b`.
-    let mut z = b.to_vec();
+    ///
-    for i in 0..n {
+    /// The point of the result is the dot product, not the vector: for two
-        let mut s = z[i];
+    /// contrasts `b` and `b'`, `y . y'` is `b^T A^-1 b'`. See the module docs.
-        for k in 0..i {
+    pub(crate) fn whiten(&self, b: &[f64]) -> Vec<f64> {
-            s -= a[i * n + k] * z[k];
+        debug_assert_eq!(b.len(), self.n);
+        let n = self.n;
+        let mut y = b.to_vec();
+        for i in 0..n {
+            // Folded from `y[i]` rather than summed and subtracted once, so the
+            // accumulation order matches a plain substitution loop exactly.
+            let row = &self.l[i * n..i * n + i];
+            let s = row
+                .iter()
+                .zip(&y[..i])
+                .fold(y[i], |acc, (l, v)| acc - l * v);
+            y[i] = s / self.l[i * n + i];
         }
-        z[i] = s / a[i * n + i];
+        y
-    }
-    for i in (0..n).rev() {
-        let mut s = z[i];
-        for k in i + 1..n {
-            s -= a[k * n + i] * z[k];
-        }
-        z[i] = s / a[i * n + i];
     }
+}
 
-    Some(z)
+/// `b^T A^-1 b'`, given the two whitened contrasts.
+pub(crate) fn bilinear(y: &[f64], y_prime: &[f64]) -> f64 {
+    y.iter().zip(y_prime).map(|(a, b)| a * b).sum()
 }
 
 #[cfg(test)]
 mod tests {
     use super::*;
 
+    /// `[[4, 1], [1, 3]] z = [1, 2]` has `z = [1/11, 7/11]`, so the quadratic
+    /// form `b^T A^-1 b` is `1 * 1/11 + 2 * 7/11 = 15/11`.
     #[test]
-    fn solves_a_known_system() {
+    fn reproduces_a_known_quadratic_form() {
-        // [[4, 1], [1, 3]] z = [1, 2]  =>  z = [1/11, 7/11]
+        let c = Cholesky::factor(vec![4.0, 1.0, 1.0, 3.0], 2).unwrap();
-        let a = vec![4.0, 1.0, 1.0, 3.0];
+        let y = c.whiten(&[1.0, 2.0]);
-        let z = solve_spd(a, &[1.0, 2.0]).unwrap();
+        assert!((bilinear(&y, &y) - 15.0 / 11.0).abs() < 1e-12);
-        assert!((z[0] - 1.0 / 11.0).abs() < 1e-12, "{z:?}");
-        assert!((z[1] - 7.0 / 11.0).abs() < 1e-12, "{z:?}");
     }
 
+    /// Whitening `e_i` recovers the inverse's diagonal, which is the variance
+    /// of a single variable.
     #[test]
     fn recovers_the_inverse_diagonal() {
         // A = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]]; inverse diagonal is
         // [0.75, 1.0, 0.75].
         let a = vec![2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0];
+        let c = Cholesky::factor(a, 3).unwrap();
         for (i, expected) in [0.75, 1.0, 0.75].into_iter().enumerate() {
             let mut e = vec![0.0; 3];
             e[i] = 1.0;
-            let z = solve_spd(a.clone(), &e).unwrap();
+            let y = c.whiten(&e);
-            assert!((z[i] - expected).abs() < 1e-12, "row {i}: {z:?}");
+            assert!((bilinear(&y, &y) - expected).abs() < 1e-12, "row {i}");
         }
     }
 
+    /// The off-diagonal bilinear form is symmetric and matches the inverse.
+    #[test]
+    fn recovers_an_off_diagonal_covariance() {
+        // Same A; (A^-1)_{0,1} = 0.5.
+        let a = vec![2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0];
+        let c = Cholesky::factor(a, 3).unwrap();
+        let y0 = c.whiten(&[1.0, 0.0, 0.0]);
+        let y1 = c.whiten(&[0.0, 1.0, 0.0]);
+        assert!((bilinear(&y0, &y1) - 0.5).abs() < 1e-12);
+        assert!((bilinear(&y1, &y0) - 0.5).abs() < 1e-12);
+    }
+
+    /// A variance can never come out negative, because it is a sum of squares.
+    #[test]
+    fn a_quadratic_form_is_never_negative() {
+        let a = vec![1e12, 1e12 - 1.0, 1e12 - 1.0, 1e12];
+        let c = Cholesky::factor(a, 2).unwrap();
+        let y = c.whiten(&[1.0, -1.0]);
+        assert!(bilinear(&y, &y) >= 0.0);
+    }
+
     #[test]
     fn rejects_a_non_positive_definite_matrix() {
         // Singular: the second row is a multiple of the first.
-        let a = vec![1.0, 2.0, 2.0, 4.0];
+        assert!(Cholesky::factor(vec![1.0, 2.0, 2.0, 4.0], 2).is_none());
-        assert!(solve_spd(a, &[1.0, 1.0]).is_none());
     }
 }

src/lib.rs

@@ -141,7 +141,7 @@ pub use event::{Event, Member, Team};
 pub use event_builder::EventBuilder;
 pub use game::{Game, GameOptions, OwnedGame};
 pub use gaussian::Gaussian;
-pub use history::{History, HistoryBuilder};
+pub use history::{History, HistoryBuilder, Joint};
 pub use key_table::KeyTable;
 use matrix::Matrix;
 pub use observer::{NullObserver, Observer};

tests/joint_handle.rs

@@ -0,0 +1,220 @@
+//! `History::joint` factorises once and answers many questions.
+//!
+//! The contract that matters is *identity*: a `Joint` must return exactly what
+//! the one-shot call returns, bit for bit. A faster path that quietly disagreed
+//! with the slow one would be worse than no fast path — a caller would get
+//! different numbers depending on how many questions they happened to ask.
+
+use smallvec::smallvec;
+use trueskill_tt::{
+    ConstantDrift, ConvergenceOptions, Event, History, InferenceError, Member, Outcome, Team,
+    UnknownKeys,
+};
+
+type H = History<i64, ConstantDrift, trueskill_tt::NullObserver, &'static str>;
+
+fn duel(a: &'static str, b: &'static str, t: i64, sa: f64, sb: f64) -> Event<i64, &'static str> {
+    Event {
+        time: t,
+        teams: smallvec![
+            Team::with_members([Member::new(a)]),
+            Team::with_members([Member::new(b)]),
+        ],
+        outcome: Outcome::scores([sa, sb]),
+    }
+}
+
+fn ranked(a: &'static str, b: &'static str, t: i64) -> Event<i64, &'static str> {
+    Event {
+        time: t,
+        teams: smallvec![
+            Team::with_members([Member::new(a)]),
+            Team::with_members([Member::new(b)]),
+        ],
+        outcome: Outcome::winner(0, 2),
+    }
+}
+
+fn history(unknown: UnknownKeys) -> H {
+    History::builder()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.5))
+        .unknown_keys(unknown)
+        .convergence(ConvergenceOptions {
+            max_iter: 20_000,
+            epsilon: 1e-13,
+            alpha: 1.0,
+        })
+        .build()
+}
+
+/// Several slices, competitors with different last appearances, so `latest`
+/// and `at_slice` both have work to do.
+fn fitted(unknown: UnknownKeys) -> H {
+    let mut h = history(unknown);
+    h.add_events(vec![
+        duel("a", "b", 1, 5.0, 2.0),
+        duel("c", "d", 1, 3.0, 3.5),
+        duel("a", "c", 2, 6.0, 1.0),
+        duel("b", "d", 3, 4.0, 3.0),
+        duel("a", "d", 4, 7.0, 2.0),
+        duel("b", "c", 5, 2.0, 4.0),
+    ])
+    .unwrap();
+    let report = h.converge().unwrap();
+    assert!(report.converged, "fixture must converge");
+    h
+}
+
+const PAIRS: [(&str, &str); 6] = [
+    ("a", "b"),
+    ("a", "c"),
+    ("a", "d"),
+    ("b", "c"),
+    ("b", "d"),
+    ("c", "d"),
+];
+
+#[test]
+fn a_joint_answers_exactly_what_the_one_shot_call_does() {
+    let h = fitted(UnknownKeys::Reject);
+    let joint = h.joint().unwrap();
+
+    for (a, b) in PAIRS {
+        let terms = [(&a, 1.0), (&b, -1.0)];
+        let one_shot = h.posterior_of(&terms).unwrap();
+        let cached = joint.posterior_of(&terms).unwrap();
+        assert_eq!(one_shot.pi(), cached.pi(), "{a} - {b}");
+        assert_eq!(one_shot.tau(), cached.tau(), "{a} - {b}");
+    }
+}
+
+#[test]
+fn a_joint_agrees_at_a_pinned_time_too() {
+    let h = fitted(UnknownKeys::Reject);
+    let joint = h.joint().unwrap();
+
+    for time in 1..=5 {
+        for (a, b) in PAIRS {
+            let terms = [(&a, 1.0), (&b, -1.0)];
+            let one_shot = h.posterior_of_at(time, &terms);
+            let cached = joint.posterior_of_at(time, &terms);
+            match (one_shot, cached) {
+                (Ok(x), Ok(y)) => {
+                    assert_eq!(x.pi(), y.pi(), "t={time} {a} - {b}");
+                    assert_eq!(x.tau(), y.tau(), "t={time} {a} - {b}");
+                }
+                (Err(x), Err(y)) => assert_eq!(x, y, "t={time} {a} - {b}"),
+                (x, y) => panic!("t={time} {a} - {b}: disagreed on success: {x:?} vs {y:?}"),
+            }
+        }
+    }
+}
+
+#[test]
+fn a_joint_scores_candidate_matchups_identically() {
+    let h = fitted(UnknownKeys::Reject);
+    let joint = h.joint().unwrap();
+    let (a, b) = ("a", "b");
+    let target = [(&a, 1.0), (&b, -1.0)];
+
+    for (x, y) in PAIRS {
+        let teams: [&[&&str]; 2] = [&[&x], &[&y]];
+        let one_shot = h.expected_variance_reduction(&teams, &target).unwrap();
+        let cached = joint.expected_variance_reduction(&teams, &target).unwrap();
+        assert_eq!(one_shot, cached, "{x} vs {y}");
+    }
+}
+
+/// The whole point: a competitor appears once per slice, so the joint is over
+/// appearances rather than competitors, and a caller sizing a batch needs to
+/// know which.
+#[test]
+fn variables_counts_appearances_not_competitors() {
+    let h = fitted(UnknownKeys::Reject);
+    let joint = h.joint().unwrap();
+    // Four competitors, twelve appearances across five slices, all with
+    // positive drift between them, so no two collapse.
+    assert_eq!(joint.variables(), 12);
+}
+
+/// With `drift = 0` consecutive appearances are the same latent variable, so
+/// the joint is smaller than the appearance count.
+#[test]
+fn pinned_competitors_collapse_consecutive_appearances() {
+    let mut h = History::builder()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.0))
+        .convergence(ConvergenceOptions {
+            max_iter: 20_000,
+            epsilon: 1e-13,
+            alpha: 1.0,
+        })
+        .build();
+    h.add_events(vec![
+        duel("a", "b", 1, 5.0, 2.0),
+        duel("a", "b", 2, 4.0, 3.0),
+        duel("a", "b", 3, 6.0, 1.0),
+    ])
+    .unwrap();
+    assert!(h.converge().unwrap().converged);
+    assert_eq!(h.joint().unwrap().variables(), 2);
+}
+
+#[test]
+fn a_ranked_history_has_no_exact_joint() {
+    let mut h = history(UnknownKeys::Reject);
+    h.add_events(vec![duel("a", "b", 1, 5.0, 2.0), ranked("a", "b", 2)])
+        .unwrap();
+    let _ = h.converge().unwrap();
+    assert!(matches!(
+        h.joint().unwrap_err(),
+        InferenceError::JointUnavailable { .. }
+    ));
+}
+
+#[test]
+fn an_empty_history_has_no_joint() {
+    let h = history(UnknownKeys::Reject);
+    assert!(matches!(
+        h.joint().unwrap_err(),
+        InferenceError::JointUnavailable { .. }
+    ));
+}
+
+/// Unknown keys are decided per query, not when the joint is factorised — the
+/// factorisation does not depend on the question.
+#[test]
+fn unknown_keys_are_rejected_per_query() {
+    let h = fitted(UnknownKeys::Reject);
+    let joint = h.joint().unwrap();
+    let (a, z) = ("a", "nobody");
+    assert!(matches!(
+        joint.posterior_of(&[(&a, 1.0), (&z, -1.0)]).unwrap_err(),
+        InferenceError::UnknownKey { .. }
+    ));
+    // The handle is still usable afterwards.
+    let b = "b";
+    assert!(joint.posterior_of(&[(&a, 1.0), (&b, -1.0)]).is_ok());
+}
+
+/// Under `Prior`, an unseen competitor is independent of everything in the
+/// history, and the cached path must add the same prior variance the one-shot
+/// path does.
+#[test]
+fn unseen_competitors_match_the_one_shot_path() {
+    let h = fitted(UnknownKeys::Prior);
+    let joint = h.joint().unwrap();
+    let (a, z) = ("a", "nobody");
+    let terms = [(&a, 1.0), (&z, -1.0)];
+    let one_shot = h.posterior_of(&terms).unwrap();
+    let cached = joint.posterior_of(&terms).unwrap();
+    assert_eq!(one_shot.pi(), cached.pi());
+    assert_eq!(one_shot.tau(), cached.tau());
+}