Regression using a piecewise kernel¶

This notebook explains how to use Piecewise kernels in Markovflow models.

We focus on the Sparse Variational Gaussian Process model.

Our probabilistic model for this data is:

$\begin{split}\begin{align} f \sim \mathcal{GP}(0, k(., .)) \\ y_i \sim f(x_i) + \mathcal{N}(0, \epsilon^2) \end{align}\end{split}$

NOTE: If you have difficulty running this notebook, consider clearing the output and then restarting the kernel.

[1]:

%load_ext autoreload
%autoreload 2

[2]:

# Setup

import warnings

# Turn off warnings
warnings.simplefilter('ignore')

from typing import Tuple

import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
import gpflow

from gpflow import default_float, set_trainable
from gpflow.ci_utils import ci_niter
from gpflow.likelihoods import Gaussian

from markovflow.models.sparse_variational import SparseVariationalGaussianProcess
from markovflow.kernels import Matern52
from markovflow.kernels import PiecewiseKernel
from markovflow.ssm_natgrad import SSMNaturalGradient
FLOAT_TYPE = default_float()

# uncomment in notebook
# try:
#     from IPython import get_ipython
#     get_ipython().run_line_magic('matplotlib', 'inline')
# except AttributeError:
#     print('Magic function can only be used in IPython environment')
#     matplotlib.use('Agg')

plt.rcParams["figure.figsize"] = [15, 8]

2022-09-17 15:51:02.054012: W tensorflow/stream_executor/platform/default/dso_loader.cc:60] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /opt/hostedtoolcache/Python/3.7.13/x64/lib
2022-09-17 15:51:02.054049: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.

Step 1: Generate training data¶

First, let’s generate a frequency modulated wave form.

[3]:

def create_observations(time_points: np.ndarray) -> Tuple[tf.Tensor, tf.Tensor]:
    """
    A helper function to create training data.
    :param time_points: Time points to generate observations for.
    :return: Tuple[x,y] Data that represents the observations' shapes:
        X = [num_points, 1],
        Y = [num_points, state_dim , 1] where state_dim is currently 1
    """
    omega_ = np.exp((time_points - 50.) / 6.) / (1. + np.exp((time_points - 50.) / 6.)) + 0.1
    y = (np.cos(time_points * omega_ / 3) * np.sin(time_points * omega_ / 3)).reshape(-1, 1)
    return time_points, y + np.random.randn(*y.shape) * .1

# Generate some observations
N = 300
time_points, observations = create_observations(np.linspace(0, 100, N))

# Plot
plt.figure(figsize=(10, 6))
plt.plot(time_points, observations, 'C0x', ms=8, mew=2)
plt.xlabel("Time")
plt.ylabel("Label")
plt.show()
plt.close()

../_images/notebooks_piecewise_kernels_5_0.png

Step 2: Build a Piecewise kernel¶

[4]:

num_inducing = 30
step_z = N // num_inducing
num_change_points = 5
step_c = num_inducing // num_change_points

# What happens if you don't choose your inducing points from your data
inducing_points = time_points[::step_z]
inducing_points = np.linspace(np.min(time_points), np.max(time_points), num_inducing)

# What happens if you don't choose your change points from your inducing points
change_points = inducing_points[::step_c]
change_points = np.linspace(np.min(time_points), np.max(time_points), num_change_points)

assert num_change_points == len(change_points)

base = Matern52
state_dim = 3
variances = np.array([1.] * (num_change_points + 1))
lengthscales = np.array([4.] * (num_change_points + 1))

ks = [base(variance=variances[l],
            lengthscale=lengthscales[l])
      for l in range(num_change_points + 1)]

# Set state mean trainable
[k.set_state_mean(tf.zeros(state_dim,), trainable=True) for k in ks]

kernel = PiecewiseKernel(
    ks, tf.convert_to_tensor(change_points, dtype=default_float()))

2022-09-17 15:51:03.764431: I tensorflow/compiler/jit/xla_cpu_device.cc:41] Not creating XLA devices, tf_xla_enable_xla_devices not set
2022-09-17 15:51:03.764627: W tensorflow/stream_executor/platform/default/dso_loader.cc:60] Could not load dynamic library 'libcuda.so.1'; dlerror: libcuda.so.1: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /opt/hostedtoolcache/Python/3.7.13/x64/lib
2022-09-17 15:51:03.764638: W tensorflow/stream_executor/cuda/cuda_driver.cc:326] failed call to cuInit: UNKNOWN ERROR (303)
2022-09-17 15:51:03.764656: I tensorflow/stream_executor/cuda/cuda_diagnostics.cc:156] kernel driver does not appear to be running on this host (fv-az178-774): /proc/driver/nvidia/version does not exist
2022-09-17 15:51:03.764900: I tensorflow/core/platform/cpu_feature_guard.cc:142] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 AVX512F FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.
2022-09-17 15:51:03.765018: I tensorflow/compiler/jit/xla_gpu_device.cc:99] Not creating XLA devices, tf_xla_enable_xla_devices not set

