Create a coregionalization kernel for several outputs

A coregionalization kernel acts on one input column that holds an output index $p \in \{1, \ldots, P\}$ and gives the covariance between outputs, $k(p, q) = B_{pq}$ with $$B = W W^\top + \operatorname{diag}(\kappa),$$ which is positive semidefinite by construction. Multiplied by a kernel on the other inputs it gives the intrinsic coregionalization model (ICM), $\operatorname{cov}(f_p(x), f_q(x')) = B_{pq}\, k(x, x')$; a sum of such products is the linear model of coregionalization (LMC).

Usage

coregionalization_kernel(n_outputs, rank = 1L, W = NULL, kappa = NULL)

Arguments

n_outputs: Number of outputs $P$.
rank: Number of columns $R$ of $W$: the rank of the shared part of $B$.
W: $P \times R$ matrix, or a vector of its $P R$ values in column-major order. The default, zero, makes the outputs independent.
kappa: Positive output-specific variances, one per output; the default is 1.

Value

A gaussianprocesses_kernel.

Details

The kernel needs exactly one input column, so in a model with other inputs it is restricted to the output-index column with select_dimensions(). With inputs in column 1 and the output index in column 2, an ICM kernel is product_kernel(select_dimensions(rbf_kernel(), 1), select_dimensions(coregionalization_kernel(P), 2)). gp_stack_outputs() builds such inputs from one response column per output, and initialize_coregionalization() gives a starting $W$ and $\kappa$.

The optimizer works on $W$ directly and on $\log \kappa$. Two cautions apply:

$W$ is identified only up to rotation and sign: $W Q$ gives the same $B$ for any orthogonal $Q$. Report $B$, with coregionalization_matrix(), not $W$.
$W = 0$ is a stationary point of the likelihood, where the gradient with respect to $W$ is zero, so an optimization started there keeps $W = 0$; start from initialize_coregionalization() instead.

The overall scale of an ICM term is shared by $B$ and the variance of the input kernel; fixing the latter at 1 removes the redundancy only if that parameter is not optimized, so expect the two to trade off.

Stability

Experimental: this interface may change in a minor release, and every change is listed in NEWS. See gaussianprocesses-package for the policy.

References

Bonilla, E. V., Chai, K. M. A., and Williams, C. K. I. (2008). Multi-task Gaussian process prediction. Advances in Neural Information Processing Systems, 20.

Alvarez, M. A., Rosasco, L., and Lawrence, N. D. (2012). Kernels for vector-valued functions: a review. Foundations and Trends in Machine Learning, 4(3), 195–266.

Examples

kernel <- coregionalization_kernel(
  n_outputs = 2,
  rank = 1,
  W = c(1, 0.5),
  kappa = c(0.1, 0.2)
)
coregionalization_matrix(kernel)
#> [[1]]
#>      [,1] [,2]
#> [1,]  1.1 0.50
#> [2,]  0.5 0.45
#> 
evaluate_kernel(kernel, matrix(c(1, 2, 2)))
#>      [,1] [,2] [,3]
#> [1,]  1.1 0.50 0.50
#> [2,]  0.5 0.45 0.45
#> [3,]  0.5 0.45 0.45

# An ICM kernel for inputs in column 1 and the output index in column 2.
icm <- product_kernel(
  select_dimensions(rbf_kernel(length_scale = 0.5), 1),
  select_dimensions(kernel, 2)
)