Merge perf/sparse-joint (#52)
This commit is contained in:
@@ -47,12 +47,15 @@ harness = false
|
||||
|
||||
[dependencies]
|
||||
approx = { version = "0.5.1", optional = true }
|
||||
feral-amd = "0.2"
|
||||
libm = "0.2.16"
|
||||
rayon = { version = "1", optional = true }
|
||||
smallvec = "1"
|
||||
|
||||
[features]
|
||||
approx = ["dep:approx"]
|
||||
# Exposes the joint sparsity pattern for the #52 measurement. Test-only.
|
||||
measure-sparsity = []
|
||||
rayon = ["dep:rayon"]
|
||||
|
||||
[dev-dependencies]
|
||||
|
||||
+33
-10
@@ -386,8 +386,8 @@ pub(crate) struct CompetitorConfig {
|
||||
/// The joint precision over a history's appearances, with the maps needed to
|
||||
/// address a competitor either at their latest appearance or at a given slice.
|
||||
struct TimeExpanded {
|
||||
/// Row-major precision matrix over appearances.
|
||||
lambda: Vec<f64>,
|
||||
/// Precision matrix over appearances, accumulated sparsely.
|
||||
lambda: crate::joint::SymmetricBuilder,
|
||||
/// `(row, slice)` of each competitor's latest appearance.
|
||||
latest: HashMap<Index, (usize, usize)>,
|
||||
/// Row of each `(competitor, slice)` appearance.
|
||||
@@ -1415,6 +1415,23 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<K, T, D
|
||||
prior_variance * SQRT_EPSILON
|
||||
}
|
||||
|
||||
/// Test-only: the joint's sparsity pattern, for the #52 measurement.
|
||||
///
|
||||
/// Returns `(n, adjacency)` where `adjacency[i]` holds the off-diagonal
|
||||
/// nonzero columns of row `i`.
|
||||
#[cfg(feature = "measure-sparsity")]
|
||||
pub fn joint_pattern_for_measurement(&self) -> (usize, Vec<std::collections::HashSet<usize>>) {
|
||||
let te = self.time_expanded_joint();
|
||||
let n = te.width;
|
||||
let mut adj = vec![std::collections::HashSet::new(); n];
|
||||
for (i, j) in te.lambda.pattern() {
|
||||
if i != j {
|
||||
adj[i].insert(j);
|
||||
}
|
||||
}
|
||||
(n, adj)
|
||||
}
|
||||
|
||||
fn time_expanded_joint(&self) -> TimeExpanded {
|
||||
let mut latest: HashMap<Index, (usize, usize)> = HashMap::new();
|
||||
let mut at_slice: HashMap<(Index, usize), usize> = HashMap::new();
|
||||
@@ -1455,17 +1472,23 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<K, T, D
|
||||
}
|
||||
}
|
||||
|
||||
let mut lambda = vec![0.0; n * n];
|
||||
// Accumulated sparsely: this matrix is ~0.19% dense at scale, and the
|
||||
// dense form was 31 MB at n = 1976 and 128 MB at ustat's ~4000
|
||||
// appearances, of which 99.8% was zeros. See `joint.rs` and #52.
|
||||
let mut lambda = crate::joint::SymmetricBuilder::new();
|
||||
|
||||
for (row, competitor) in first_rows {
|
||||
lambda[row * n + row] +=
|
||||
1.0 / self.competitors[competitor].rating.prior.sigma().powi(2);
|
||||
lambda.add(
|
||||
row,
|
||||
row,
|
||||
1.0 / self.competitors[competitor].rating.prior.sigma().powi(2),
|
||||
);
|
||||
}
|
||||
for (a, b, drift) in drift_links {
|
||||
lambda[a * n + a] += 1.0 / drift;
|
||||
lambda[b * n + b] += 1.0 / drift;
|
||||
lambda[a * n + b] -= 1.0 / drift;
|
||||
lambda[b * n + a] -= 1.0 / drift;
|
||||
lambda.add(a, a, 1.0 / drift);
|
||||
lambda.add(b, b, 1.0 / drift);
|
||||
lambda.add(a, b, -1.0 / drift);
|
||||
lambda.add(b, a, -1.0 / drift);
|
||||
}
|
||||
for (slice_idx, slice) in self.time_slices.iter().enumerate() {
|
||||
for (contrast, noise) in slice.scored_contrasts(&self.competitors) {
|
||||
@@ -1473,7 +1496,7 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<K, T, D
|
||||
let ra = at_slice[&(*ia, slice_idx)];
|
||||
for (ib, cb) in &contrast {
|
||||
let rb = at_slice[&(*ib, slice_idx)];
|
||||
lambda[ra * n + rb] += ca * cb / noise;
|
||||
lambda.add(ra, rb, ca * cb / noise);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
+406
-48
@@ -1,4 +1,4 @@
|
||||
//! Cholesky factorisation of a joint precision matrix.
|
||||
//! Sparse Cholesky factorisation of a joint precision matrix.
|
||||
//!
|
||||
//! 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
|
||||
@@ -18,72 +18,324 @@
|
||||
//! 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.
|
||||
//! # Why this is sparse (#52)
|
||||
//!
|
||||
//! A time-expanded joint is *extremely* sparse and gets sparser as the history
|
||||
//! grows: a row couples only to its own previous and next appearance through
|
||||
//! the drift link, and to whoever co-appeared in its slice. Measured on a
|
||||
//! 76-slice, 988-duel, 200-competitor history: `n = 1976`, `nnz = 7504`,
|
||||
//! **0.19% dense**.
|
||||
//!
|
||||
//! This used to store all `n^2` entries and run a dense `O(n^3)` factorisation
|
||||
//! over them. Two measurements decided the replacement:
|
||||
//!
|
||||
//! - **Ordering alone does nothing to a dense factorisation.** Its inner loops
|
||||
//! run over every `k` whether or not the entry is zero. A 700x700 banded
|
||||
//! matrix at 0.43% density factorised in 30.196 ms in band order and
|
||||
//! 29.544 ms under a scramble that destroyed the band — identical, as the
|
||||
//! flop count says it must be. Fill-reducing order is worth nothing until
|
||||
//! the factorisation skips zeros.
|
||||
//! - **Together they are worth four orders of magnitude.** On that `n = 1976`
|
||||
//! fixture, against `n^3/3 = 2.572e9` flops dense: sparse in the natural
|
||||
//! order needs `5.597e7` (46x better), and sparse under an AMD fill-reducing
|
||||
//! order needs `8.656e4` — **29,710x**. AMD is worth 646x *on top of*
|
||||
//! sparsity and nothing without it.
|
||||
//!
|
||||
//! Natural ordering fills in badly here for the reason #52 predicted: a
|
||||
//! competitor who appears in slice 0 and not again until slice 75 creates a
|
||||
//! drift link spanning nearly the whole matrix. `nnz(L)` is 292,437 under the
|
||||
//! natural order against 11,583 under AMD, from an `A` with 7,504.
|
||||
//!
|
||||
//! The ordering comes from `feral-amd`. The factorisation is the up-looking
|
||||
//! sparse Cholesky of Davis's *Direct Methods for Sparse Linear Systems*,
|
||||
//! written here rather than taken from a crate: the sparse solvers on
|
||||
//! crates.io either pull SIMD dispatch (`faer`, and `feral` itself, both
|
||||
//! through `pulp`), which would make results differ between an AVX-512 host
|
||||
//! and an AVX2 one — the same class of drift the `libm`-over-`std` decision
|
||||
//! was made to avoid — or are LGPL, or disclaim fill-reduction in their own
|
||||
//! docs.
|
||||
|
||||
use std::collections::BTreeMap;
|
||||
|
||||
/// A symmetric matrix accumulated entry by entry, before factorisation.
|
||||
///
|
||||
/// A `BTreeMap` rather than a hash map because the iteration order becomes the
|
||||
/// factorisation's summation order, and a hash map's order varies per process.
|
||||
/// `tests/cross_process_determinism.rs` exists because that has bitten before.
|
||||
#[derive(Default)]
|
||||
pub(crate) struct SymmetricBuilder {
|
||||
entries: BTreeMap<(usize, usize), f64>,
|
||||
}
|
||||
|
||||
impl SymmetricBuilder {
|
||||
pub(crate) fn new() -> Self {
|
||||
Self::default()
|
||||
}
|
||||
|
||||
/// Add `value` to entry `(row, col)`. Both triangles must be supplied.
|
||||
pub(crate) fn add(&mut self, row: usize, col: usize, value: f64) {
|
||||
*self.entries.entry((row, col)).or_insert(0.0) += value;
|
||||
}
|
||||
|
||||
/// The `(row, col)` positions that hold a nonzero. For the #52 measurement.
|
||||
#[cfg(feature = "measure-sparsity")]
|
||||
pub(crate) fn pattern(&self) -> impl Iterator<Item = (usize, usize)> + '_ {
|
||||
self.entries
|
||||
.iter()
|
||||
.filter(|(_, v)| **v != 0.0)
|
||||
.map(|(&rc, _)| rc)
|
||||
}
|
||||
}
|
||||
|
||||
/// 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,
|
||||
/// `inv[old] = new`: where each original row sits after the AMD reorder.
|
||||
inv: Vec<usize>,
|
||||
/// `L` in compressed-column form, permuted. Within a column the diagonal
|
||||
/// is first and the rest ascend by row.
|
||||
col_ptr: Vec<usize>,
|
||||
row_idx: Vec<usize>,
|
||||
val: Vec<f64>,
|
||||
}
|
||||
|
||||
impl Cholesky {
|
||||
/// Factorise `a` (row-major, `n * n`, symmetric) into `L L^T`.
|
||||
///
|
||||
/// `a` is consumed as scratch.
|
||||
/// Factorise the accumulated matrix into `L L^T`, under a fill-reducing
|
||||
/// permutation.
|
||||
///
|
||||
/// 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 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];
|
||||
/// neither a proper prior nor any evidence — or if the ordering fails.
|
||||
pub(crate) fn factor(built: SymmetricBuilder, n: usize) -> Option<Self> {
|
||||
if n == 0 {
|
||||
return Some(Self {
|
||||
n: 0,
|
||||
inv: Vec::new(),
|
||||
col_ptr: vec![0],
|
||||
row_idx: Vec::new(),
|
||||
val: Vec::new(),
|
||||
});
|
||||
}
|
||||
|
||||
let inv = Self::amd_permutation(n, &built)?;
|
||||
|
||||
// Upper triangle of the permuted matrix, column-major: column `c`
|
||||
// holds the rows `r <= c`. Exactly one of a symmetric pair survives
|
||||
// the `r <= c` filter, so nothing is double-counted.
|
||||
let mut cols: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];
|
||||
for (&(old_r, old_c), &v) in &built.entries {
|
||||
if v == 0.0 {
|
||||
continue;
|
||||
}
|
||||
let (r, c) = (inv[old_r], inv[old_c]);
|
||||
if r <= c {
|
||||
cols[c].push((r, v));
|
||||
}
|
||||
}
|
||||
let mut a_ptr = Vec::with_capacity(n + 1);
|
||||
let mut a_row = Vec::new();
|
||||
let mut a_val = Vec::new();
|
||||
a_ptr.push(0usize);
|
||||
for col in &mut cols {
|
||||
col.sort_unstable_by_key(|&(r, _)| r);
|
||||
for &(r, v) in col.iter() {
|
||||
a_row.push(r);
|
||||
a_val.push(v);
|
||||
}
|
||||
a_ptr.push(a_row.len());
|
||||
}
|
||||
|
||||
let parent = Self::etree(n, &a_ptr, &a_row);
|
||||
|
||||
// Symbolic pass: how many entries each column of L will hold. Running
|
||||
// `ereach` per column costs O(nnz(L)) in total, which is the same order
|
||||
// as the numeric pass it sizes.
|
||||
let mut counts = vec![0usize; n];
|
||||
let mut stack = vec![0usize; n];
|
||||
let mut mark = vec![false; n];
|
||||
for k in 0..n {
|
||||
let top = Self::ereach(k, &a_ptr, &a_row, &parent, &mut stack, &mut mark);
|
||||
for &i in &stack[top..] {
|
||||
counts[i] += 1;
|
||||
}
|
||||
counts[k] += 1; // the diagonal
|
||||
}
|
||||
|
||||
let mut col_ptr = Vec::with_capacity(n + 1);
|
||||
col_ptr.push(0usize);
|
||||
for &c in &counts {
|
||||
col_ptr.push(col_ptr[col_ptr.len() - 1] + c);
|
||||
}
|
||||
let nnz = col_ptr[n];
|
||||
let mut row_idx = vec![0usize; nnz];
|
||||
let mut val = vec![0.0f64; nnz];
|
||||
|
||||
// `next[i]` is the slot column `i` will fill next. Column `i`'s
|
||||
// diagonal lands first, at `col_ptr[i]`, because nothing is written to
|
||||
// a column before its own iteration.
|
||||
let mut next: Vec<usize> = col_ptr[..n].to_vec();
|
||||
let mut x = vec![0.0f64; n];
|
||||
|
||||
for k in 0..n {
|
||||
let top = Self::ereach(k, &a_ptr, &a_row, &parent, &mut stack, &mut mark);
|
||||
|
||||
for p in a_ptr[k]..a_ptr[k + 1] {
|
||||
if a_row[p] <= k {
|
||||
x[a_row[p]] = a_val[p];
|
||||
}
|
||||
}
|
||||
|
||||
let mut d = x[k];
|
||||
x[k] = 0.0;
|
||||
|
||||
for &i in &stack[top..] {
|
||||
let lki = x[i] / val[col_ptr[i]];
|
||||
x[i] = 0.0;
|
||||
for p in col_ptr[i] + 1..next[i] {
|
||||
x[row_idx[p]] -= val[p] * lki;
|
||||
}
|
||||
d -= lki * lki;
|
||||
let p = next[i];
|
||||
next[i] += 1;
|
||||
row_idx[p] = k;
|
||||
val[p] = lki;
|
||||
}
|
||||
|
||||
// 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;
|
||||
}
|
||||
let p = next[k];
|
||||
next[k] += 1;
|
||||
row_idx[p] = k;
|
||||
val[p] = d.sqrt();
|
||||
}
|
||||
|
||||
Some(Self { l: a, n })
|
||||
Some(Self {
|
||||
n,
|
||||
inv,
|
||||
col_ptr,
|
||||
row_idx,
|
||||
val,
|
||||
})
|
||||
}
|
||||
|
||||
/// Whiten a contrast: `y = L^-1 b`.
|
||||
/// AMD fill-reducing order, as `inv[old] = new`.
|
||||
fn amd_permutation(n: usize, built: &SymmetricBuilder) -> Option<Vec<usize>> {
|
||||
let mut cols: Vec<Vec<i32>> = vec![Vec::new(); n];
|
||||
for (&(r, c), &v) in &built.entries {
|
||||
if v != 0.0 {
|
||||
cols[c].push(i32::try_from(r).ok()?);
|
||||
}
|
||||
}
|
||||
let mut col_ptr = Vec::with_capacity(n + 1);
|
||||
let mut row_idx = Vec::new();
|
||||
col_ptr.push(0i32);
|
||||
for (j, col) in cols.iter_mut().enumerate() {
|
||||
col.push(i32::try_from(j).ok()?);
|
||||
col.sort_unstable();
|
||||
col.dedup();
|
||||
row_idx.extend_from_slice(col);
|
||||
col_ptr.push(i32::try_from(row_idx.len()).ok()?);
|
||||
}
|
||||
|
||||
let pattern = feral_amd::CscPattern::new(n, &col_ptr, &row_idx)?;
|
||||
// `perm[new] = old`; we want the inverse.
|
||||
let perm = feral_amd::amd_order(&pattern).ok()?;
|
||||
let mut inv = vec![0usize; n];
|
||||
for (new, &old) in perm.iter().enumerate() {
|
||||
inv[usize::try_from(old).ok()?] = new;
|
||||
}
|
||||
Some(inv)
|
||||
}
|
||||
|
||||
/// Elimination tree of the upper-triangular pattern. `usize::MAX` is "no
|
||||
/// parent", i.e. a root.
|
||||
fn etree(n: usize, col_ptr: &[usize], row_idx: &[usize]) -> Vec<usize> {
|
||||
let mut parent = vec![usize::MAX; n];
|
||||
let mut ancestor = vec![usize::MAX; n];
|
||||
for k in 0..n {
|
||||
for &row in &row_idx[col_ptr[k]..col_ptr[k + 1]] {
|
||||
let mut i = row;
|
||||
while i != usize::MAX && i < k {
|
||||
let next = ancestor[i];
|
||||
ancestor[i] = k;
|
||||
if next == usize::MAX {
|
||||
parent[i] = k;
|
||||
}
|
||||
i = next;
|
||||
}
|
||||
}
|
||||
}
|
||||
parent
|
||||
}
|
||||
|
||||
/// Nonzero pattern of row `k` of `L`, written into `stack[top..n]` in
|
||||
/// topological order. Returns `top`.
|
||||
///
|
||||
/// `stack` is used from both ends — a scratch region from `0` while walking
|
||||
/// each path up the tree, and the result from `n` downwards. They cannot
|
||||
/// collide because every node is pushed at most once across the whole call.
|
||||
fn ereach(
|
||||
k: usize,
|
||||
col_ptr: &[usize],
|
||||
row_idx: &[usize],
|
||||
parent: &[usize],
|
||||
stack: &mut [usize],
|
||||
mark: &mut [bool],
|
||||
) -> usize {
|
||||
let n = mark.len();
|
||||
let mut top = n;
|
||||
mark[k] = true;
|
||||
for &row in &row_idx[col_ptr[k]..col_ptr[k + 1]] {
|
||||
let mut i = row;
|
||||
if i > k {
|
||||
continue;
|
||||
}
|
||||
let mut len = 0usize;
|
||||
while i != usize::MAX && !mark[i] {
|
||||
stack[len] = i;
|
||||
len += 1;
|
||||
mark[i] = true;
|
||||
i = parent[i];
|
||||
}
|
||||
// Reverse the path onto the output end, so the result stays in
|
||||
// topological order overall.
|
||||
while len > 0 {
|
||||
len -= 1;
|
||||
top -= 1;
|
||||
stack[top] = stack[len];
|
||||
}
|
||||
}
|
||||
for &i in &stack[top..] {
|
||||
mark[i] = false;
|
||||
}
|
||||
mark[k] = false;
|
||||
top
|
||||
}
|
||||
|
||||
/// Whiten a contrast: `y = L^-1 P 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.
|
||||
///
|
||||
/// The result is in the permuted order, and stays there — a dot product
|
||||
/// does not care, as long as both operands were permuted the same way.
|
||||
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];
|
||||
let mut y = vec![0.0f64; n];
|
||||
for (old, &v) in b.iter().enumerate() {
|
||||
y[self.inv[old]] = v;
|
||||
}
|
||||
for j in 0..n {
|
||||
y[j] /= self.val[self.col_ptr[j]];
|
||||
let yj = y[j];
|
||||
for p in self.col_ptr[j] + 1..self.col_ptr[j + 1] {
|
||||
y[self.row_idx[p]] -= self.val[p] * yj;
|
||||
}
|
||||
}
|
||||
y
|
||||
}
|
||||
@@ -98,11 +350,24 @@ pub(crate) fn bilinear(y: &[f64], y_prime: &[f64]) -> f64 {
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
/// Factorise a dense row-major matrix, for the goldens below.
|
||||
fn dense(a: &[f64], n: usize) -> Option<Cholesky> {
|
||||
let mut b = SymmetricBuilder::new();
|
||||
for i in 0..n {
|
||||
for j in 0..n {
|
||||
if a[i * n + j] != 0.0 {
|
||||
b.add(i, j, a[i * n + j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
Cholesky::factor(b, n)
|
||||
}
|
||||
|
||||
/// `[[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 reproduces_a_known_quadratic_form() {
|
||||
let c = Cholesky::factor(vec![4.0, 1.0, 1.0, 3.0], 2).unwrap();
|
||||
let c = dense(&[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);
|
||||
}
|
||||
@@ -113,8 +378,8 @@ mod tests {
|
||||
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();
|
||||
let a = [2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0];
|
||||
let c = dense(&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;
|
||||
@@ -127,8 +392,8 @@ mod tests {
|
||||
#[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 a = [2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0];
|
||||
let c = dense(&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);
|
||||
@@ -138,15 +403,108 @@ mod tests {
|
||||
/// 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 a = [1e12, 1e12 - 1.0, 1e12 - 1.0, 1e12];
|
||||
let c = dense(&a, 2).unwrap();
|
||||
let y = c.whiten(&[1.0, -1.0]);
|
||||
assert!(bilinear(&y, &y) >= 0.0);
|
||||
}
|
||||
|
||||
/// Against an independent dense reference, on random sparse SPD matrices.
|
||||
///
|
||||
/// The goldens above are 2x2 and 3x3 — small enough that AMD does nothing
|
||||
/// and no fill-in occurs, so they cannot catch a symbolic-pass bug. This
|
||||
/// builds matrices big enough to permute and fill in, and checks every
|
||||
/// bilinear form against a textbook dense factorisation of the *same*
|
||||
/// matrix in its original order.
|
||||
#[test]
|
||||
fn agrees_with_a_dense_reference_on_random_sparse_systems() {
|
||||
/// Dense Cholesky and quadratic form, deliberately naive: this is the
|
||||
/// reference, so it must not share code with what it is checking.
|
||||
fn dense_quadratic_form(a: &[f64], n: usize, b: &[f64], c: &[f64]) -> f64 {
|
||||
let mut l = a.to_vec();
|
||||
for j in 0..n {
|
||||
let mut d = l[j * n + j];
|
||||
for k in 0..j {
|
||||
d -= l[j * n + k] * l[j * n + k];
|
||||
}
|
||||
let d = d.sqrt();
|
||||
l[j * n + j] = d;
|
||||
for i in j + 1..n {
|
||||
let mut sum = l[i * n + j];
|
||||
for k in 0..j {
|
||||
sum -= l[i * n + k] * l[j * n + k];
|
||||
}
|
||||
l[i * n + j] = sum / d;
|
||||
}
|
||||
}
|
||||
let solve = |rhs: &[f64]| -> Vec<f64> {
|
||||
let mut y = rhs.to_vec();
|
||||
for i in 0..n {
|
||||
for k in 0..i {
|
||||
y[i] -= l[i * n + k] * y[k];
|
||||
}
|
||||
y[i] /= l[i * n + i];
|
||||
}
|
||||
y
|
||||
};
|
||||
let (yb, yc) = (solve(b), solve(c));
|
||||
yb.iter().zip(&yc).map(|(x, y)| x * y).sum()
|
||||
}
|
||||
|
||||
// A cheap deterministic generator; no dependency, and reproducible.
|
||||
let mut seed = 0x2545_F491_4F6C_DD1Du64;
|
||||
let mut rand = move || {
|
||||
seed ^= seed << 13;
|
||||
seed ^= seed >> 7;
|
||||
seed ^= seed << 17;
|
||||
(seed >> 11) as f64 / (1u64 << 53) as f64
|
||||
};
|
||||
|
||||
for n in [7usize, 23, 60] {
|
||||
let mut a = vec![0.0f64; n * n];
|
||||
// A chain plus scattered long-range couplings: the shape of a
|
||||
// time-expanded joint, where a competitor's drift link can span
|
||||
// the whole matrix.
|
||||
for i in 0..n {
|
||||
a[i * n + i] = 4.0 + rand();
|
||||
if i + 1 < n {
|
||||
let v = -(0.5 + rand() * 0.5);
|
||||
a[i * n + i + 1] = v;
|
||||
a[(i + 1) * n + i] = v;
|
||||
}
|
||||
}
|
||||
for step in 0..n / 3 {
|
||||
let i = (step * 7) % n;
|
||||
let j = (step * 29 + 3) % n;
|
||||
if i != j {
|
||||
let v = -(0.1 + rand() * 0.2);
|
||||
a[i * n + j] = v;
|
||||
a[j * n + i] = v;
|
||||
// Keep it diagonally dominant, hence positive-definite.
|
||||
a[i * n + i] += 0.6;
|
||||
a[j * n + j] += 0.6;
|
||||
}
|
||||
}
|
||||
|
||||
let sparse = dense(&a, n).expect("spd");
|
||||
|
||||
for trial in 0..8 {
|
||||
let b: Vec<f64> = (0..n).map(|_| rand() * 2.0 - 1.0).collect();
|
||||
let c: Vec<f64> = (0..n).map(|_| rand() * 2.0 - 1.0).collect();
|
||||
let got = bilinear(&sparse.whiten(&b), &sparse.whiten(&c));
|
||||
let want = dense_quadratic_form(&a, n, &b, &c);
|
||||
assert!(
|
||||
(got - want).abs() <= 1e-10 * want.abs().max(1.0),
|
||||
"n={n} trial={trial}: sparse {got} vs dense {want}"
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// A permutation must not change which matrices are rejected.
|
||||
#[test]
|
||||
fn rejects_a_non_positive_definite_matrix() {
|
||||
// Singular: the second row is a multiple of the first.
|
||||
assert!(Cholesky::factor(vec![1.0, 2.0, 2.0, 4.0], 2).is_none());
|
||||
assert!(dense(&[1.0, 2.0, 2.0, 4.0], 2).is_none());
|
||||
}
|
||||
}
|
||||
|
||||
@@ -0,0 +1,178 @@
|
||||
//! What a sparse factorisation of the joint would actually buy (#52).
|
||||
//!
|
||||
//! Run explicitly:
|
||||
//!
|
||||
//! ```text
|
||||
//! cargo test --release --features approx,measure-sparsity \
|
||||
//! --test sparsity_measurement -- --ignored --nocapture
|
||||
//! ```
|
||||
//!
|
||||
//! The whole file is gated: it reaches for the joint's sparsity pattern, which
|
||||
//! is exposed only under `measure-sparsity`.
|
||||
#![cfg(feature = "measure-sparsity")]
|
||||
|
||||
use std::collections::HashSet;
|
||||
|
||||
use trueskill_tt::{ConvergenceOptions, History};
|
||||
|
||||
/// A history shaped like the issue's fixture: many slices, scored duels,
|
||||
/// competitors reappearing across slices so the drift links are long.
|
||||
fn fitted(slices: i64, duels: usize, competitors: usize) -> History<String> {
|
||||
let mut h: History<String> = History::builder()
|
||||
.key_type::<String>()
|
||||
.mu(0.0)
|
||||
.sigma(6.0)
|
||||
.beta(1.0)
|
||||
.score_sigma(2.0)
|
||||
.gamma(0.05)
|
||||
.convergence(ConvergenceOptions {
|
||||
max_iter: trueskill_tt::ITERATIONS,
|
||||
epsilon: 1e-8,
|
||||
alpha: 1.0,
|
||||
})
|
||||
.build();
|
||||
|
||||
let mut k = 0usize;
|
||||
for t in 0..slices {
|
||||
for _ in 0..duels {
|
||||
k += 1;
|
||||
h.event(t)
|
||||
.team([format!("p{}", k % competitors)])
|
||||
.team([format!("p{}", (k + 37) % competitors)])
|
||||
.scores([
|
||||
(k as f64 * 0.3).sin().abs() * 20.0,
|
||||
(k as f64 * 0.3).cos().abs() * 20.0,
|
||||
])
|
||||
.commit()
|
||||
.expect("ingests");
|
||||
}
|
||||
}
|
||||
h.converge().expect("converges");
|
||||
h
|
||||
}
|
||||
|
||||
/// Symbolic Cholesky by row-merge: returns (nnz(L), flops).
|
||||
///
|
||||
/// Fill-in is simulated directly — for each column, the set of rows below the
|
||||
/// diagonal that are nonzero — which is exact and easily checked, at the cost
|
||||
/// of being O(n * nnz(L)) rather than the linear elimination-tree method.
|
||||
fn symbolic(n: usize, adj: &[HashSet<usize>], perm_of: &[usize]) -> (usize, f64) {
|
||||
// `perm_of[old] = new`. Build the permuted lower-triangle pattern.
|
||||
let mut cols: Vec<HashSet<usize>> = vec![HashSet::new(); n];
|
||||
for (old, nbrs) in adj.iter().enumerate() {
|
||||
let i = perm_of[old];
|
||||
for &old_j in nbrs {
|
||||
let j = perm_of[old_j];
|
||||
if j < i {
|
||||
cols[j].insert(i);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
let mut nnz = 0usize;
|
||||
let mut flops = 0.0f64;
|
||||
for j in 0..n {
|
||||
// Column j's pattern is final once every earlier column has merged in.
|
||||
let rows: Vec<usize> = cols[j].iter().copied().collect();
|
||||
let c = rows.len();
|
||||
nnz += c + 1; // below-diagonal entries plus the diagonal
|
||||
// Cholesky work for this column: one outer product over its pattern.
|
||||
flops += (c as f64 + 1.0) * (c as f64 + 1.0);
|
||||
// Fill-in: every pair in column j becomes an edge in the remaining graph.
|
||||
for (a_idx, &a) in rows.iter().enumerate() {
|
||||
for &b in &rows[a_idx + 1..] {
|
||||
let (lo, hi) = if a < b { (a, b) } else { (b, a) };
|
||||
cols[lo].insert(hi);
|
||||
}
|
||||
}
|
||||
}
|
||||
(nnz, flops)
|
||||
}
|
||||
|
||||
#[test]
|
||||
#[ignore = "measurement, run explicitly"]
|
||||
fn what_sparsity_would_buy() {
|
||||
for (slices, duels, competitors) in [(30, 8, 100), (76, 13, 200)] {
|
||||
let h = fitted(slices, duels, competitors);
|
||||
let (n, pattern) = h.joint_pattern_for_measurement();
|
||||
|
||||
let nnz_a: usize = pattern.iter().map(HashSet::len).sum::<usize>() + n;
|
||||
let dense_flops = (n as f64).powi(3) / 3.0;
|
||||
|
||||
let natural: Vec<usize> = (0..n).collect();
|
||||
let (nnz_nat, flops_nat) = symbolic(n, &pattern, &natural);
|
||||
|
||||
// AMD returns `perm[new] = old`; invert it.
|
||||
let (col_ptr, row_idx) = csc(n, &pattern);
|
||||
let p = feral_amd::amd_order(
|
||||
&feral_amd::CscPattern::new(n, &col_ptr, &row_idx).expect("valid pattern"),
|
||||
)
|
||||
.expect("amd");
|
||||
let mut perm_of = vec![0usize; n];
|
||||
for (new, &old) in p.iter().enumerate() {
|
||||
perm_of[old as usize] = new;
|
||||
}
|
||||
let (nnz_amd, flops_amd) = symbolic(n, &pattern, &perm_of);
|
||||
|
||||
println!(
|
||||
"\n=== {slices} slices x {duels} duels, {competitors} competitors ===\n\
|
||||
n = {n}\n\
|
||||
nnz(A) = {nnz_a} ({:.4}% dense)\n\
|
||||
dense flops = {:.3e}\n\
|
||||
nnz(L) natural = {nnz_nat} flops = {:.3e} ({:.1}x vs dense)\n\
|
||||
nnz(L) AMD = {nnz_amd} flops = {:.3e} ({:.1}x vs dense)",
|
||||
100.0 * nnz_a as f64 / (n * n) as f64,
|
||||
dense_flops,
|
||||
flops_nat,
|
||||
dense_flops / flops_nat,
|
||||
flops_amd,
|
||||
dense_flops / flops_amd,
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
/// Full symmetric pattern to CSC, as `feral-amd` wants it.
|
||||
fn csc(n: usize, adj: &[HashSet<usize>]) -> (Vec<i32>, Vec<i32>) {
|
||||
let mut col_ptr = Vec::with_capacity(n + 1);
|
||||
let mut row_idx = Vec::new();
|
||||
col_ptr.push(0i32);
|
||||
for (j, nbrs) in adj.iter().enumerate() {
|
||||
let mut rows: Vec<i32> = nbrs.iter().map(|&i| i as i32).collect();
|
||||
rows.push(j as i32);
|
||||
rows.sort_unstable();
|
||||
rows.dedup();
|
||||
row_idx.extend_from_slice(&rows);
|
||||
col_ptr.push(row_idx.len() as i32);
|
||||
}
|
||||
(col_ptr, row_idx)
|
||||
}
|
||||
|
||||
/// End-to-end factorisation time at the scale #52 was opened about.
|
||||
#[test]
|
||||
#[ignore = "measurement, run explicitly"]
|
||||
fn factorisation_time_at_scale() {
|
||||
use std::time::Instant;
|
||||
|
||||
for (slices, duels, competitors) in [(30, 8, 100), (76, 13, 200), (150, 26, 400)] {
|
||||
let h = fitted(slices, duels, competitors);
|
||||
let (n, _) = h.joint_pattern_for_measurement();
|
||||
|
||||
// Warm, then time.
|
||||
let _ = h.joint().expect("scored history");
|
||||
let t = Instant::now();
|
||||
let joint = h.joint().expect("scored history");
|
||||
let factor = t.elapsed();
|
||||
|
||||
let a = "p0".to_string();
|
||||
let b = "p1".to_string();
|
||||
let t = Instant::now();
|
||||
let _ = joint.posterior_of(&[(&a, 1.0), (&b, -1.0)]).expect("known");
|
||||
let query = t.elapsed();
|
||||
|
||||
println!(
|
||||
"n = {n:5} factorise = {factor:>12?} query = {query:>10?} \
|
||||
(dense was O(n^3): {:.3e} flops)",
|
||||
(n as f64).powi(3) / 3.0
|
||||
);
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user