Module gp 

Gaussian process regression with a constant mean and a stationary kernel, its hyperparameters fitted by maximum marginal likelihood: the surrogate model of Bayesian optimization.

Unstable for one release. This module is public from genoxide 0.13 so that a Gaussian process can be fitted and queried on its own, but its API and the bits of its fits may still change in 0.14, as the surrogate-assisted methods planned after it use it.

A GaussianProcess models a function f of the genes of a Real genome as f(x) ~ GP(m, σ_f² k(x, x′)) (Rasmussen and Williams, 2006, eq. 2.37), observed with independent normal noise of variance σ_n² (eq. 2.20), and gives the posterior mean and variance of f at any point (eq. 2.25 and 2.26), with their gradients. Built with GaussianProcess::builder, fitted with fit:

use genoxide::model::gp::GaussianProcess;
use genoxide::prelude::*;

// a smooth function of one gene, from 8 evaluations
let f = |x: f64| x.sin() + 0.1 * x * x;
let points: Vec<Reals> = (0..8).map(|i| Reals::from(vec![f64::from(i)])).collect();
let values: Vec<f64> = points.iter().map(|x| f(x[0])).collect();
let gp = GaussianProcess::builder(Real::uniform(1, 0.0..=7.0)?).fit(&points, &values)?;
// close to f between the points, and sure of itself there
let prediction = gp.predict(&[3.5]);
assert!((prediction.mean() - f(3.5)).abs() < 0.05);
assert!(prediction.sd() < 0.1);

The model

• Inputs are scaled to the unit cube by the bounds of the Real genome: each gene’s range becomes [0, 1], so the length scales of every gene start alike. A gene whose bounds are equal takes no part in the model.
• Outputs are standardized: the model fits (y − ȳ) / s, with ȳ and s the mean and the standard deviation of the values, and its predictions are mapped back. With given hyperparameters, that changes nothing but the rounding: the predictions are those of the model in the values’ own units.
• The kernel (Kernel) is stationary with a length scale per gene (automatic relevance determination, Rasmussen and Williams, eq. 5.1 and 5.2): k(x, x′) = k(r) of the scaled distance r² = Σᵢ (xᵢ − x′ᵢ)² / ℓᵢ². The Matérn kernel with ν = 5/2 by default (eq. 4.17), twice differentiable, as Snoek, Larochelle and Adams (2012) advise for Bayesian optimization against the squared exponential’s infinitely smooth functions (eq. 4.9; Stein, 1999).
• The noise (Noise) is none by default: the model interpolates the values, as suits the deterministic functions genoxide optimizes, with the jitter below as the only nugget. For a noisy function, it’s learned with the other hyperparameters, or fixed.
• The constant mean m is, for given kernel hyperparameters, the one that maximizes the marginal likelihood: the generalized least squares estimate m = 1ᵀK⁻¹y / 1ᵀK⁻¹1 (setting the derivative of eq. 2.30 with y − m for y to zero, as eq. 2.38 models a fixed mean). The hyperparameters therefore maximize the likelihood over the mean too, and the mean’s derivative being zero there, eq. 5.9 is the gradient of the likelihood so maximized.

Fitting

fit maximizes the log marginal likelihood (eq. 2.30 and 5.8) ln p(y | X, θ) = −½ (y − m)ᵀ K_y⁻¹ (y − m) − ½ ln |K_y| − (n/2) ln 2π, K_y = σ_f² K + σ_n² I, over the logarithms of the length scales, of σ_f² and of σ_n², with genoxide’s Lbfgsb and the analytic gradient of eq. 5.9, ∂/∂θⱼ ln p = ½ tr((ααᵀ − K_y⁻¹) ∂K_y/∂θⱼ), α = K_y⁻¹(y − m). The search runs from several starts: the first from fixed values (length scales of 0.5 of each range, σ_f² the values’ variance, a learned σ_n² 1e-4 of it or its least), the others from random points of the box below, drawn from streams derived from the seed, independent of each other. The likelihood’s best wins, the earlier start on ties, so the fit is the same on any number of threads (with the parallel feature, the starts run on rayon). The box, in the scaled units: length scales from 0.01 to 100, σ_f² from 1e-3 to 1e3, σ_n² from the noise’s least to 1.

Every matrix is factored by Cholesky (Rasmussen and Williams, Algorithm 2.1). A kernel matrix that rounding makes indefinite, e.g. with points very close together and little noise, gets a jitter on its diagonal: none, then 1e-10 of the diagonal’s scale, raised tenfold until it factors (GaussianProcess::jitter).

Every operation is a sum, product, quotient or square root in a fixed order, with math’s exp and ln, so a fit gives the same bits on every platform.

References: Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press, ch. 2, 4 and 5. Snoek, J., Larochelle, H. and Adams, R. P. (2012). Practical Bayesian optimization of machine learning algorithms. NeurIPS 25, arXiv:1206.2944. Stein, M. L. (1999). Interpolation of Spatial Data. Springer.

Structs

GaussianProcess

A Gaussian process fitted to evaluations of a function: its posterior mean and variance at any point, with their gradients. See the module for the model and how it’s fitted.

GaussianProcessBuilder

A builder for a GaussianProcess, from GaussianProcess::builder.

Hyperparameters

The hyperparameters of a GaussianProcess, in the units of the genes and of the values.

Prediction

A GaussianProcess’s prediction at a point: the posterior mean and variance of f there (Rasmussen and Williams, 2006, eq. 2.25 and 2.26), without the noise.

Enums

Kernel

The kernel of a GaussianProcess: the correlation of f at two points as a function of their scaled distance r, with r² = Σᵢ (xᵢ − x′ᵢ)² / ℓᵢ², a length scale ℓᵢ per gene.

Noise

The observation noise of a GaussianProcess: its variance σ_n², as a fraction of the values’ variance (the model’s outputs are standardized).