In [1]:
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import scipy.io as spio

Bayesian linear regression¶

In this exercise session, we consider the supervised regression problem of finding a function $f(x)$ that describes the relationship between a scalar input $x$ and a scalar output $y$:

$$ y = f(x) + \epsilon, \qquad \epsilon \sim \mathcal{N}(0, \beta^{-1}). $$

We model this with a Bayesian linear regression model

$$ f(x) = \boldsymbol{\phi}(x)^{\mathsf{T}} \mathbf{w}, \qquad \mathbf{w} \sim \mathcal{N}(\mathbf{m}_0, \mathbf{S}_0), $$

where $\boldsymbol{\phi}(x)$ is a vector of the input features. Note that we used the notation $\mathbf{x}$ for the input features in the lecture. We changed this notation to $\boldsymbol{\phi}(x)$ here in order to not mix it up with the scalar input $x$.

The Bayesian linear regression model is then given by

$$ \begin{aligned} p(\mathbf{y} \,|\, \mathbf{w}) &= \mathcal{N}(\mathbf{y}; \boldsymbol{\Phi}\mathbf{w}, \beta^{-1}\mathbf{I}_N) \qquad && \text{(likelihood)}, \\ p(\mathbf{w}) &= \mathcal{N}(\mathbf{w}; \mathbf{m}_0, \mathbf{S}_0) \qquad && \text{(prior)}, \end{aligned} $$

where

$$ \boldsymbol{\Phi} = \begin{bmatrix} \boldsymbol{\phi}(x_1)^{\mathsf{T}} \\ \vdots \\ \boldsymbol{\phi}(x_N)^{\mathsf{T}} \end{bmatrix} \qquad \text{and} \qquad \mathbf{y} = \begin{bmatrix} y_1 \\ \vdots \\ y_N \end{bmatrix}. $$

Given a set of training data of inputs and outputs $\mathcal{D} = \{(x_i, y_i)\}_{i=1}^N$, we are interested in finding the posterior of the weights $p(\mathbf{w} \,|\, \mathbf{y})$ and also the predictive distribution $p(f(x_{\star}) \,|\, \mathbf{y})$ of an unseen input $x_{\star}$. For further information about the Bayesian linear regression model, see Lecture 3 and/or Christopher Bishop's book "Pattern recognition and machine learning".

Exercise 3.1: Understanding the code¶

Download the files `lindata.mat` and `nlindata.mat` and save them to the folder of this notebook. These datasets are borrowed from Philipp Hennig's course "Probabilistic machine learning", given at the University of Tübingen.

The following code cell loads inputs, outputs, and the precision parameter from lindata.mat and plots the feature vector

$$ \boldsymbol{\phi}(x)^{\mathsf{T}} = [1, x]. $$
In [2]:
# Load data from disk
# File should be in the same folder as the Jupyter notebook,
# otherwise you have to adjust the path
data = spio.loadmat("lindata.mat")
x = data["X"].flatten() # inputs
y = data["Y"].flatten() # outputs
beta = float(data["sigma"].item())**(-2) # measurement noise precision
N = x.size

# Define the feature vector
def Phi(a):  # Phi(a) = [1, a]
    return np.power(np.reshape(a, (-1, 1)), range(2))

# Plot the features
plt.plot(x, Phi(x), '-o')
plt.title('features')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

We then compute the posterior distribution of the weights $\mathbf{w}$ of a Bayesian linear regression model using these features.

In [3]:
# Define the prior on the weights
# p(w) = N(w; m0, S0)
D = Phi(0).size  # number of features
m0 = np.zeros(D)
S0 = 10*np.eye(D) / D

# Compute the posterior distribution of the Bayesian linear regression model
# p(w | y) = N(w; mN, SN)
SN = np.linalg.inv(np.linalg.inv(S0) + beta * Phi(x).T @ Phi(x))
mN = SN @ (np.linalg.inv(S0) @ m0 + beta * Phi(x).T @ y)

We visualize the posterior distribution by plotting the functions $f$ corresponding to different samples of $\mathbf{w}$.

In [4]:
# Generate grid of new inputs x* for plotting
n = 100  # number of grid-points
xs = np.linspace(-8, 8, n)

# Visualize the posterior p(w | y) = N(w; mN, SN)
# For samples of w, f(x) = phi(x)^T w is evaluated at inputs xs
# Draw samples of w from the posterior
samples = 5
seed = 100
ws = stats.multivariate_normal(mean=mN, cov=SN, allow_singular=True).rvs(samples, random_state=seed)

# Compute corresponding values f(x*)
fs = Phi(xs) @ ws.T

# Plot the samples
plt.plot(xs, fs, 'gray') # samples
plt.scatter(x, y, zorder=3)
plt.title('posterior - samples')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Next, we plot samples from and credibility regions of the predictive distribution.

