Covariance Kernels and Composition

Covariance as model structure

A Gaussian-process kernel specifies how values of the latent function co-vary.

For any finite input set X, the matrix

K_{ij} = k(x_i,x_j)

must be positive semidefinite.

This package keeps the numerical kernel functions and reusable kernel specifications separate.

For example:

x <- seq(-2, 2, length.out = 5)

kernel_rbf(
  x,
  variance = 2,
  length_scale = 0.8
)
#>              [,1]        [,2]       [,3]        [,4]         [,5]
#> [1,] 2.000000e+00 0.915666724 0.08787387 0.001767653 7.453306e-06
#> [2,] 9.156667e-01 2.000000000 0.91566672 0.087873867 1.767653e-03
#> [3,] 8.787387e-02 0.915666724 2.00000000 0.915666724 8.787387e-02
#> [4,] 1.767653e-03 0.087873867 0.91566672 2.000000000 9.156667e-01
#> [5,] 7.453306e-06 0.001767653 0.08787387 0.915666724 2.000000e+00

specification <- rbf_kernel(
  variance = 2,
  length_scale = 0.8
)

evaluate_kernel(specification, x)
#>              [,1]        [,2]       [,3]        [,4]         [,5]
#> [1,] 2.000000e+00 0.915666724 0.08787387 0.001767653 7.453306e-06
#> [2,] 9.156667e-01 2.000000000 0.91566672 0.087873867 1.767653e-03
#> [3,] 8.787387e-02 0.915666724 2.00000000 0.915666724 8.787387e-02
#> [4,] 1.767653e-03 0.087873867 0.91566672 2.000000000 9.156667e-01
#> [5,] 7.453306e-06 0.001767653 0.08787387 0.915666724 2.000000e+00

Squared-exponential kernel

The squared-exponential kernel is

k(x,x') = \sigma_f^2 \exp\left( -\frac{1}{2} \sum_{d=1}^D \frac{(x_d-x_d')^2}{\ell_d^2} \right).

A single \ell gives an isotropic kernel. A vector (\ell_1,\ldots,\ell_D) gives ARD.

In code:

rbf_kernel(
  variance = 1.5,
  length_scale = 0.7
)
#> RBF(variance=1.5, length_scale=0.7)

Matérn family

The package supplies \nu=1/2, 3/2, and 5/2.

Let the scaled distance be

r = \sqrt{ \sum_{d=1}^D \frac{(x_d-x_d')^2}{\ell_d^2} }.

Then

k_{1/2}(r)= \sigma_f^2 e^{-r},

k_{3/2}(r)= \sigma_f^2(1+\sqrt{3}r)e^{-\sqrt{3}r},

and

k_{5/2}(r)= \sigma_f^2 \left( 1+\sqrt{5}r+\frac{5r^2}{3} \right) e^{-\sqrt{5}r}.

The smoothness assumption is therefore visible in the kernel choice rather than hidden inside the fitting function.

Rational quadratic

The rational-quadratic kernel is

k(r) = \sigma_f^2 \left( 1+\frac{r^2}{2\alpha} \right)^{-\alpha}.

It can be interpreted as a scale mixture of squared-exponential kernels.

rational_quadratic_kernel(
  variance = 1,
  length_scale = 0.5,
  alpha = 1.2
)
#> RationalQuadratic(variance=1, length_scale=0.5, alpha=1.2)

Periodic structure

For one input dimension,

k(x,x') = \sigma_f^2 \exp\left( -\frac{ 2\sin^2\left(\pi(x-x')/p\right) }{\ell^2} \right).

period is p, while length_scale controls smoothness within the cycle.

periodic_kernel(
  variance = 1,
  length_scale = 0.8,
  period = 12
)
#> Periodic(variance=1, length_scale=0.8, period=12)

Sums and products

If k_1 and k_2 are valid kernels, then so are

k(x,x') = k_1(x,x') + k_2(x,x')

and

k(x,x') = k_1(x,x')k_2(x,x').

The package mirrors this directly:

trend <- linear_kernel(variance = 0.05)

seasonal <- periodic_kernel(
  variance = 1,
  length_scale = 0.8,
  period = 12
)

local <- matern32_kernel(
  variance = 0.3,
  length_scale = 2
)

combined <- sum_kernel(
  trend,
  product_kernel(seasonal, local)
)

combined
#> SumKernel(
#>   Linear(variance=0.05)
#>   ProductKernel(
#>     Periodic(variance=1, length_scale=0.8, period=12)
#>     Matern-3/2(variance=0.3, length_scale=2)
#>   )
#> )

The corresponding covariance matrix is exactly the element-wise algebra:

time <- 0:6

K_combined <- evaluate_kernel(combined, time)

K_manual <-
  evaluate_kernel(trend, time) +
  evaluate_kernel(seasonal, time) *
  evaluate_kernel(local, time)

max(abs(K_combined - K_manual))
#> [1] 0

Parameter paths

Composite kernels retain explicit parameter paths:

kernel_parameters(
  combined,
  flatten = TRUE
)
#>             kernel1.variance     kernel2.kernel1.variance 
#>                         0.05                         1.00 
#> kernel2.kernel1.length_scale       kernel2.kernel1.period 
#>                         0.80                        12.00 
#>     kernel2.kernel2.variance kernel2.kernel2.length_scale 
#>                         0.30                         2.00

These paths are used directly by the hyperparameter optimizer.

Partial updates preserve the rest of the structure:

updated <- update_kernel_parameters(
  combined,
  c("kernel2.kernel1.period" = 24)
)

kernel_parameters(
  updated,
  flatten = TRUE
)
#>             kernel1.variance     kernel2.kernel1.variance 
#>                         0.05                         1.00 
#> kernel2.kernel1.length_scale       kernel2.kernel1.period 
#>                         0.80                        24.00 
#>     kernel2.kernel2.variance kernel2.kernel2.length_scale 
#>                         0.30                         2.00

Kernel composition is therefore both mathematical and operational: the expression defining the covariance remains inspectable after fitting.