TreeKernelSum

Fast approximate all-pairs kernel summation on weighted trees via centroid decomposition and sum-of-exponentials approximation.

Given a weighted tree T on n nodes with node values mu and a kernel K, computes

Phi(v) = sum_u K(d_T(v,u)) * mu(u)   for all v

in O(r n log n) time instead of the naive O(n^2), where r is the number of exponential terms in a sum-of-exponentials (SOE) approximation of K.

Key results

• Single exponential (OU kernel): exact to machine epsilon, 972x speedup at n=50k
• Gaussian kernel (r=13 terms): 495x speedup at n=50k, scales to 1M nodes in 12s
• Real phylogenies: machine-precision OU kernel sums on Fritz (5,020 species) and Upham (5,912 species) mammal trees
• Kernel prediction: Nadaraya-Watson LOO cross-validation of mammalian body mass across four kernel families (OU, Matern-1.5, Matern-2.5, Gaussian), R^2 = 0.955

Prerequisites

• Rust (2021 edition, stable toolchain)
• No other dependencies required; all crates are fetched by Cargo

For the Python prototype (optional):

• Python 3.8+ with NumPy and SciPy

Build

cargo build --release

Run

cargo run --release

This runs the full validation and benchmark suite in sequence:

1. Correctness validation — single-exponential exactness, general kernel error bounds, edge cases
2. SOE comparison — empirical term count vs Braess-Hackbusch theoretical bounds for four kernel types
3. Performance benchmarks — naive vs TreeKernelSum at n = 1k to 1M for single-exponential and Gaussian kernels
4. Phylogenetics demo — OU kernel summation on simulated coalescent trees (1k-500k tips)
5. Real tree benchmarks — OU kernel on Fritz and Upham mammal phylogenies with comparison to published methods
6. Kernel prediction — body mass prediction on Fritz tree with OU, Matern-1.5, Matern-2.5, and Gaussian kernels

The full run takes approximately 5-10 minutes. The correctness validation sections print PASS/FAIL for each test. All benchmarks use deterministic ChaCha8 RNG with fixed seeds for exact reproducibility.

Project structure

src/
  main.rs                 # Entry point: runs all validation and benchmarks
  tree.rs                 # WeightedTree, random tree generator (Prufer sequence),
                          #   path/star trees, all-pairs distances
  kernels.rs              # Kernel definitions (Gaussian, exponential, Matern, InvMQ),
                          #   naive O(n^2) kernel sum
  centroid.rs             # CentroidDecomposition::build, tree_kernel_sum (core algorithm),
                          #   single_exp_sum (r=1 wrapper)
  soe.rs                  # SOE fitting: dictionary OMP (legacy) + sinc-quadrature
                          #   (preferred). auto_fit_sinc for automatic fitting
  soe_comparison.rs       # Formal comparison: OMP vs sinc vs Braess-Hackbusch bounds
  bench.rs                # validate_correctness + run_benchmarks
  phylo.rs                # Phylogenetics OU demo on coalescent trees
  newick.rs               # Newick/Nexus tree file parsers
  real_tree_bench.rs      # Benchmarks on published mammal phylogenies
  prediction.rs           # Kernel prediction of mammal body mass (Nadaraya-Watson LOO)

data/
  Fritz_mammal_tree.nwk           # Fritz et al. (2009) mammal supertree, 5020 tips
  Upham_mammal_5911sp_MCC.nex     # Upham et al. (2019) MCC mammal tree, 5912 tips
  PanTHERIA_WR05.txt              # PanTHERIA mammal trait database (Jones et al. 2009)

py_prototype/             # Archived Python prototype (NumPy/SciPy)
  tree.py                 #   Same structure as Rust, used for initial development
  kernels.py
  centroid.py
  test_correctness.py

paper/
  paper.tex               # Manuscript draft (LaTeX, article class)

TreeKernelSum.md          # Full theory: problem statement, algorithm, proofs
RESULTS.md                # All empirical results and benchmark data

The algorithm

1. Phase 0 (SOE fit): Approximate K(d) as a sum of r exponentials on [0, D].
2. Phase 1 (Centroid decomposition): Build centroid hierarchy. Each node stores O(log n) ancestor centroids with distances.
3. Phase 2 (Aggregation): For each exponential term, compute weighted sums at each centroid and per-child partial sums.
4. Phase 3 (Assembly): For each node, walk centroid ancestors and accumulate contributions with inclusion-exclusion.

Total: O(r n log n) time, O(n log n) space. For standard kernels, r = O(log(D/eps) log(1/eps)).

Validated claims

Data sources

• Fritz tree: Fritz et al. (2009) "Geographical variation in predictors of mammalian extinction risk." Ecology Letters 12:538-549.
• Upham tree: Upham et al. (2019) "Inferring the mammal tree." PLoS Biology 17:e3000494.
• PanTHERIA: Jones et al. (2009) "PanTHERIA: a species-level database of life history, ecology, and geography of extant and recently extinct mammals." Ecology 90:2648.

License

MIT