# A Simple Demonstration of Coregionalisation

This notebook will show how to construct a multi-output GP model using GPflow. We will consider a regression problem for functions $f: \mathbb{R}^D \rightarrow \mathbb{R}^P$. We assume that the dataset is of the form $(X_1, f_1), \dots, (X_P, f_P)$, i.e. we do not necessarily observe all the outputs for a given input location. In case we have fully observed outputs for each input, please refer to the [multioutput notebook](./multioutput.ipynb) for a more efficient implementation. For this problem, we model $f$ as a *coregionalised* Gaussian process, which assumes a kernel of the form:

$$\textrm{cov}(f_i(x), f_j(y)) = k_1(x, y) \times B[i, j].$$

The covariance of the $i$th function at $x$ and the $j$th function at $y$ is a kernel applied at $x$ and $y$, times the $i, j$th entry of a positive definite matrix $B$. This is known as the **intrinsic model of coregionalization (ICM)** _(Bonilla and Williams, 2008)_.

To make sure that B is positive-definite, we parameterise it as

$$B = W W^\top + \textrm{diag}(\kappa).$$

To build such model in GPflow, we need to perform the two following steps:

 * Create the kernel function defined above, using the `Coregion` kernel class.
 * Augment the training data X with an extra column containing an integer index to indicate which output an observation is associated with. This is essential to make the data work with the `Coregion` kernel.


```python
import gpflow
import numpy as np
import matplotlib
%matplotlib inline

matplotlib.rcParams['figure.figsize'] = (12, 6)
plt = matplotlib.pyplot
from gpflow.test_util import notebook_niter
np.random.seed(123)
```

## Data preparation
We start of by generating some training data to fit the model with. For this example, we choose the following two correlated functions for our outputs:
$$
\begin{align}
y_1 &= \sin(6x) + \epsilon_1, \qquad \epsilon_1 \sim \mathcal{N}(0, 0.009) \\
y_2 &= \sin(6x + 0.7) + \epsilon_2, \qquad \epsilon_2 \sim \mathcal{N}(0, 0.01) \\
\end{align}
$$


```python
# make a dataset with two outputs, correlated, heavy-tail noise. One has more noise than the other.
X1 = np.random.rand(100, 1) # Observed locations for first output
X2 = np.random.rand(50, 1) * 0.5 # Observed locations for second output

Y1 = np.sin(6*X1) + np.random.randn(*X1.shape) * 0.03
Y2 = np.sin(6*X2+ 0.7) + np.random.randn(*X2.shape) * 0.1

plt.figure(figsize=(8, 4))
plt.plot(X1, Y1, 'x', mew=2)
plt.plot(X2, Y2, 'x', mew=2);
```

## Data formatting for the Coregionalised model
To our training dataset, we add an extra column containing an index that specifies which output is observed.


```python
# Augment the input with ones or zeros to indicate the required output dimension
X_augmented = np.vstack((np.hstack((X1, np.zeros_like(X1))), np.hstack((X2, np.ones_like(X2)))))

# Augment the Y data
Y_augmented = np.vstack((np.hstack((Y1, np.zeros_like(Y1))), np.hstack((Y2, np.ones_like(Y2)))))
```

## Building the Coregionalisation kernel:
We build a Coregionalisation kernel with the Matern 3/2 kernel as the base kernel. This acts on the leading ([0]) data dimension of the augmented X values. The 'Coregion' kernel indexes the outputs, and acts on the last ([1]) data dimension (indices) of the augmented X values. To specify these dimensions, we use the built-in `active_dims` argument in the kernel constructor. To construct the full multi-output kernel, we take the product of the two kernels (For a more in-depth tutorial on kernel combination see this [notebook](./kernels.ipynb)).