In [5]:
# Compute the predictive distribution of the outputs y*
# p(y* | y) = N(y*; m*, S*)
mstar = Phi(xs) @ mN
Sstar = Phi(xs) @ SN @ Phi(xs).T + beta**(-1) * np.eye(n)

# Extract standard deviation of predictive distribution
stdpred= np.sqrt(np.diag(Sstar))

# Plot credibility regions
plt.plot(xs, mstar, 'black') # predictive mean
plt.fill_between(xs, mstar + 3*stdpred, mstar - 3*stdpred, color='lightgray')
plt.fill_between(xs, mstar + 2*stdpred, mstar - 2*stdpred, color='darkgray')
plt.fill_between(xs, mstar + 1*stdpred, mstar - 1*stdpred, color='gray')
plt.scatter(x, y, zorder=3)
plt.title('predictive distribution')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

(a)¶

Go through the code and make sure that you can map the code to the model and regression method explained in Lecture 3. Also run the code. Make sure you understand the figures.

(b)¶

Reduce the training data to only the first 5 data points in the training data. What impact does this have on the predictive distribution?

In [6]:
# Define the inputs and outputs
N = 5
x = data["X"][:N].flatten() # inputs
y = data["Y"][:N].flatten() # outputs

# We use the same feature vector and prior as above

# Compute the posterior distribution of the Bayesian linear regression model
# p(w | y) = N(w; mN, SN)
SN = np.linalg.inv(np.linalg.inv(S0) + beta * Phi(x).T @ Phi(x))
mN = SN @ (np.linalg.inv(S0) @ m0 + beta * Phi(x).T @ y)

# Compute the predictive distribution of the outputs y*
# p(y* | y) = N(y*; m*, S*)
mstar = Phi(xs) @ mN
Sstar = Phi(xs) @ SN @ Phi(xs).T + beta**(-1) * np.eye(n)

# Extract standard deviation of predictive distribution
stdpred= np.sqrt(np.diag(Sstar))

# Plot credibility regions
plt.plot(xs, mstar, 'black') # predictive mean
plt.fill_between(xs, mstar + 3*stdpred, mstar - 3*stdpred, color='lightgray')
plt.fill_between(xs, mstar + 2*stdpred, mstar - 2*stdpred, color='darkgray')
plt.fill_between(xs, mstar + 1*stdpred, mstar - 1*stdpred, color='gray')
plt.scatter(x, y, zorder=3)
plt.title('predictive distribution')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Now the model becomes very uncertain in regions where we don't have any measurements.

Exercise 3.2: Feature vectors¶

(a)¶

Load nlindata.mat instead of lindata.mat and run the code for this data. Use all data, not only the first five data points as in Exercise 3.1 (b). Do you think the model performs well on this data?

In [7]:
# Load data from disk
# File should be in the same folder as the Jupyter notebook,
# otherwise you have to adjust the path
data = spio.loadmat("nlindata.mat")
x = data["X"].flatten() # inputs
y = data["Y"].flatten() # outputs
beta = float(data["sigma"].item())**(-2) # measurement noise precision
N = x.size

# Plot the features
# We use the same feature vector as above
plt.plot(x, Phi(x), '-o')
plt.title('features')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
In [8]:
# Compute the posterior distribution of the Bayesian linear regression model
# p(w | y) = N(w; mN, SN)
# We use the same prior as above
SN = np.linalg.inv(np.linalg.inv(S0) + beta * Phi(x).T @ Phi(x))
mN = SN @ (np.linalg.inv(S0) @ m0 + beta * Phi(x).T @ y)

# Visualize the posterior p(w | y) = N(w; mN, SN)
# For samples of w, f(x) = phi(x)^T w is evaluated at inputs xs
# Draw samples of w from the posterior
ws = stats.multivariate_normal(mean=mN, cov=SN, allow_singular=True).rvs(samples, random_state=seed)

# Compute corresponding values f(x*)
fs = Phi(xs) @ ws.T

# Plot the samples
plt.plot(xs, fs, 'gray') # samples
plt.scatter(x, y, zorder=3)
plt.title('posterior - samples')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
In [9]:
# Compute the predictive distribution of the outputs y*
# p(y* | y) = N(y*; m*, S*)
mstar = Phi(xs) @ mN
Sstar = Phi(xs) @ SN @ Phi(xs).T + beta**(-1) * np.eye(n)

# Extract standard deviation of predictive distribution
stdpred= np.sqrt(np.diag(Sstar))

# Plot credibility regions
plt.plot(xs, mstar, 'black') # predictive mean
plt.fill_between(xs, mstar + 3*stdpred, mstar - 3*stdpred, color='lightgray')
plt.fill_between(xs, mstar + 2*stdpred, mstar - 2*stdpred, color='darkgray')
plt.fill_between(xs, mstar + 1*stdpred, mstar - 1*stdpred, color='gray')
plt.scatter(x, y, zorder=3)
plt.title('predictive distribution')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

For this data, a feature vector with only a constant and linear term does not perform that well.

(b)¶

In order to improve the performance, consider instead a feature vector with an additional quadratic term

$$ \boldsymbol{\phi}(x)^{\mathsf{T}} = [1, x, x^2]. $$

Change the code accordingly and run it.

Hint: Only a very minor modification in the code is required to accommodate this change.

Only the definition of the features Phi has to be changed. An additional column for the quadratic term is required.

In [10]:
# Define the feature vector with an additional quadratic term
def Phi(a):  # Phi(a) = [1, a, a**2]
    return np.power(np.reshape(a, (-1, 1)), range(3))

# Plot the features
plt.plot(x, Phi(x), '-o')
plt.title('features')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
In [11]:
# We have to redefine the prior on the weights since the
# dimension of the feature vector changed
# p(w) = N(w; m0, S0)
D = Phi(0).size  # number of features
m0 = np.zeros(D)
S0 = 10*np.eye(D) / D

# Compute the posterior distribution of the Bayesian linear regression model
# p(w | y) = N(w; mN, SN)
SN = np.linalg.inv(np.linalg.inv(S0) + beta * Phi(x).T @ Phi(x))
mN = SN @ (np.linalg.inv(S0) @ m0 + beta * Phi(x).T @ y)

# Visualize the posterior p(w | y) = N(w; mN, SN)
# For samples of w, f(x) = phi(x)^T w is evaluated at inputs xs
# Draw samples of w from the posterior
ws = stats.multivariate_normal(mean=mN, cov=SN, allow_singular=True).rvs(samples, random_state=seed)

# Compute corresponding values f(x*)
fs = Phi(xs) @ ws.T

# Plot the samples
plt.plot(xs, fs, 'gray') # samples
plt.scatter(x, y, zorder=3)
plt.title('posterior - samples')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
In [12]:
# Compute the predictive distribution of the outputs y*
# p(y* | y) = N(y*; m*, S*)
mstar = Phi(xs) @ mN
Sstar = Phi(xs) @ SN @ Phi(xs).T + beta**(-1) * np.eye(n)

# Extract standard deviation of predictive distribution
stdpred= np.sqrt(np.diag(Sstar))

# Plot credibility regions
plt.plot(xs, mstar, 'black') # predictive mean
plt.fill_between(xs, mstar + 3*stdpred, mstar - 3*stdpred, color='lightgray')
plt.fill_between(xs, mstar + 2*stdpred, mstar - 2*stdpred, color='darkgray')
plt.fill_between(xs, mstar + 1*stdpred, mstar - 1*stdpred, color='gray')
plt.scatter(x, y, zorder=3)
plt.title('predictive distribution')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

(c)¶

Use step functions

$$ h(x) = \begin{cases} 1, \qquad x \geq 0,\\ 0, \qquad x < 0, \end{cases} $$

as features and change the code accordingly. Place a total of 9 of these features, with the steps two units apart, between $x = -8$ and $x = 8$. The feature vector is then

$$ \boldsymbol{\phi}(x)^{\mathsf{T}} = [h(x - 8), h(x - 6), \ldots, h(x + 8)]. $$

Again only the definition of the features Phi has to be changed.

In [13]:
# Define the feature vector with step functions
def Phi(a):  # Phi(a) = [h(a - 8), h(a - 6), ..., h(a + 8)]
    return (np.reshape(a, (-1, 1)) > np.linspace(-8, 8, 9))

# Plot the features
plt.plot(x, Phi(x), '-o')
plt.title('features')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
In [14]:
# We have to redefine the prior on the weights since the
# dimension of the feature vector changed
# p(w) = N(w; m0, S0)
D = Phi(0).size  # number of features
m0 = np.zeros(D)
S0 = 10*np.eye(D) / D

# Compute the posterior distribution of the Bayesian linear regression model
# p(w | y) = N(w; mN, SN)
SN = np.linalg.inv(np.linalg.inv(S0) + beta * Phi(x).T @ Phi(x))
mN = SN @ (np.linalg.inv(S0) @ m0 + beta * Phi(x).T @ y)

# Visualize the posterior p(w | y) = N(w; mN, SN)
# For samples of w, f(x) = phi(x)^T w is evaluated at inputs xs
# Draw samples of w from the posterior
ws = stats.multivariate_normal(mean=mN, cov=SN, allow_singular=True).rvs(samples, random_state=seed)

# Compute corresponding values f(x*)
fs = Phi(xs) @ ws.T

