Matérn Gaussian Processes on Graphs

In plain words

A Gaussian process is a way to guess an unknown function and say how sure the guess is. Its behaviour is set by a "kernel" that says how similar two inputs are. For inputs laid out in ordinary space, a well-known family of kernels (Matérn) covers smooth, rough and in-between functions. For inputs that are nodes of a network, such as road junctions or papers in a citation graph, there were few such kernels to choose from.

This paper takes the Matérn idea from its differential-equation description and rewrites it on a graph, using the graph's Laplacian, a matrix that records how nodes connect. The result is a covariance between nodes that falls off with graph distance and that can be tuned in the same way as the usual kernel: a length scale, a variance and a smoothness setting. The authors show how to train these models at scale with a few standard tricks, and how to use them for classification as well as regression.

The payoff is a kernel family that fits graph data as naturally as Matérn fits ordinary coordinates, and that works with existing Gaussian-process tooling.

TL;DR

• Graph Matérn GPs are defined by the Matérn stochastic-PDE form, with the graph Laplacian in place of the Euclidean Laplacian §3.
• The covariance is built from Laplacian eigenpairs. The smoothness ν is the exponent, the length scale κ shifts the eigenvalues, and the prior variance varies from vertex to vertex §3 fig. 1.
• Scalable training uses truncated eigenpairs, sparse precision matrices for integer ν, or sparse variational inference with graph inducing points §3.1.
• Author-reported: on traffic-speed interpolation and Cora classification, the graph Matérn kernel is similar to or slightly better than the graph diffusion kernel tab. 1 §4.3.

Problem

Replacing the Euclidean distance in a stationary kernel with graph distance generally does not give a valid covariance §1. Existing graph GPs are Markov models (GMRFs), which limit flexibility, or a few hand-built kernels such as the diffusion and random-walk kernels, whose choice is narrow §1. The paper assumes the graph is finite and undirected with non-negative edge weights §3.

Method

Setting. Let G=(V,E) be a weighted undirected graph with non-negative weights, adjacency W, degree D and Laplacian L = D - W. The sign is the positive-semidefinite convention, the analyst's rather than the geometer's §3, footnote. L is symmetric and positive semidefinite, and for a connected graph has a one-dimensional null space §3. Write L = U\Lambda U^\top with U orthogonal §3.

Functional calculus. For a scalar function \Phi, set \Phi(L) = U\Phi(\Lambda)U^\top §3.

Graph Matérn GP. Replace white noise in the Matérn SPDE by standard Gaussian noise and apply \Phi(\lambda) = (2\nu/\kappa^2 + \lambda)^{\nu/2}. The resulting prior is f \sim \mathcal{N}(0, \Psi(L)) with \Psi(\lambda) = (2\nu/\kappa^2 + \lambda)^{-\nu} §3. Note that the constant 2\nu/\kappa^2 is added to the eigenvalues, not to the entries of L §3. The paper writes \Phi with a plus sign because of the Laplacian convention §3, footnote.

Graph diffusion GP (limit \nu \to \infty, scaled): f \sim \mathcal{N}(0, e^{-\kappa^2 L/2}), with \Psi(\lambda) = e^{-\kappa^2\lambda/2} §3. The Matérn kernel converges to it as \nu \to \infty when scaled by (2\nu/\kappa^2)^\nu, because the eigenvectors are the same §3.

Random-walk limit. The kernel (I - (1-\alpha)D^{-1/2}LD^{-1/2})^p tends to the normalised diffusion kernel as p\to\infty with \alpha = 1 - \kappa^2/(2p) §3.

Normalised Laplacian variants. The same formulas with L replaced by D^{-1/2} L D^{-1/2} define the symmetric normalised Matérn and diffusion GPs. The choice between the two Laplacians depends on the application §3.