```python
output_dim = 2 # Number of outputs
rank = 1 # Rank of W
k1 = gpflow.kernels.Matern32(1, active_dims=[0]) # Base kernel
coreg = gpflow.kernels.Coregion(1, output_dim=output_dim, rank=rank, active_dims=[1]) # Coregion kernel
coreg.W = np.random.rand(output_dim, rank) # Initialise the W matrix
kern = k1 * coreg
```

    WARNING:tensorflow:From /home/skitimoon/miniconda3/envs/gpflow/lib/python3.6/site-packages/tensorflow/python/framework/op_def_library.py:263: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.
    Instructions for updating:
    Colocations handled automatically by placer.


    WARNING:tensorflow:From /home/skitimoon/miniconda3/envs/gpflow/lib/python3.6/site-packages/tensorflow/python/framework/op_def_library.py:263: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.
    Instructions for updating:
    Colocations handled automatically by placer.


Note that, the default initialisation on W is zeroes; however, this is a saddle point in the objective so the value of W will not be optimised to fit the data. Hence, re-initializing the matrix to random entries should give a more accurate result.

## Constructing the model
The final ingredient to building the model is to specify the likelihood for each output dimension. To do this, we make use of the `SwitchedLikelihood` object in GPflow.


```python
kern
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>class</th>
      <th>prior</th>
      <th>transform</th>
      <th>trainable</th>
      <th>shape</th>
      <th>fixed_shape</th>
      <th>value</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>VGP/kern/kernels/0/lengthscales</th>
      <td>Parameter</td>
      <td>None</td>
      <td>+ve</td>
      <td>True</td>
      <td>()</td>
      <td>True</td>
      <td>0.7373221525197209</td>
    </tr>
    <tr>
      <th>VGP/kern/kernels/0/variance</th>
      <td>Parameter</td>
      <td>None</td>
      <td>+ve</td>
      <td>True</td>
      <td>()</td>
      <td>True</td>
      <td>0.3615150473405289</td>
    </tr>
    <tr>
      <th>VGP/kern/kernels/1/W</th>
      <td>Parameter</td>
      <td>None</td>
      <td>(none)</td>
      <td>True</td>
      <td>(2, 1)</td>
      <td>True</td>
      <td>[[-1.4482598862206812], [-2.323905683042864]]</td>
    </tr>
    <tr>
      <th>VGP/kern/kernels/1/kappa</th>
      <td>Parameter</td>
      <td>None</td>
      <td>+ve</td>
      <td>True</td>
      <td>(2,)</td>
      <td>True</td>
      <td>[1.3934457220562013, 0.14046778925072353]</td>
    </tr>
  </tbody>
</table>
</div>




```python
# This likelihood switches between Gaussian noise with different variances for each f_i:
lik = gpflow.likelihoods.SwitchedLikelihood([gpflow.likelihoods.Gaussian(), gpflow.likelihoods.StudentT()])

# now build the GP model as normal
m = gpflow.models.VGP(X_augmented, Y_augmented, kern=kern, likelihood=lik, num_latent=1)
# Here we specify num_latent=1 to avoid getting two outputs when predicting as Y_augmented is 2-dimensional

# fit the covariance function parameters
gpflow.train.ScipyOptimizer().minimize(m, maxiter=1000)
# gpflow.train.AdamOptimizer(0.01).minimize(m, maxiter=notebook_niter(1000))
```

    WARNING:tensorflow:From /home/skitimoon/miniconda3/envs/gpflow/lib/python3.6/site-packages/tensorflow/python/ops/math_ops.py:3066: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.
    Instructions for updating:
    Use tf.cast instead.


    WARNING:tensorflow:From /home/skitimoon/miniconda3/envs/gpflow/lib/python3.6/site-packages/tensorflow/python/ops/math_ops.py:3066: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.
    Instructions for updating:
    Use tf.cast instead.
    /home/skitimoon/miniconda3/envs/gpflow/lib/python3.6/site-packages/tensorflow/python/ops/gradients_impl.py:110: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.
      "Converting sparse IndexedSlices to a dense Tensor of unknown shape. "


    INFO:tensorflow:Optimization terminated with:
      Message: b'STOP: TOTAL NO. of ITERATIONS REACHED LIMIT'
      Objective function value: -224.699314
      Number of iterations: 1000
      Number of functions evaluations: 1099


    INFO:tensorflow:Optimization terminated with:
      Message: b'STOP: TOTAL NO. of ITERATIONS REACHED LIMIT'
      Objective function value: -224.699314
      Number of iterations: 1000
      Number of functions evaluations: 1099


That's it: the model has trained. Let's plot the model fit to see what's happened.


```python
def plot_gp(x, mu, var, color, label):
    plt.plot(x, mu, color=color, lw=2, label=label)
    plt.fill_between(x[:, 0], 
                     (mu - 2*np.sqrt(var))[:, 0], 
                     (mu + 2*np.sqrt(var))[:, 0], 
                     color=color, alpha=0.4)

def plot(m):
    plt.figure(figsize=(8, 4))
    xtest = np.linspace(0, 1, 100)[:,None]
    line, = plt.plot(X1, Y1, 'x', mew=2)
    mu, var = m.predict_f(np.hstack((xtest, np.zeros_like(xtest))))
    plot_gp(xtest, mu, var, line.get_color(), 'Y1')

    line, = plt.plot(X2, Y2, 'x', mew=2)
    mu, var = m.predict_f(np.hstack((xtest, np.ones_like(xtest))))
    plot_gp(xtest, mu, var, line.get_color(), 'Y2')
    
    plt.legend()

plot(m)
```

From the plots we see:

 - The first function (blue) has low posterior variance everywhere because there are so many observations, and the noise variance is small. 
 - The second function (orange) has higher posterior variance near the data, because the data are more noisy, and very high posterior variance where there are no observations (x > 0.5). 
 - The model has done a reasonable job of estimating the noise variance and lengthscales.
 - The model recognises the correlation between the two functions and is able to suggest (with uncertainty) that as x > 0.5 the orange curve should follow the blue curve (which we know to be the truth from the data generating procedure above).
 
The covariance matrix between outputs is as follows:


```python
kern
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>class</th>
      <th>prior</th>
      <th>transform</th>
      <th>trainable</th>
      <th>shape</th>
      <th>fixed_shape</th>
      <th>value</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>VGP/kern/kernels/0/lengthscales</th>
      <td>Parameter</td>
      <td>None</td>
      <td>+ve</td>
      <td>True</td>
      <td>()</td>
      <td>True</td>
      <td>0.7373221525197209</td>
    </tr>
    <tr>
      <th>VGP/kern/kernels/0/variance</th>
      <td>Parameter</td>
      <td>None</td>
      <td>+ve</td>
      <td>True</td>
      <td>()</td>
      <td>True</td>
      <td>0.3615150473405289</td>
    </tr>
    <tr>
      <th>VGP/kern/kernels/1/W</th>
      <td>Parameter</td>
      <td>None</td>
      <td>(none)</td>
      <td>True</td>
      <td>(2, 1)</td>
      <td>True</td>
      <td>[[-1.4482598862206812], [-2.323905683042864]]</td>
    </tr>
    <tr>
      <th>VGP/kern/kernels/1/kappa</th>
      <td>Parameter</td>
      <td>None</td>
      <td>+ve</td>
      <td>True</td>
      <td>(2,)</td>
      <td>True</td>
      <td>[1.3934457220562013, 0.14046778925072353]</td>
    </tr>
  </tbody>
</table>
</div>




```python
B = coreg.W.value @ coreg.W.value.T + np.diag(coreg.kappa.value)
print('B =', B)
plt.imshow(B);
```

## Further Reading:
- [Kernels](../advanced/kernels.ipynb), explaining advanced usage of GPflow kernels
- [Multioutput](./multioutput.ipynb), which details other GPflow features for multi-output prediction.

## References:

Bonilla, Edwin V., Kian M. Chai, and Christopher Williams. "Multi-task Gaussian process prediction." Advances in neural information processing systems. 2008.

00df87203c930dd2e345da1b94cd78b61226b232
ipynb
Jupyter Notebook
coregionalisation.ipynb
skitimoon/gpflow-examples
0e9eee0b6d68b9f811882919b2bf01ecae77a285
MIT
coregionalisation.ipynb
skitimoon/gpflow-examples
0e9eee0b6d68b9f811882919b2bf01ecae77a285
MIT
coregionalisation.ipynb
skitimoon/gpflow-examples
0e9eee0b6d68b9f811882919b2bf01ecae77a285
MIT
Qwen/Qwen-72B
1. YES
2. YES