moreo
moreo

Reputation: 81

How to write a custom kernel in GPflow for the covariance matrix RBF plus noise only on the main block diagonal?

The required covariance matrix is

enter image description here

where t is 1D time and k={0, 1}

A sample from the kernel should look like:

enter image description here

with the orange sequence corresponding to k=0, and the blue one to k=1.

Upvotes: 0

Views: 463

Answers (2)

STJ
STJ

Reputation: 1537

It sounds like you're looking for a kernel for different discrete outputs. You can achieve this in GPflow for example with the Coregion kernel, for which there is a tutorial notebook.

To construct a coregionalization kernel that is block-diagonal (all off-diagonal entries are zero), you can set rank=0. Note that you need to explicitly specify which kernel should act on which dimensions:

import gpflow
k_time = gpflow.kernels.SquaredExponential(active_dims=[0])
k_coreg = gpflow.kernels.Coregion(output_dim=2, rank=0, active_dims=[1])

You can combine them with * as in the notebook, or with + as specified in the question:

k = k_time + k_coreg

You can see that the k_coreg term is block-diagonal as you specified: Evaluating

test_inputs = np.array([
    [0.1, 0.0],
    [0.5, 0.0],
    [0.7, 1.0],
    [0.1, 1.0],
])
k_coreg(test_inputs)

returns

<tf.Tensor: shape=(4, 4), dtype=float64, numpy=
array([[1., 1., 0., 0.],
       [1., 1., 0., 0.],
       [0., 0., 1., 1.],
       [0., 0., 1., 1.]])>

And you can get samples as in the graph in the question by running

import numpy as np
num_inputs = 51
num_outputs = 2
X = np.linspace(0, 5, num_inputs)
Q = np.arange(num_outputs)
XX, QQ = np.meshgrid(X, Q, indexing='ij')
pts = np.c_[XX.flatten(), QQ.flatten()]

K = k(pts)
L = np.linalg.cholesky(K + 1e-8 * np.eye(len(K)))

num_samples = 3
v = np.random.randn(len(L), num_samples)
f = L @ v

import matplotlib.pyplot as plt
for i in range(num_samples):
    plt.plot(X, f[:, i].reshape(num_inputs, num_outputs))

Upvotes: 1

Vincent Dutordoir
Vincent Dutordoir

Reputation: 276

In GPflow, you can construct this kernel using a sum kernel consisting of a Squared Exponential (RBF) and a White kernel.

import gpflow
kernel = gpflow.kernels.SquaredExponential() + gpflow.kernels.White()

Upvotes: 0

Related Questions