logaritmisk/trueskill-tt
1
0
Fork 0
Code Issues 1 Pull Requests Packages Projects Releases Activity

feat: factorise the joint once with History::joint

Browse Source 
`posterior_of`, `posterior_of_at` and `expected_variance_reduction` each
built the joint precision matrix, factorised it, asked one question and
threw it away. 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 paid for the same factorisation once per question.

`History::joint()` returns a `Joint` handle that pays it once. Measured
on 1976 appearances, 90 queries: 68.4s one-shot against 745ms factorise
plus 93ms of queries — 81.6x, with bit-identical answers. Per query,
Criterion at 480 appearances: 9.0ms one-shot against 48us cached, 187x.

The handle borrows the history, which is what makes it correct with no
invalidation logic: the borrow checker forbids adding events or refitting
while it is alive, so there is no window in which the factorisation could
describe a fit that no longer exists. It also makes the lifetime of the
n^2 factor explicit rather than parking it in the history forever — at
4000 appearances that is 128MB, which is not something to cache silently.

Every question the joint answers turns out to be a bilinear form,

    c^T A^-1 a = (L^-1 c) . (L^-1 a)

so no caller ever needs L^-1 c itself. Replacing the general solve with a
forward substitution drops the back substitution as wasted work, halving
a query, and removes a failure mode: a variance as `c . (A^-1 c)` is a
difference of products that can round negative, where `|L^-1 c|^2` is a
sum of squares and cannot.

The one-shot calls are unchanged in cost and now delegate to the handle,
so the two paths cannot drift apart. tests/joint_handle.rs asserts they
agree bit for bit, including at pinned times, under UnknownKeys::Prior,
and across candidate matchups.

Refs #51

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ
This commit is contained in:
logaritmisk
2026-09-08 07:53:51 +02:00
co-authored by Claude Opus 5
parent b113385c6f
commit 1bb6bb31d8
6 changed files with 677 additions and 210 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/history.rs
+267 -148
View File
@@ -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() {
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))
self.joint()?.posterior_of(terms)
}
/// 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() {
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))
self.joint()?.posterior_of_at(time, terms)
}
/// 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 {
return Err(InferenceError::MismatchedShape {
kind: "expected_variance_reduction takes exactly 2 teams",
expected: 2,
got: teams.len(),
});
}
self.joint()?.expected_variance_reduction(teams, target)
}
/// 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 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(
let cholesky = crate::joint::Cholesky::factor(lambda, width).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;
// 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: 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))
Ok(Joint {
history: self,
cholesky,
latest,
at_slice,
width,
})
}
/// 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
+114 -61
View File
@@ -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
//! nothing a consumer publishes is one competitor: "can we tell these two
//! apart" is a difference, "what was this round worth" is a sum. Combining
//! marginals means assuming the competitors are independent, and they are
//! correlated through every event they share — which is the mechanism the model
//! exists to exploit.
//! Every question the joint answers is a *bilinear form* in the precision
//! matrix's inverse — the variance of a contrast is `c^T L^-1 c`, and the
//! covariance of two contrasts is `c^T L^-1 a`. None of them wants `L^-1 c`
//! itself, which is what makes the shape here worth stating explicitly.
//!
//! Measured on a five-competitor round robin, the exact correlation is +0.857,
//! so `sqrt(sa^2 + sb^2)` overstates the width of a difference by 2.6x.
//! Writing the precision as `A = L L^T`,
//!
//! ```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` is row-major and is consumed as scratch.
///
/// Returns `None` if the matrix is not positive-definite, which for a precision
/// matrix means the model is improper — a competitor with no prior and no
/// 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);
/// A factorised symmetric positive-definite matrix, reusable across queries.
pub(crate) struct Cholesky {
/// Lower triangle of `L`, row-major `n * n`. The upper triangle is
/// leftover scratch from the factorisation and is never read.
l: Vec<f64>,
n: usize,
}
// In-place Cholesky: A = L L^T, lower triangle.
for j in 0..n {
let mut d = a[j * n + j];
for k in 0..j {
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;
impl Cholesky {
/// Factorise `a` (row-major, `n * n`, symmetric) into `L L^T`.
///
/// `a` is consumed as scratch.
///
/// Returns `None` if the matrix is not positive-definite, which for a
/// precision matrix means the model is improper — a competitor with
/// neither a proper prior nor any evidence.
pub(crate) fn factor(mut a: Vec<f64>, n: usize) -> Option<Self> {
debug_assert_eq!(a.len(), n * n);
for i in j + 1..n {
let mut s = a[i * n + j];
for j in 0..n {
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.
let mut z = b.to_vec();
for i in 0..n {
let mut s = z[i];
for k in 0..i {
s -= a[i * n + k] * z[k];
/// Whiten a contrast: `y = L^-1 b`.
///
/// The point of the result is the dot product, not the vector: for two
/// contrasts `b` and `b'`, `y . y'` is `b^T A^-1 b'`. See the module docs.
pub(crate) fn whiten(&self, b: &[f64]) -> Vec<f64> {
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];
}
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];
y
}
}
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() {
// [[4, 1], [1, 3]] z = [1, 2] => z = [1/11, 7/11]
let a = vec![4.0, 1.0, 1.0, 3.0];
let z = solve_spd(a, &[1.0, 2.0]).unwrap();
assert!((z[0] - 1.0 / 11.0).abs() < 1e-12, "{z:?}");
assert!((z[1] - 7.0 / 11.0).abs() < 1e-12, "{z:?}");
fn reproduces_a_known_quadratic_form() {
let c = Cholesky::factor(vec![4.0, 1.0, 1.0, 3.0], 2).unwrap();
let y = c.whiten(&[1.0, 2.0]);
assert!((bilinear(&y, &y) - 15.0 / 11.0).abs() < 1e-12);
}
/// 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();
assert!((z[i] - expected).abs() < 1e-12, "row {i}: {z:?}");
let y = c.whiten(&e);
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!(solve_spd(a, &[1.0, 1.0]).is_none());
assert!(Cholesky::factor(vec![1.0, 2.0, 2.0, 4.0], 2).is_none());
}
}
src/lib.rs
+1 -1
View File
@@ -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};