Step 3: Build and optimise a model¶

[5]:

# Create a likelihood object
bernoulli_likelihood = Gaussian()

s2vgp = SparseVariationalGaussianProcess(
    kernel=kernel,
    inducing_points=tf.convert_to_tensor(inducing_points, dtype=default_float()),
    likelihood=bernoulli_likelihood
)

# equivalent to loss = -vgpc.elbo()
input_data = (time_points, observations)
loss = s2vgp.loss(input_data)

# Before optimisation, calculate the loss of the observations given the current kernel parameters
print("Loss before optimisation: ", loss.numpy())

Loss before optimisation:  447.3254994314457

[6]:

# Start at a small learning rate
adam_learning_rate = 0.005
natgrad_learning_rate = .9
max_iter = ci_niter(500)

adam_opt = tf.optimizers.Adam(learning_rate=adam_learning_rate)
natgrad_opt = SSMNaturalGradient(gamma=natgrad_learning_rate, momentum=False)
set_trainable(s2vgp.dist_q, False)
adam_var_list = s2vgp.trainable_variables
set_trainable(s2vgp.dist_q, True)


@tf.function
def loss(input_data):
    return -s2vgp.elbo(input_data)


@tf.function
def opt_step(input_data):
    natgrad_opt.minimize(lambda : loss(input_data), s2vgp.dist_q)
    adam_opt.minimize(lambda : loss(input_data), adam_var_list)

def plot_model(s2vgp):
    pred = s2vgp.posterior
    latent_mean, latent_var = pred.predict_f(tf.constant(time_points))
    predicted_mean, predicted_cov = latent_mean.numpy(), latent_var.numpy()
    # Plot the means and covariances for these future time points
    fig, ax = plt.subplots(1, 1)
    ax.plot(time_points, observations, 'C0x', ms=8, mew=2)

    ax.plot(time_points, predicted_mean, 'C0', lw=2)
    ax.fill_between(time_points,
                     predicted_mean[:, 0] - 2 * np.sqrt(predicted_cov[:, 0]),
                     predicted_mean[:, 0] + 2 * np.sqrt(predicted_cov[:, 0]),
                     color='C0', alpha=0.2)

    cp = s2vgp.kernel.change_points.numpy()
    ax.vlines(cp, ymin=-1, ymax=1, colors='blue', label='change points')
    z_ = s2vgp.inducing_inputs.numpy()
    ax.vlines(z_, ymin=-1, ymax=-.9, colors='red', label='inducing points')

    ax.set_xlim((0., 100.))
    ax.set_ylim((-1.1, 1.1))
    return fig
    # plt.show()


for i in range(max_iter):
    opt_step(input_data)
    if i % 10 == 0:
        print("Iteration:", i, ", Loss:", s2vgp.loss(input_data).numpy())
        fig = plot_model(s2vgp)
        plt.show()

print("Loss after optimisation: ", s2vgp.loss(input_data).numpy())

# Save our trained hyperparamters (these will be used in Step 8)
saved_hyperparams = kernel.trainable_variables

2022-09-17 15:51:04.174618: W tensorflow/python/util/util.cc:348] Sets are not currently considered sequences, but this may change in the future, so consider avoiding using them.
2022-09-17 15:51:41.857911: I tensorflow/compiler/mlir/mlir_graph_optimization_pass.cc:116] None of the MLIR optimization passes are enabled (registered 2)
2022-09-17 15:51:41.858373: I tensorflow/core/platform/profile_utils/cpu_utils.cc:112] CPU Frequency: 2793435000 Hz

Iteration: 0 , Loss: 316.48121407177354

../_images/notebooks_piecewise_kernels_10_2.png

Iteration: 10 , Loss: 311.41395165313884

../_images/notebooks_piecewise_kernels_10_4.png

Iteration: 20 , Loss: 306.5121307926455

../_images/notebooks_piecewise_kernels_10_6.png

Iteration: 30 , Loss: 301.60126848384016

