Build a phylogenetic graph representation from a tree

Computes a symmetric-normalised Gaussian kernel adjacency matrix and Laplacian spectral node features from the cophenetic distance matrix. Results are optionally cached to an .rds file.

Usage

build_phylo_graph(tree, k_eigen = "auto", sigma_mult = 0.5, cache_path = NULL)

Arguments

tree

object of class "phylo".

k_eigen

integer or "auto". Number of Laplacian eigenvectors to use as node features. When "auto" (default), scales with tree size using min(max(ceiling(n/20), 4), 32), then clamps the result to at most max(n - 1, 1) for the tree size. For example, the formula gives 5 for 100 tips, 8 for 160, 15 for 300, and 32 for 640 or more. For a one-tip tree the clamp yields one, but the eigendecomposition rank guard still stops because no non-trivial eigenvector is available. An explicitly supplied integer is not reduced by this automatic clamp.

sigma_mult

numeric. Bandwidth multiplier (call it \(s\)): \(\sigma = \mathrm{median}(D) \times s\) (default 0.5).

cache_path

character or NULL. Path to an .rds cache file. If the file exists and dimensions match, it is loaded instead of recomputing; matching dimensions do not establish that the tree tip identities/order or bandwidth are the same. Use a distinct cache path for each unchanged tree, tip order, and sigma_mult combination, and remove or refresh the cache after changing any of them. NULL disables caching (default).

Value

A list with:

adj

Numeric matrix (n x n). Symmetric-normalised adjacency.

coords

Numeric matrix (n x k_eigen). Spectral node features.

n

Integer. Number of tips.

sigma

Numeric. Bandwidth used.

D

Numeric matrix (n x n). Cophenetic (patristic) distances between tips, row/column-ordered by tree$tip.label. Returned so downstream functions can reuse it instead of calling ape::cophenetic.phylo() again.

D_sq

Numeric matrix (n x n). Element-wise squared cophenetic distances (D * D). Used to add a learnable per-head Gaussian distance bias to attention scores. During impute() training, D_sq remains allocated while D and R_phy are cleared before model fitting.

R_phy

Numeric matrix (n x n). Phylogenetic correlation matrix cov2cor(ape::vcv(tree)) (diagonal = 1). Used by fit_baseline for BM conditional imputation.

Details

The graph is built in three steps: (1) cophenetic distances between all pairs of tips; (2) Gaussian kernel \(A_{ij} = \exp(-d_{ij}^2 / (2\sigma^2))\) with \(\sigma = \mathrm{median}(D) \times s\) (where \(s\) is the user-supplied sigma_mult); (3) symmetric normalisation \(\tilde{A} = D^{-1/2} A D^{-1/2}\) with self-loops added before normalisation.

Spectral node features are the k_eigen smallest non-zero eigenvectors of the unnormalised Laplacian, which encode the broad cluster structure of the phylogeny.

Scaling

The implementation stores the dense cophenetic-distance, adjacency, and Laplacian matrices, so memory use remains \(O(n^2)\). For trees with more than 7500 tips, RSpectra's iterative Lanczos eigensolver (eigs_sym) requests a partial spectrum from the dense Laplacian; it reduces the eigen-decomposition work, but does not make graph storage sparse. The dense base::eigen() path is used below this threshold and as a fallback when RSpectra is unavailable or Lanczos fails. Runtime depends on tree structure, hardware, and convergence. A benchmark on one 9,993-tip AVONET tree is a fixture-specific result, not a general timing guarantee. The cophenetic distance matrix is computed once per call and reused for both the adjacency and the spectral features; it is also returned in the result so that downstream consumers (notably fit_baseline) can skip recomputation.

Examples

set.seed(1)
tree <- ape::rtree(30)
g <- build_phylo_graph(tree, k_eigen = 4)
dim(g$adj)    # 30 x 30
#> [1] 30 30
dim(g$coords) # 30 x 4
#> [1] 30  4
dim(g$D)      # 30 x 30
#> [1] 30 30