Sparsity. For integer \nu and sparse graphs, the precision matrix (2\nu/\kappa^2 + L)^\nu is sparse, so GMRF machinery applies §3. The paper notes that the d/4 term present in Euclidean and Riemannian cases is absent here, which changes the scaling and makes integer \nu the Markov case §3.

Output cross-covariance. For multi-output functions indexed by graph nodes, the paper uses an intrinsic coregionalisation model with a separable kernel k_{\mathbb{R}^d}(x,x')\,k_G(i,i') §3.

Scalable training (§3.1).

1. Graph Fourier features. Compute the l smallest eigenpairs of L, for example with Lanczos. Because \Psi is decreasing, these give the largest kernel eigenvalues. Treat the approximate GP as Bayesian linear regression on the eigenvectors §3.1 app. B.
2. Woodbury solve. With \tilde U(x) the rows of the retained eigenvectors at the training nodes, (K + \sigma^2 I)^{-1} \approx \sigma^{-2} I - \sigma^{-4}\tilde U(\Psi^{-1} + \sigma^{-2}\tilde U^\top\tilde U)^{-1}\tilde U^\top. The paper calls the inner matrix diagonal app. B; see Pitfalls.
3. Non-conjugate learning. Use sparse variational inference with inducing inputs z drawn from the graph and q(u) = \mathcal{N}(\mu,\Sigma). Minimise the KL divergence by stochastic optimisation on mini-batches §3.1. For classification, the expected log-likelihood reduces to cross-entropy.
4. Interdomain inducing variables. Use the kernel eigenvectors as inducing functions §3.1.

Hyperparameter learning. The kernel parameters are \kappa, \sigma^2 and \nu. Eigenvector features make the derivative with respect to \nu simple, so \nu can be optimised like \kappa app. B.

Pseudocode for the paper's kernel on l retained eigenpairs §3 app. B:

# Graph Matérn covariance from l Laplacian eigenpairs  [§3] [app. B]
# Input:  eigenvalues (l,), eigenvectors (N, l), nu, kappa
# Output: covariance (N, N)

def graph_matern_covariance(eigenvalues, eigenvectors, nu, kappa):
    # Spectral weight per eigenpair  [§3] [app. B]
    # Diffusion variant: exp(-kappa**2 / 2 * eigenvalues)  [§3]
    weights = (2 * nu / kappa**2 + eigenvalues) ** (-nu)
    # Sum over retained pairs: K = U diag(weights) U^T  [app. B]
    return (eigenvectors * weights) @ eigenvectors.T

Pseudocode for the eigen-feature solve in app. B:

# Solve (K + noise I) alpha = y with eigen features  [app. B]
# Input:  features (n, l) = U_tilde(x), weights (l,), noise, targets (n,)
# Output: alpha (n,)

def woodbury_solve(features, weights, noise, targets):
    # Inner matrix Psi^-1 + noise^-1 U^T U; the paper calls it diagonal
    # (see Pitfalls)  [app. B]
    inner = diag(1.0 / weights) + (features.T @ features) / noise
    # Projected targets  [app. B]
    projected = features.T @ targets
    # Woodbury identity  [app. B]
    correction = features @ solve(inner, projected) / noise**2
    return targets / noise - correction

Hyperparameters

Values for the two experiments app. A, with the repository's notebook values where they differ. Values "not stated" are not given in the paper.

Results

Author-reported, on two tasks §4.

• Traffic: interpolation of Monday 17:30 flow speed on a road graph, 250 training and 75 test nodes, 10 repeats. Metric: RMSE. The standard deviation of the speed across the whole graph is about 17.1 mph tab. 1.
• Cora: seven-class topic classification on the largest connected component of the citation graph, 140 training nodes, 1000 test nodes, 10 repeats. Metric: accuracy §4.2 tab. 1.

Source: tab. 1. Mean with standard deviation in parentheses. Random walk and inverse cosine are shown only with the normalised Laplacian tab. 1. The paper's conclusion is that the Matérn class gives similar or slightly better results than these alternatives §4.3.

Erratum. The published version labelled the traffic values as MSE, and they were scaled by 1/\sqrt{75}. The current version corrects both points. The scaling was applied to the proposed method and to every baseline, so relative comparisons are unaffected § Erratum.

Figures (described, not reproduced). fig. 1: prior variance on a complete graph is uniform; on a star graph the centre has lower variance; on two random graphs the variance follows the graph structure, not the individual degrees §3. fig. 2: a Matérn kernel on increasingly fine triangle meshes converges to the Riemannian kernel on the sphere §3. fig. 3: predicted mean and standard deviation of traffic speed over the graph; uncertainty grows with distance from sensors §4.1. fig. 4: predictions at a junction differ between the two opposite carriageways §4.1. fig. 5: Cora ground truth, and errors in red §4.2.

Paper vs code

Repository: spbu-math-cs/Graph-Gaussian-Processes at 781be70c8913b96c365475889f1f39629a5c9d0f, MIT, with the LICENSE naming all six authors. The package is graph_matern (kernels, inducing variables, a GraphSVGP class). The examples are two notebooks.

• Spectral normalisation (not in the paper). Both kernels rescale the weights so that their sum over the N vertices equals N\sigma_f: S = S * N / sum(S), then S * sigma_f graph_matern_kernel.py graph_diffusion_kernel.py. Because the retained eigenvectors have unit norm, the average prior variance over vertices equals \sigma_f exactly for the truncated kernel. The paper's formulas \Psi(\lambda) §3 app. B have no such factor, and §3 does not define \sigma^2 for graph kernels. The two therefore differ by a constant that depends on \nu and \kappa.
• Cora \nu. Paper initialises \nu = 3 app. A; the notebook uses \nu = 3/2 classification.ipynb.
• Traffic \kappa. Paper initialises \kappa = 3 app. A; the notebook uses \kappa = 5 regression.ipynb.
• Iterations. Paper: 20000 Adam steps for both tasks app. A. Notebooks: 12000 for traffic regression.ipynb and 1000 for Cora classification.ipynb.
• Edge weights. Paper: inverse travel distance app. A. Code: each weight is the mean edge length divided by the edge length preprocessing.py. That is the paper's weight times a constant. Scaling L also scales its eigenvalues while the 2\nu/\kappa^2 offset stays fixed, so the learned \kappa is not the paper's \kappa on the same graph.
• Traffic data. Paper: CPMS sensor data and OpenStreetMap highways in San Jose, Monday 2 January 2017 at 17:30 app. A. Code: the PEMS-Bay file pems-bay.h5 and a Bay Area bounding box with a motorway-only filter preprocessing.py preprocessing.py. The code selects all Monday 17:30 records and keeps the first rows preprocessing.py preprocessing.py. Which timestamp that is, and whether it matches the paper's date, was not checked.
• Two directions. Paper: the two travel directions of a highway are separate edges §4.1 app. A. Code: calls osmnx.get_undirected, first before sensor insertion and again after it preprocessing.py preprocessing.py. How osmnx merges the two directions was not checked, so it is unclear whether the paper's separate edges survive.
• Parallel edges. Paper: keep one edge between two nodes, weight equal to the inverse mean travel distance app. A. No such step was found in the preprocessing code.
• Likelihood. Paper: robust-max with \varepsilon = 10^{-3} app. A. The notebook passes gpflow.likelihoods.MultiClass(7) with no epsilon classification.ipynb. Any default robust-max lives in GPflow, which is not in this repository, so this was not verified.
• Cora model. The notebook uses GPflow's SVGP. The GraphSVGP class (Fourier-feature model) is commented out in the notebook classification.ipynb, and the eigen-feature inducing variables (inducing_variables=[0]*num_eigenpairs) are not used.
• Laplacian. Both notebooks build nx.laplacian_matrix(G), which is the unnormalised Laplacian regression.ipynb classification.ipynb. The normalised rows of Table 1 are therefore not produced by the notebooks as shipped.
• Eigendecomposition. Paper: for large graphs, use Lanczos to compute the l smallest eigenpairs §3.1. Code: densifies the whole Laplacian and calls tf.linalg.eigh regression.ipynb, then slices.
• Sparse and Woodbury paths. The GMRF sparse-precision path §3 and the Woodbury solve app. B have no implementation in the repository.
• Metric. The regression notebook's printed "relative l2 error" is not an RMSE; the README says multiply it by \sqrt{N_\text{test}} to get the RMSE README.md. The notebook's saved single-run output is 1.266 before that conversion regression.ipynb.
• Status. The README says the project is no longer actively maintained and points to GeometricKernels README.md.

The weight functions (2\nu/\kappa^2+\lambda)^{-\nu} and e^{-\kappa^2\lambda/2} match §3 and the code, apart from the normalisation above. The Adam learning rate matches app. A.

Pitfalls

• Normalisation changes the scale. Implementing the paper's formula directly gives a different prior variance from the repository, because the repository rescales so the average vertex variance is \sigma_f graph_matern_kernel.py. Pick one convention, and say which, before comparing with paper numbers.
• Diagonal claim in appendix B. The paper says the inner matrix \Psi^{-1} + \sigma^{-2}\tilde U(x)^\top\tilde U(x) is diagonal app. B. Orthogonality of U over all vertices makes U^\top U diagonal only when x is the full vertex set. For a training subset, \tilde U(x)^\top\tilde U(x) is in general not diagonal. Do not rely on a diagonal solve without checking.
• Dense eigendecomposition. The code builds the full dense Laplacian and calls tf.linalg.eigh regression.ipynb. Memory and time grow with N^2 and N^3. The paper's Lanczos suggestion §3.1 is what scales.
• Truncation error not quantified. The paper uses 500 eigenpairs for both tasks app. A without a stated error bound. Keeping the smallest-eigenvalue pairs keeps the largest kernel weights §3.1, but how well 500 pairs approximate the full kernel on another graph is not stated.
• Null space and connectivity. The one-dimensional null space holds for a connected graph §3. The Cora code keeps only the largest component preprocessing.py. The traffic graph was not checked for connectivity.
• Sign convention. The paper's L is positive semidefinite §3. If you use the negative (geometer's) Laplacian, the spectral weights must be written with the opposite sign on the eigenvalues §3, footnote.
• Learned \nu varies. Across runs the learned \nu varies strongly: traffic mean 1.8 (sd 0.9) and Cora mean 3.2 (sd 2.7) app. A. Results from a single run may not carry over.
• Non-uniform prior variance. The prior variance differs per vertex and is not a function of degree alone §3 fig. 1. Do not assume a constant prior variance over the graph.
• Laplacian choice matters. The normalised and unnormalised Laplacians give different traffic RMSE (for example 11.86 vs 12.82 for Matérn) tab. 1. The paper says the choice is application-dependent §3.
• Mini-batch and whitening settings. For Cora the batch is the whole training set, the covariance is diagonal and the inducing covariance is whitened app. A classification.ipynb classification.ipynb. The paper says the diagonal form worked best in its setting app. A. Other settings were not compared.
• Erratum. The traffic values in tab. 1 are the corrected RMSE §Erratum. Earlier tables that report MSE or the 1/\sqrt{75} scaling are wrong.

Lineage

• Builds on: the Matérn SPDE characterisation in Euclidean and Riemannian settings §2.1 §3, the GMRF framework §3, graph random-walk and diffusion kernels §1 §3, and inducing-point variational inference §3.1.
• Compared against: graph diffusion kernel, random-walk kernel and inverse-cosine kernel §4.3 tab. 1.
• Superseded by: not stated in the paper. The repository README recommends GeometricKernels README.md.