../_images/notebooks_piecewise_kernels_10_8.png

Iteration: 40 , Loss: 296.66110968144653

../_images/notebooks_piecewise_kernels_10_10.png

Iteration: 50 , Loss: 291.67944022933386

../_images/notebooks_piecewise_kernels_10_12.png

Iteration: 60 , Loss: 286.65316373860423

../_images/notebooks_piecewise_kernels_10_14.png

Iteration: 70 , Loss: 281.5812215891556

../_images/notebooks_piecewise_kernels_10_16.png

Iteration: 80 , Loss: 276.46211834870815

../_images/notebooks_piecewise_kernels_10_18.png

Iteration: 90 , Loss: 271.29412214535876

../_images/notebooks_piecewise_kernels_10_20.png

Iteration: 100 , Loss: 266.0758733983721

../_images/notebooks_piecewise_kernels_10_22.png

Iteration: 110 , Loss: 260.80700595251994

../_images/notebooks_piecewise_kernels_10_24.png

Iteration: 120 , Loss: 255.48827315589872

../_images/notebooks_piecewise_kernels_10_26.png

Iteration: 130 , Loss: 250.12139105157087

../_images/notebooks_piecewise_kernels_10_28.png

Iteration: 140 , Loss: 244.7087704962197

../_images/notebooks_piecewise_kernels_10_30.png

Iteration: 150 , Loss: 239.25333065700448

../_images/notebooks_piecewise_kernels_10_32.png

Iteration: 160 , Loss: 233.75838202678153

../_images/notebooks_piecewise_kernels_10_34.png

Iteration: 170 , Loss: 228.22756254888898

../_images/notebooks_piecewise_kernels_10_36.png

Iteration: 180 , Loss: 222.66480143562507

../_images/notebooks_piecewise_kernels_10_38.png

Iteration: 190 , Loss: 217.07429836183192

../_images/notebooks_piecewise_kernels_10_40.png

Iteration: 200 , Loss: 211.46050808930363

../_images/notebooks_piecewise_kernels_10_42.png

Iteration: 210 , Loss: 205.82812713588945

../_images/notebooks_piecewise_kernels_10_44.png

Iteration: 220 , Loss: 200.1964052960365

../_images/notebooks_piecewise_kernels_10_46.png

Iteration: 230 , Loss: 194.5274992441604

../_images/notebooks_piecewise_kernels_10_48.png

Iteration: 240 , Loss: 188.8697076756959

../_images/notebooks_piecewise_kernels_10_50.png

Iteration: 250 , Loss: 183.21418993696676

../_images/notebooks_piecewise_kernels_10_52.png

Iteration: 260 , Loss: 177.56660146415135

../_images/notebooks_piecewise_kernels_10_54.png

Iteration: 270 , Loss: 171.9342211297267

../_images/notebooks_piecewise_kernels_10_56.png

Iteration: 280 , Loss: 166.3179655675469

../_images/notebooks_piecewise_kernels_10_58.png

Iteration: 290 , Loss: 160.72855958812212

../_images/notebooks_piecewise_kernels_10_60.png

Iteration: 300 , Loss: 155.17140266422572

../_images/notebooks_piecewise_kernels_10_62.png

Iteration: 310 , Loss: 149.65367797690644

../_images/notebooks_piecewise_kernels_10_64.png

Iteration: 320 , Loss: 144.16811541140817

../_images/notebooks_piecewise_kernels_10_66.png

Iteration: 330 , Loss: 138.73535029367162

../_images/notebooks_piecewise_kernels_10_68.png

Iteration: 340 , Loss: 133.35654309660904

../_images/notebooks_piecewise_kernels_10_70.png

Iteration: 350 , Loss: 128.03368053668186

../_images/notebooks_piecewise_kernels_10_72.png

Iteration: 360 , Loss: 122.76734239484284

../_images/notebooks_piecewise_kernels_10_74.png

Iteration: 370 , Loss: 117.56711170428748

../_images/notebooks_piecewise_kernels_10_76.png

Iteration: 380 , Loss: 112.44075175911198

../_images/notebooks_piecewise_kernels_10_78.png

Iteration: 390 , Loss: 107.36388144914186

../_images/notebooks_piecewise_kernels_10_80.png

Iteration: 400 , Loss: 102.35963533996237

../_images/notebooks_piecewise_kernels_10_82.png

Iteration: 410 , Loss: 97.42035593446765

../_images/notebooks_piecewise_kernels_10_84.png