# Plot the samples
plt.plot(xs, fs, 'gray') # samples
plt.scatter(x, y, zorder=3)
plt.title('posterior - samples')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
In [15]:
# Compute the predictive distribution of the outputs y*
# p(y* | y) = N(y*; m*, S*)
mstar = Phi(xs) @ mN
Sstar = Phi(xs) @ SN @ Phi(xs).T + beta**(-1) * np.eye(n)

# Extract standard deviation of predictive distribution
stdpred= np.sqrt(np.diag(Sstar))

# Plot credibility regions
plt.plot(xs, mstar, 'black') # predictive mean
plt.fill_between(xs, mstar + 3*stdpred, mstar - 3*stdpred, color='lightgray')
plt.fill_between(xs, mstar + 2*stdpred, mstar - 2*stdpred, color='darkgray')
plt.fill_between(xs, mstar + 1*stdpred, mstar - 1*stdpred, color='gray')
plt.scatter(x, y, zorder=3)
plt.title('predictive distribution')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

(d)¶

Can you come up with any other features that improve performance even further?

One can try, e.g., piecewise linear features:

In [16]:
# Define the feature vector with piecewise linear functions
def Phi(a):  # Phi(a) = [|a + 8| + 8, |a + 7| + 7, ..., |a - 8| - 8]
    return np.abs(np.reshape(a, (-1, 1)) - np.linspace(-8, 8, 17)) - np.linspace(-8, 8, 17)

# Plot the features
plt.plot(x, Phi(x), '-o')
plt.title('features')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
In [17]:
# We have to redefine the prior on the weights since the
# dimension of the feature vector changed
# p(w) = N(w; m0, S0)
D = Phi(0).size  # number of features
m0 = np.zeros(D)
S0 = 10*np.eye(D) / D

# Compute the posterior distribution of the Bayesian linear regression model
# p(w | y) = N(w; mN, SN)
SN = np.linalg.inv(np.linalg.inv(S0) + beta * Phi(x).T @ Phi(x))
mN = SN @ (np.linalg.inv(S0) @ m0 + beta * Phi(x).T @ y)

# Visualize the posterior p(w | y) = N(w; mN, SN)
# For samples of w, f(x) = phi(x)^T w is evaluated at inputs xs
# Draw samples of w from the posterior
ws = stats.multivariate_normal(mean=mN, cov=SN, allow_singular=True).rvs(samples, random_state=seed)

# Compute corresponding values f(x*)
fs = Phi(xs) @ ws.T

# Plot the samples
plt.plot(xs, fs, 'gray') # samples
plt.scatter(x, y, zorder=3)
plt.title('posterior - samples')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
In [18]:
# Compute the predictive distribution of the outputs y*
# p(y* | y) = N(y*; m*, S*)
mstar = Phi(xs) @ mN
Sstar = Phi(xs) @ SN @ Phi(xs).T + beta**(-1) * np.eye(n)

# Extract standard deviation of predictive distribution
stdpred= np.sqrt(np.diag(Sstar))

# Plot credibility regions
plt.plot(xs, mstar, 'black') # predictive mean
plt.fill_between(xs, mstar + 3*stdpred, mstar - 3*stdpred, color='lightgray')
plt.fill_between(xs, mstar + 2*stdpred, mstar - 2*stdpred, color='darkgray')
plt.fill_between(xs, mstar + 1*stdpred, mstar - 1*stdpred, color='gray')
plt.scatter(x, y, zorder=3)
plt.title('predictive distribution')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Exercise 3.3: Marginal likelihood¶

To get a quantitative measure of the performance of the proposed feature vectors, we want to compare them by computing the marginal likelihood $p(\mathbf{y})$ for each of the models. Refer to Exercise 2.9(a) for the expression of the marginal likelihood of the Bayesian linear regression model.

(a)¶

Extend the code to also compute the logarithm of the marginal likelihood. Which one of the four feature vectors in Exercise 3.2 gives the largest log marginal likelihood on the data nlindata.mat?

The following function computes the log marginal likelihood:

In [19]:
# Compute the log-likelihood of the marginal distribution
# p(y) = N(y; my, Sy)
# for the inputs X and outputs y in a Bayesian regression model
# p(y | w, X) = N(y; Xw, I/beta)
# p(w) = N(w; m0, S0)
def marginal_loglik(m0, S0, beta, X, y):
    N = X.shape[0]
    my = X @ m0
    Sy = X @ S0 @ X.T + beta**(-1) * np.eye(N)
    rv = stats.multivariate_normal(mean=my, cov=Sy, allow_singular=True)
    return rv.logpdf(y)
In [20]:
# Load data from disk
# File should be in the same folder as the Jupyter notebook,
# otherwise you have to adjust the path
data = spio.loadmat("nlindata.mat")
x = data["X"].flatten() # inputs
y = data["Y"].flatten() # outputs
beta = float(data["sigma"].item())**(-2) # measurement noise precision

# Feature vector and prior used in 3.2 (a)
def Phi(a):  # Phi(a) = [1, a]
    return np.power(np.reshape(a, (-1, 1)), range(2))
D = Phi(0).size  # number of features
m0 = np.zeros(D)
S0 = 10*np.eye(D) / D
print("Linear (3.2 (a)):\t\t{:8.1f}".format(marginal_loglik(m0, S0, beta, Phi(x), y)))

# Feature vector and prior used in 3.2 (b)
def Phi(a):  # Phi(a) = [1, a, a**2]
    return np.power(np.reshape(a, (-1, 1)), range(3))
D = Phi(0).size  # number of features
m0 = np.zeros(D)
S0 = 10*np.eye(D) / D
print("Quadratic (3.2 (b)):\t\t{:8.1f}".format(marginal_loglik(m0, S0, beta, Phi(x), y)))

# Feature vector and prior used in 3.2 (c)
def Phi(a):  # Phi(a) = [h(a - 8), h(a - 6), ..., h(a + 8)]
    return (np.reshape(a, (-1, 1)) > np.linspace(-8, 8, 9))
D = Phi(0).size  # number of features
m0 = np.zeros(D)
S0 = 10*np.eye(D) / D
print("Step function (3.2 (c)):\t{:8.1f}".format(marginal_loglik(m0, S0, beta, Phi(x), y)))

# Define the feature vector with piecewise linear functions
def Phi(a):  # Phi(a) = [|a + 8| + 8, |a + 7| + 7, ..., |a - 8| - 8]
    return np.abs(np.reshape(a, (-1, 1)) - np.linspace(-8, 8, 17)) - np.linspace(-8, 8, 17)
D = Phi(0).size  # number of features
m0 = np.zeros(D)
S0 = 10*np.eye(D) / D
print("Piecewise linear (3.2 (d)):\t{:8.1f}".format(marginal_loglik(m0, S0, beta, Phi(x), y)))

Linear (3.2 (a)):		  -861.3
Quadratic (3.2 (b)):		  -202.5
Step function (3.2 (c)):	  -382.6
Piecewise linear (3.2 (d)):	   -87.1

We see that linear features perform worst in terms of marginal likelihood and that the piecewise linear features show the largest marginal likelihood.

(b)¶

Perform the same comparison on the data lindata.mat. What are your conclusions?

For this data, linear features perform best in terms of marginal likelihood. Note that the model with the linear features is a special case of the models with the quadratic and with the piecewise linear features. Hence we see that the marginal likelihood also penalizes model complexity.

In [21]:
# Load data from disk
# File should be in the same folder as the Jupyter notebook,
# otherwise you have to adjust the path
data = spio.loadmat("lindata.mat")
x = data["X"].flatten() # inputs
y = data["Y"].flatten() # outputs
beta = float(data["sigma"].item())**(-2) # measurement noise precision

# Feature vector and prior used in 3.2 (a)
def Phi(a):  # Phi(a) = [1, a]
    return np.power(np.reshape(a, (-1, 1)), range(2))
D = Phi(0).size  # number of features
m0 = np.zeros(D)
S0 = 10*np.eye(D) / D
print("Linear (3.2 (a)):\t\t{:8.1f}".format(marginal_loglik(m0, S0, beta, Phi(x), y)))

# Feature vector and prior used in 3.2 (b)
def Phi(a):  # Phi(a) = [1, a, a**2]
    return np.power(np.reshape(a, (-1, 1)), range(3))
D = Phi(0).size  # number of features
m0 = np.zeros(D)
S0 = 10*np.eye(D) / D
print("Quadratic (3.2 (b)):\t\t{:8.1f}".format(marginal_loglik(m0, S0, beta, Phi(x), y)))

# Feature vector and prior used in 3.2 (c)
def Phi(a):  # Phi(a) = [h(a - 8), h(a - 6), ..., h(a + 8)]
    return (np.reshape(a, (-1, 1)) > np.linspace(-8, 8, 9))
D = Phi(0).size  # number of features
m0 = np.zeros(D)
S0 = 10*np.eye(D) / D
print("Step function (3.2 (c)):\t{:8.1f}".format(marginal_loglik(m0, S0, beta, Phi(x), y)))

# Define the feature vector with piecewise linear functions
def Phi(a):  # Phi(a) = [|a + 8| + 8, |a + 7| + 7, ..., |a - 8| - 8]
    return np.abs(np.reshape(a, (-1, 1)) - np.linspace(-8, 8, 17)) - np.linspace(-8, 8, 17)
D = Phi(0).size  # number of features
m0 = np.zeros(D)
S0 = 10*np.eye(D) / D
print("Piecewise linear (3.2 (d)):\t{:8.1f}".format(marginal_loglik(m0, S0, beta, Phi(x), y)))

Linear (3.2 (a)):		   -57.9
Quadratic (3.2 (b)):		   -61.1
Step function (3.2 (c)):	   -71.1
Piecewise linear (3.2 (d)):	   -62.5

(c)¶

Can you come up with any other feature vectors and/or values for prior/likelihood precisions that give an even larger log marginal likelihood?

Open-ended question...

Exercise 3.4: Challenge Problem: From weights to functions¶

In the previous exercises, we described a random function through Gaussian weights. Can we make the same predictions by working directly with the function values? We will first recover the quadratic model this way, then use the new representation to go beyond a finite feature vector.

Use nlindata.mat and the zero-mean weight prior below.

In [22]:
# Load the nonlinear data
data = spio.loadmat("nlindata.mat")
x = data["X"].flatten()
y = data["Y"].flatten()
beta = float(np.asarray(data["sigma"]).item())**(-2)

# Define the quadratic features and the weight prior
def Phi_quad(a):
    return np.asarray(a).reshape(-1, 1)**np.arange(3)

S0 = (10 / 3) * np.eye(3)
rng = np.random.default_rng(7)
x_grid = np.linspace(-8, 8, 120)

(a) The same prior viewed in two ways¶

Previously, we sampled weights and used them to draw a curve. Consider instead the values of that curve at $M$ fixed grid inputs $x_{g,1},\ldots,x_{g,M}$,

$$ \mathbf f_g = \begin{bmatrix} f(x_{g,1})\\ \vdots\\ f(x_{g,M}) \end{bmatrix} = \underbrace{\begin{bmatrix} \boldsymbol\phi(x_{g,1})^{\mathsf T}\\ \vdots\\ \boldsymbol\phi(x_{g,M})^{\mathsf T} \end{bmatrix}}_{\boldsymbol\Phi_g\,\in\,\mathbb R^{M\times D}} \mathbf w, \qquad \mathbf w\sim\mathcal N(\mathbf 0,\mathbf S_0). $$

One draw of $\mathbf f_g\in\mathbb R^M$ gives $M$ values from one random curve. They depend on the same weights and are generally dependent.

Find the distribution of $\mathbf f_g$, including its mean and covariance. Write its covariance entries as $k(x,x')=\operatorname{Cov}(f(x),f(x'))$ and complete k_quadratic. For arrays a and b, this function should return the matrix with entries $k(a_i, b_j)$.

Compare the two sampling procedures below. The right panel samples a vector with 120 entries without sampling any weights. Why do its entries still trace out a quadratic curve?

In [23]:
def k_quadratic(a, b):
    # Use Phi_quad and S0
    # Phi_quad handles the reshaping of each input array
    K = ...
    K = Phi_quad(a) @ S0 @ Phi_quad(b).T
    return K

# Sample four curves through the weights
weights = rng.multivariate_normal(np.zeros(3), S0, size=4)
f_from_weights = (Phi_quad(x_grid) @ weights.T).T

# Sample four vectors of function values directly
K_grid = k_quadratic(x_grid, x_grid)
# Tiny diagonal jitter protects against numerical round-off
f_from_covariance = rng.multivariate_normal(
    np.zeros(len(x_grid)), K_grid + 1e-10 * np.eye(len(x_grid)), size=4
)

fig, axes = plt.subplots(1, 2, figsize=(10, 3.5), sharey=True)
axes[0].plot(x_grid, f_from_weights.T)
axes[0].set_title("Sample weights, then evaluate")
axes[1].plot(x_grid, f_from_covariance.T)
axes[1].set_title("Sample function values directly")
for ax in axes:
    ax.set_xlabel("x")
axes[0].set_ylabel("f(x)")
plt.tight_layout()
plt.show()

Solution to (a)

A linear transformation of a Gaussian vector is Gaussian, so

$$ \mathbf f_g\sim\mathcal N(\mathbf 0,\mathbf K_g), \qquad \mathbf K_g=\boldsymbol\Phi_g\mathbf S_0\boldsymbol\Phi_g^{\mathsf T}, \qquad k(x,x')=\boldsymbol\phi(x)^{\mathsf T}\mathbf S_0\boldsymbol\phi(x'). $$

The covariance carries the restrictions previously expressed through the weights. Although $\mathbf f_g$ has 120 entries, it lies in the span of the three columns of $\boldsymbol\Phi_g$. Its entries must therefore be evaluations of a quadratic curve. Sampling each entry independently would discard these restrictions.

The two procedures generate the same distribution over grid values, though these independent draws need not look identical. The derivation establishes their equivalence and the plots illustrate it. The tiny jitter in the code adds a negligible numerical perturbation to the second distribution.

(b) Prediction without explicit weights¶

The calculation in (a) holds for any finite collection of inputs. A random function whose values at every such collection are jointly Gaussian is, by definition, a Gaussian process. Our quadratic model is therefore already a Gaussian process. We now use its covariance function to make predictions.

Let $X$ contain the $N$ training inputs and $X_\star$ contain $M$ test inputs. Collect the covariances into matrices,

$$ \begin{aligned} (\mathbf K_{XX})_{ij} &= k(x_i,x_j), &\mathbf K_{XX}&\in\mathbb R^{N\times N},\\ (\mathbf K_{\star X})_{ij} &= k(x_{\star,i},x_j), &\mathbf K_{\star X}&\in\mathbb R^{M\times N},\\ (\mathbf K_{\star\star})_{ij} &= k(x_{\star,i},x_{\star,j}), &\mathbf K_{\star\star}&\in\mathbb R^{M\times M}. \end{aligned} $$

The weight posterior has mean

$$ \mathbf m_N=\beta(\mathbf S_0^{-1}+\beta\boldsymbol\Phi^{\mathsf T}\boldsymbol\Phi)^{-1} \boldsymbol\Phi^{\mathsf T}\mathbf y. $$

Use the following form of the matrix inversion lemma to rewrite the predictive mean $\boldsymbol\Phi_\star\mathbf m_N$ using only these covariance matrices, $\beta$, and $\mathbf y$,

$$ (\mathbf A^{-1}+\mathbf B^{\mathsf T}\mathbf C^{-1}\mathbf B)^{-1} \mathbf B^{\mathsf T}\mathbf C^{-1} =\mathbf A\mathbf B^{\mathsf T} (\mathbf B\mathbf A\mathbf B^{\mathsf T}+\mathbf C)^{-1}. $$

For the code, you may use the corresponding predictive covariance,

$$ \operatorname{Cov}(\mathbf f_\star\mid\mathbf y) =\mathbf K_{\star\star} -\mathbf K_{\star X}(\mathbf K_{XX}+\beta^{-1}\mathbf I_N)^{-1}\mathbf K_{\star X}^{\mathsf T}. $$

Complete predict_from_kernel and check that it reproduces the quadratic model's mean and covariance. Use np.linalg.solve(A, b) to compute $\mathbf A^{-1}\mathbf b$.

In [24]:
def predict_from_kernel(k, x_train, y_train, beta, x_test):
    # K has shape (N, N), Ks has shape (M, N), Kss has shape (M, M)
    K = ...
    Ks = ...
    Kss = ...
    A = ...
    mean = ...
    cov = ...
    K = k(x_train, x_train)
    Ks = k(x_test, x_train)
    Kss = k(x_test, x_test)
    A = K + np.eye(len(x_train)) / beta
    mean = Ks @ np.linalg.solve(A, y_train)
    cov = Kss - Ks @ np.linalg.solve(A, Ks.T)
    return mean, (cov + cov.T) / 2

# Compute the familiar weight posterior and its predictions
Phi_X = Phi_quad(x)
Phi_s = Phi_quad(x_grid)
SN = np.linalg.inv(np.linalg.inv(S0) + beta * Phi_X.T @ Phi_X)
mN = beta * SN @ Phi_X.T @ y
mean_weight = Phi_s @ mN
cov_weight = Phi_s @ SN @ Phi_s.T

# Check the same predictions computed through covariances
mean_kernel, cov_kernel = predict_from_kernel(k_quadratic, x, y, beta, x_grid)
print("Largest mean difference", np.max(np.abs(mean_weight - mean_kernel)))
print("Largest covariance difference", np.max(np.abs(cov_weight - cov_kernel)))

Largest mean difference 1.4480860954790842e-11
Largest covariance difference 3.5517144780783383e-12

Did this rewrite change our model? Compare the sizes of the systems solved in weight space and function space. What would we need to evaluate in order to make predictions without constructing the $D$ features?

Solution to (b)

Set $\mathbf A=\mathbf S_0$, $\mathbf B=\boldsymbol\Phi$, and $\mathbf C=\beta^{-1}\mathbf I_N$. The identity gives

$$ \mathbf m_N=\mathbf S_0\boldsymbol\Phi^{\mathsf T} (\boldsymbol\Phi\mathbf S_0\boldsymbol\Phi^{\mathsf T}+\beta^{-1}\mathbf I_N)^{-1}\mathbf y. $$

Multiplying by $\boldsymbol\Phi_\star$ and identifying the covariance matrices yields

$$ \mathbb E[\mathbf f_\star\mid\mathbf y] =\mathbf K_{\star X}(\mathbf K_{XX}+\beta^{-1}\mathbf I_N)^{-1}\mathbf y. $$

The model has not changed. Its prior, observations, and predictions are the same. We have replaced a $D\times D$ system involving features by an $N\times N$ system involving covariances between observations. If $k(x,x')$ can be evaluated directly, we can predict without constructing the features. This is the kernel trick. It need not be faster when $D$ is small, but it permits models whose feature vectors would be too large to construct.

The supplied covariance follows from Gaussian conditioning. In the joint distribution

$$ \begin{bmatrix}\mathbf f_\star\\\mathbf y\end{bmatrix} \sim\mathcal N\!\left(\mathbf 0, \begin{bmatrix} \mathbf K_{\star\star} & \mathbf K_{\star X}\\ \mathbf K_{\star X}^{\mathsf T} & \mathbf K_{XX}+\beta^{-1}\mathbf I_N \end{bmatrix}\right), $$

the conditional covariance is $\boldsymbol\Sigma_{11}-\boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1}\boldsymbol\Sigma_{21}$. Substitution gives the formula in the question. These are covariances of the latent function values. Predicting new noisy outputs requires adding $\beta^{-1}\mathbf I_M$.

(c) Choose a covariance without constructing its features¶

The Gaussian kernel, also called the RBF or squared exponential kernel, is

$$ k_{\mathrm{RBF}}(x,x')=\theta\exp\!\left(-\frac{(x-x')^2}{2\ell^2}\right). $$

Here $\theta$ is the prior variance at each input and $\ell$ sets the distance over which function values are strongly correlated. This kernel admits an infinite-dimensional feature representation. You may take this fact as given. Explain why part (b) still lets us compute its predictions.

Implement k_rbf and compare its predictions with the quadratic model over the wider interval below. Far from all training inputs, what happens to $k(x_\star,x_i)$ and $k(x_\star,x_\star)$ for the RBF kernel? Use the predictive formulas to explain the mean and uncertainty you see. Why does the quadratic model behave differently?

In [25]:
def k_rbf(a, b, theta=64.0, ell=1.0):
    a = np.asarray(a).reshape(-1, 1)
    b = np.asarray(b).reshape(-1, 1)
    K = ...
    K = theta * np.exp(-(a - b.T)**2 / (2 * ell**2))
    return K

def plot_posterior(ax, x_test, mean, cov, title):
    # Show pointwise uncertainty about the latent function f
    std = np.sqrt(np.maximum(np.diag(cov), 0))
    ax.fill_between(x_test, mean - 2 * std, mean + 2 * std,
                    color="tab:blue", alpha=0.2, label="Mean plus or minus 2 sd")
    ax.plot(x_test, mean, color="tab:blue", label="Posterior mean")
    ax.scatter(x, y, color="black", s=20, zorder=3, label="Noisy observations")
    for edge in [x.min(), x.max()]:
        ax.axvline(edge, color="gray", linestyle="dotted")
    ax.set_title(title)
    ax.set_xlabel("x")
    ax.set_ylabel("f(x)")

x_far = np.linspace(-10, 10, 400)
mean_quad, cov_quad = predict_from_kernel(k_quadratic, x, y, beta, x_far)
mean_rbf, cov_rbf = predict_from_kernel(k_rbf, x, y, beta, x_far)

fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)
plot_posterior(axes[0], x_far, mean_quad, cov_quad, "Quadratic covariance")
plot_posterior(axes[1], x_far, mean_rbf, cov_rbf, "Gaussian covariance")
axes[0].legend(fontsize=8)
plt.tight_layout()
plt.show()

Solution to (c)

The prediction routine only requires kernel evaluations at finitely many training and test inputs. With the Gaussian kernel, all these numbers can be computed without representing its infinite feature vector. Changing the covariance changes the model, while the same prediction routine still applies. Both the quadratic model and the RBF model are Gaussian processes.

Far from every training input, $k(x_\star,x_i)\to 0$ for the RBF kernel, but $k(x_\star,x_\star)=\theta$ remains constant. The predictive mean therefore tends to zero and the latent variance tends to $\theta=64$. The observations no longer constrain this location, so the prediction returns to the prior mean and uncertainty.

For the quadratic model, the same three coefficients determine the function everywhere. Its posterior mean continues as a global quadratic, and its variance is $\boldsymbol\phi(x_\star)^{\mathsf T}\mathbf S_N\boldsymbol\phi(x_\star)$. With these features and nonzero observation noise, the leading term is $x_\star^4$ times the posterior variance of the quadratic coefficient. The uncertainty eventually grows with distance, though it remains relatively narrow over the plotted interval.

The contrast reflects different prior assumptions. A known global polynomial relationship could justify the quadratic model. The RBF prior expresses local correlation and retains uncertainty at distant inputs. That behavior follows from its covariance, rather than from infinite dimension alone. Neither extrapolation is established as correct by the observed data.