Source code for gpflow.kernels.misc

# Copyright 2018-2020 The GPflow Contributors. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from typing import Optional

import numpy as np
import tensorflow as tf

from ..base import Parameter, TensorType
from ..utilities import positive, to_default_float
from .base import ActiveDims, Kernel


[docs]class ArcCosine(Kernel): """ The Arc-cosine family of kernels which mimics the computation in neural networks. The order parameter specifies the assumed activation function. The Multi Layer Perceptron (MLP) kernel is closely related to the ArcCosine kernel of order 0. The key reference is :: @incollection{NIPS2009_3628, title = {Kernel Methods for Deep Learning}, author = {Youngmin Cho and Lawrence K. Saul}, booktitle = {Advances in Neural Information Processing Systems 22}, year = {2009}, url = {http://papers.nips.cc/paper/3628-kernel-methods-for-deep-learning.pdf} } """ implemented_orders = {0, 1, 2} def __init__( self, order: int = 0, variance: TensorType = 1.0, weight_variances: TensorType = 1.0, bias_variance: TensorType = 1.0, *, active_dims: Optional[ActiveDims] = None, name: Optional[str] = None, ) -> None: """ :param order: specifies the activation function of the neural network the function is a rectified monomial of the chosen order :param variance: the (initial) value for the variance parameter :param weight_variances: the (initial) value for the weight_variances parameter, to induce ARD behaviour this must be initialised as an array the same length as the the number of active dimensions e.g. [1., 1., 1.] :param bias_variance: the (initial) value for the bias_variance parameter defaults to 1.0 :param active_dims: a slice or list specifying which columns of X are used """ super().__init__(active_dims=active_dims, name=name) if order not in self.implemented_orders: raise ValueError("Requested kernel order is not implemented.") self.order = order self.variance = Parameter(variance, transform=positive()) self.bias_variance = Parameter(bias_variance, transform=positive()) self.weight_variances = Parameter(weight_variances, transform=positive()) self._validate_ard_active_dims(self.weight_variances) @property def ard(self) -> bool: """ Whether ARD behaviour is active. """ return self.weight_variances.shape.ndims > 0 def _weighted_product(self, X: TensorType, X2: Optional[TensorType] = None) -> tf.Tensor: if X2 is None: return tf.reduce_sum(self.weight_variances * tf.square(X), axis=1) + self.bias_variance return ( tf.linalg.matmul((self.weight_variances * X), X2, transpose_b=True) + self.bias_variance ) def _J(self, theta: TensorType) -> TensorType: """ Implements the order dependent family of functions defined in equations 4 to 7 in the reference paper. """ if self.order == 0: return np.pi - theta elif self.order == 1: return tf.sin(theta) + (np.pi - theta) * tf.cos(theta) else: assert self.order == 2, f"Don't know how to handle order {self.order}." return 3.0 * tf.sin(theta) * tf.cos(theta) + (np.pi - theta) * ( 1.0 + 2.0 * tf.cos(theta) ** 2 ) def K(self, X: TensorType, X2: Optional[TensorType] = None) -> tf.Tensor: X_denominator = tf.sqrt(self._weighted_product(X)) if X2 is None: X2 = X X2_denominator = X_denominator else: X2_denominator = tf.sqrt(self._weighted_product(X2)) numerator = self._weighted_product(X, X2) cos_theta = numerator / X_denominator[:, None] / X2_denominator[None, :] jitter = 1e-15 theta = tf.acos(jitter + (1 - 2 * jitter) * cos_theta) return ( self.variance * (1.0 / np.pi) * self._J(theta) * X_denominator[:, None] ** self.order * X2_denominator[None, :] ** self.order ) def K_diag(self, X: TensorType) -> tf.Tensor: X_product = self._weighted_product(X) const = (1.0 / np.pi) * self._J(to_default_float(0.0)) return self.variance * const * X_product ** self.order
[docs]class Coregion(Kernel): """ A Coregionalization kernel. The inputs to this kernel are _integers_ (we cast them from floats as needed) which usually specify the *outputs* of a Coregionalization model. The kernel function is an indexing of a positive-definite matrix: K(x, y) = B[x, y] . To ensure that B is positive-definite, it is specified by the two parameters of this kernel, W and kappa: B = W Wᵀ + diag(kappa) . We refer to the size of B as "output_dim x output_dim", since this is the number of outputs in a coregionalization model. We refer to the number of columns on W as 'rank': it is the number of degrees of correlation between the outputs. NB. There is a symmetry between the elements of W, which creates a local minimum at W=0. To avoid this, it is recommended to initialize the optimization (or MCMC chain) using a random W. """ def __init__( self, output_dim: int, rank: int, *, active_dims: Optional[ActiveDims] = None, name: Optional[str] = None, ) -> None: """ :param output_dim: number of outputs expected (0 <= X < output_dim) :param rank: number of degrees of correlation between outputs """ # assert input_dim == 1, "Coregion kernel in 1D only" super().__init__(active_dims=active_dims, name=name) self.output_dim = output_dim self.rank = rank W = 0.1 * np.ones((self.output_dim, self.rank)) kappa = np.ones(self.output_dim) self.W = Parameter(W) self.kappa = Parameter(kappa, transform=positive()) def output_covariance(self) -> tf.Tensor: B = tf.linalg.matmul(self.W, self.W, transpose_b=True) + tf.linalg.diag(self.kappa) return B def output_variance(self) -> tf.Tensor: B_diag = tf.reduce_sum(tf.square(self.W), 1) + self.kappa return B_diag def K(self, X: TensorType, X2: Optional[TensorType] = None) -> tf.Tensor: shape_constraints = [ (X, [..., "N", 1]), ] if X2 is not None: shape_constraints.append((X2, [..., "M", 1])) tf.debugging.assert_shapes(shape_constraints) X = tf.cast(X[..., 0], tf.int32) if X2 is None: X2 = X else: X2 = tf.cast(X2[..., 0], tf.int32) B = self.output_covariance() return tf.gather(tf.transpose(tf.gather(B, X2)), X) def K_diag(self, X: TensorType) -> tf.Tensor: tf.debugging.assert_shapes([(X, [..., "N", 1])]) X = tf.cast(X[..., 0], tf.int32) B_diag = self.output_variance() return tf.gather(B_diag, X)