Iteration: 420 , Loss: 92.54382335086257

../_images/notebooks_piecewise_kernels_10_86.png

Iteration: 430 , Loss: 87.70388786433914

../_images/notebooks_piecewise_kernels_10_88.png

Iteration: 440 , Loss: 82.91239664840502

../_images/notebooks_piecewise_kernels_10_90.png

Iteration: 450 , Loss: 78.16035532466704

../_images/notebooks_piecewise_kernels_10_92.png

Iteration: 460 , Loss: 73.39565087175617

../_images/notebooks_piecewise_kernels_10_94.png

Iteration: 470 , Loss: 68.63150700289388

../_images/notebooks_piecewise_kernels_10_96.png

Iteration: 480 , Loss: 63.83577330777081

../_images/notebooks_piecewise_kernels_10_98.png

Iteration: 490 , Loss: 59.01067313197379

../_images/notebooks_piecewise_kernels_10_100.png

Loss after optimisation:  54.536794585358486

[7]:

fig = plot_model(s2vgp)

../_images/notebooks_piecewise_kernels_11_0.png

We can see how our kernel parameters have changed from our initial values.

[8]:

gpflow.utilities.print_summary(s2vgp._kernel)

╒═════════════════════════════════════════╤═══════════╤═════════════╤═════════╤═════════════╤═════════╤═════════╤═══════════════════════════════════════╕
│ name                                    │ class     │ transform   │ prior   │ trainable   │ shape   │ dtype   │ value                                 │
╞═════════════════════════════════════════╪═══════════╪═════════════╪═════════╪═════════════╪═════════╪═════════╪═══════════════════════════════════════╡
│ PiecewiseKernel.kernels[0]._state_mean  │ Parameter │ Identity    │         │ True        │ (3,)    │ float64 │ [0. 0. 0.]                            │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[0]._lengthscale │ Parameter │ Softplus    │         │ True        │ ()      │ float64 │ 4.000000000008629                     │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[0]._variance    │ Parameter │ Softplus    │         │ True        │ ()      │ float64 │ 0.9999999995359001                    │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[1]._state_mean  │ Parameter │ Identity    │         │ True        │ (3,)    │ float64 │ [ 0.4964595   0.00954242 -0.0249826 ] │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[1]._lengthscale │ Parameter │ Softplus    │         │ True        │ ()      │ float64 │ 6.051791212765446                     │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[1]._variance    │ Parameter │ Softplus    │         │ True        │ ()      │ float64 │ 0.12410792782563587                   │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[2]._state_mean  │ Parameter │ Identity    │         │ True        │ (3,)    │ float64 │ [-0.16640493 -0.0189977   0.48716328] │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[2]._lengthscale │ Parameter │ Softplus    │         │ True        │ ()      │ float64 │ 3.806297312973849                     │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[2]._variance    │ Parameter │ Softplus    │         │ True        │ ()      │ float64 │ 0.2243925503862007                    │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[3]._state_mean  │ Parameter │ Identity    │         │ True        │ (3,)    │ float64 │ [-0.17769348 -0.22427976  0.68339092] │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[3]._lengthscale │ Parameter │ Softplus    │         │ True        │ ()      │ float64 │ 3.8330465978884303                    │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[3]._variance    │ Parameter │ Softplus    │         │ True        │ ()      │ float64 │ 0.2098822578898005                    │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[4]._state_mean  │ Parameter │ Identity    │         │ True        │ (3,)    │ float64 │ [ 0.00424717 -0.03013539 -0.23389754] │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[4]._lengthscale │ Parameter │ Softplus    │         │ True        │ ()      │ float64 │ 2.9559960686860784                    │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[4]._variance    │ Parameter │ Softplus    │         │ True        │ ()      │ float64 │ 0.2510613742421177                    │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[5]._state_mean  │ Parameter │ Identity    │         │ True        │ (3,)    │ float64 │ [0. 0. 0.]                            │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[5]._lengthscale │ Parameter │ Softplus    │         │ True        │ ()      │ float64 │ 4.0                                   │
├─────────────────────────────────────────┼───────────┼─────────────┼─────────┼─────────────┼─────────┼─────────┼───────────────────────────────────────┤
│ PiecewiseKernel.kernels[5]._variance    │ Parameter │ Softplus    │         │ True        │ ()      │ float64 │ 1.0                                   │
╘═════════════════════════════════════════╧═══════════╧═════════════╧═════════╧═════════════╧═════════╧═════════╧═══════════════════════════════════════╛

[ ]: