Equivalent Kernel - Bishop Chapter 3 I've been struggling to understand the Equivalent Kernel in Bishop's Pattern Recognition and Machine Learning book. Can somebody explain the following Figure (3.10) from chapter 3.3.3?

(image taken from online slides)
The caption reads:
The equivalent kernel $k(x, x')$ for the Gaussian basis functions in Figure 3.1, shown as a plot of $x$ versus $x'$, together with three slices through this matrix corresponding to three different values of $x$. The data set used to generate this kernel comprised 200 values of $x$ equally spaced over the interval $(−1, 1)$.
I don't understand what he's trying to show in the plot, and I have no idea what the image in the right side is supposed to mean - I didn't see explanation about the right image at all...
 A: Maybe a coded example (Python) with some accompanything theory can help. The image on the right is is the equivalent kernel; which is a (symmetric) covariance matrix over the inputs. The image on the left shows, for a given value of $x$ (e.g. $x_i$), how $x_i$ covaries with all other values of $x$. From this we can see that as we move away from $x_i$, the covariance decays in a squared exponential fashion (i.e. points close to $x_i$ are more similar than those further away).
Recall the standard Bayesian linear regression problem, where $\pmb{\phi} = \phi(x)$ may represent a basis function projection of the inputs:
\begin{equation}
    \textbf{y} = f(\pmb{\phi}) + \epsilon = \pmb{\phi}^\text{T}\pmb{\beta} + \epsilon \quad , \quad \epsilon \sim \mathcal{N}(0,\sigma^2).
\end{equation}
Here $\textbf{y}$ are the training outputs, $\pmb{\phi}$ are the training inputs, $\pmb{\beta}$ are the set of regression parameters, and $\epsilon$ is i.i.d Gaussian noise. For this model one obtains a likelihood over the outputs $\textbf{y} \sim \mathcal{N}(\pmb{\phi}^\text{T}\pmb{\beta},\pmb{V})$, where $\pmb{V} = \sigma^2\pmb{I}$.
Now, for a zero-mean prior over the parameters $\pmb{\beta} \sim \mathcal{N}(0,\pmb{\Sigma})$, the posterior over the parameters, with mean $\tilde{\pmb{\mu}}$ and covariance $\tilde{\pmb{\Sigma}}$, is given by:
\begin{equation}
    \tilde{\pmb{\mu}} = (\pmb{\phi}\pmb{V}^{-1}\pmb{\phi}^\text{T} + \pmb{\Sigma}^{-1})^{-1}\pmb{\phi}\pmb{V}^{-1}\textbf{y} \\
\tilde{\pmb{\Sigma}} = (\pmb{\phi}\pmb{V}^{-1}\pmb{\phi}^\text{T} + \pmb{\Sigma}^{-1})^{-1}.
\end{equation}
If we now revist our model $f(\pmb{\phi})$ with our derived estimate of the posterior mean, we arrive at $f(\pmb{\phi}) = \pmb{\phi}^\text{T}\tilde{\pmb{\mu}}$. Explicitly, $\pmb{\phi}^\text{T}\tilde{\pmb{\mu}}$ is the posterior predictive mean, evaluated at the locations $\pmb{\phi}$ (in this case evaluated at the location of the training inputs). If we expand the notation we arrive at:
\begin{equation}
    \pmb{\phi}^\text{T}\tilde{\pmb{\mu}} = \pmb{\phi}^\text{T}(\pmb{\phi}\pmb{V}^{-1}\pmb{\phi}^\text{T} + \pmb{\Sigma}^{-1})^{-1}\pmb{\phi}\pmb{V}^{-1}\textbf{y}.
\end{equation}
Notice here that $\pmb{\phi}^\text{T}(\pmb{\phi}\pmb{V}^{-1}\pmb{\phi}^\text{T} + \pmb{\Sigma}^{-1})^{-1}\pmb{\phi}$ is actually a valid covariance matrix, i.e. the equivalent kernel. Hence $k(\pmb{\phi},\pmb{\phi}') = \pmb{\phi}^\text{T}(\pmb{\phi}\pmb{V}^{-1}\pmb{\phi}^\text{T} + \pmb{\Sigma}^{-1})^{-1}\pmb{\phi}$. This means that the posterior predictive mean can actually be computed as a linear combination of the outputs, for each column of $k(\pmb{\phi},\pmb{\phi}')$:
\begin{equation}
    \pmb{\phi}^\text{T}\tilde{\pmb{\mu}} = \sigma^{-2}\sum_{i=1}^{n}k(\pmb{\phi},\pmb{\phi}_i)\text{y}_i.
\end{equation}
Now for the coded example:
import numpy as np
from numpy.linalg import multi_dot as mdot
from numpy.linalg import inv as inv
import matplotlib.pyplot as plt

np.random.seed(2)

def f(x):
    np.random.seed(2)
    return np.random.multivariate_normal(np.zeros(n),np.exp(-squareform(pdist(np.atleast_2d(x).T/0.3,'sqeuclidean'))))

#Generate some synthetic data
n = 100 #number of observations
sigma2 = 0.5 #noise variance
x = np.linspace(-1,1,n) #inputs
phi = np.array([np.ones(n),x,x**2,x**3,x**4,x**5,x**6]) #polynomial basis function projection
N = phi.shape[0] #number of parameters
fx = f(x) #true function values drawn from a Gaussian process
y = fx + np.random.normal(0,sigma2,n) #add random Gaussian noise to produce outputs

Vi = np.eye(n)*(1/sigma2) #inverse covariance of the likelihood distribution
Sigmai = np.eye(N)*(1/100) #inverse prior covariance over the regression parameters

#Generate mean of predictive distribution
phi_mu_tilde = mdot([phi.T,inv(mdot([phi,Vi,phi.T]) + Sigmai),phi,Vi,y]) #posterior predictive distribution at phi

kx = mdot([phi.T,inv(mdot([phi,Vi,phi.T]) + Sigmai),phi,Vi]) #equivalent kernel

#plot kernel and extract profiles 
plt.imshow(kx,cmap='jet',vmin=0,vmax=0.1)
plt.hlines(25,100,0,color='k',alpha=0.25)
plt.hlines(50,100,0,color='k',alpha=0.5)
plt.hlines(75,100,0,color='k',alpha=1)
plt.xlim(n,0)
plt.ylim(n,0)
plt.show()

plt.plot(x,kx[:,25],color='k',lw=2,alpha=1)
plt.plot(x,kx[:,50],color='k',lw=2,alpha=0.5)
plt.plot(x,kx[:,75],color='k',lw=2,alpha=0.25)
plt.show()

#show that linear combination of outputs is equivalent to mean of predictive distribution
linear_comb = [kx[:,i] * y[i] for i in range(n)] #linear combination of the outputs
print(np.isclose(phi_mu_tilde,np.sum(linear_comb,0),atol=1e-15)) #check equivalence



