> For the complete documentation index, see [llms.txt](https://zavalab.gitbook.io/statschemebook/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://zavalab.gitbook.io/statschemebook/chapter-6-statistical-data-analysis-and-learning/01-classification-models/02-kernel-models.md).
# Kernel Models
* Kernel models are generalizations of predictive models that we have explored.
* These models have surprising properties:
* Require a handful of parameters to make predictions (regardless of number of inputs)
* Capture complex nonlinearities (without need to postulate input-output model)
* Concepts have implications in techniques that we have discussed.
* For instance, kernel concepts can be used to capture nonlinearities in PCA.
* Key principle of kernel methods is that covariance of our data can be represented as a matrix function (rather than a fixed matrix).
* Such function (a.k.a. kernel function) captures complex nonlinear relationships present in our data.
* To explain origin of kernel methods, we explore their use in regression.
* A flexible approach to construct models consists of *mixing* different model types.
* Assume that $$Y$$ and $$X=(X\_1,...,X\_n)$$ follow relationship:
$$
Y=\sum\_{j=1}^m\theta\_j \phi\_j(X)+\epsilon
$$
where $$\phi\_j(X),; j=1,...,m$$ is collection of basic models (known as basis functions).
* Here, output and input have continuous domains $$\mathcal{D}\_Y\subseteq \mathbb{R}$$ and $$\mathcal{D}\_X\subseteq \mathbb{R}^n$$.
* Predictive model is mixture $$m(X,\theta)=\sum\_{j=1}^m\theta\_j \phi\_j(X)$$.
* Define basis vector $$\phi(X)=(\phi\_1(X),...,\phi\_m(X))$$ and write $$m(X,\theta)=\theta^T\phi(X)$$.
* Basis $$\phi(X)$$ can be linear or nonlinear (e.g., polynomial, sigmoidal, log, exponential).
* Parameters $$\theta\in \mathbb{R}^m$$ are mixing coeffs for basis functions.
* Determine $$\hat{\theta}$$ by solving:
$$
\min\_{\theta}; \frac{1}{2}\sum\_{\omega \in \mathcal{S}}(y\_\omega-\theta^T\phi(x\_\omega))^2+\lambda \theta^T\theta
$$
* This is Bayesian estimation with noise $$\epsilon \sim \mathcal{N}(0,\sigma^2\mathbf{I})$$ and prior $$\theta\sim \mathcal{N}(0,\lambda \mathbf{I})$$.
* Express problem in compact form:
$$
\min\_{\theta}; \frac{1}{2}(\mathbf{y}-\mathbf{\boldsymbol{\Phi}}\theta)^T(\mathbf{y}-\boldsymbol{\Phi}\theta)+\lambda\theta^T\theta
$$
where $$\boldsymbol{\Phi}\in \mathbb{R}^{S\times m}$$ is data matrix with entries $$\boldsymbol{\Phi}*{\omega,j}=\phi*{j}(x\_\omega)$$ and $$\mathbf{y}$$ is output vector.
* Data matrix $$\boldsymbol{\Phi}$$ is transformation of original input data.
* If basis functions are linear then $$\boldsymbol{\Phi}=\mathbf{X}$$ (get a standard linear estimation prob).
* We define objective as $$S(\theta)$$; optimal estimate must satisfy $$\nabla\_\theta S(\theta)=0$$:
$$
\theta=\frac{1}{\lambda}\boldsymbol{\Phi}^T(\mathbf{y}-\boldsymbol{\Phi}\theta).
$$
* Define residuals $$\mathbf{r}=\frac{1}{\lambda}(\mathbf{y}-\boldsymbol{\Phi}\theta)$$ and thus $$\theta=\boldsymbol{\Phi}^T\mathbf{r}$$.
* By substituting estimate in $$S(\theta)$$:
$$
S(\theta)= \frac{1}{2}\mathbf{r}^T\boldsymbol{\Phi}\boldsymbol{\Phi}^T\boldsymbol{\Phi}\boldsymbol{\Phi}^T\mathbf{r}-\mathbf{r}^T\boldsymbol{\Phi}\boldsymbol{\Phi}^T\mathbf{y}+\frac{1}{2}\mathbf{y}^T\mathbf{y}+\frac{\lambda}{2}\mathbf{r}^T\boldsymbol{\Phi}\boldsymbol{\Phi}^T\mathbf{r}.
$$
* Data matrix $$\boldsymbol{\Phi}$$ always appears as $$\boldsymbol{\Phi}\boldsymbol{\Phi}^T$$; we define $$\mathbf{K}=\boldsymbol{\Phi}\boldsymbol{\Phi}^T\in \mathbb{R}^{S\times S}$$ and thus:
$$
S(\theta)= \frac{1}{2}\mathbf{r}^T\mathbf{K}\mathbf{K}\mathbf{r}-\mathbf{r}^T\mathbf{K}\mathbf{y}+\frac{1}{2}\mathbf{y}^T\mathbf{y}+\frac{\lambda}{2}\mathbf{r}^T\mathbf{K}\mathbf{r}
$$
* $$S(\theta)$$ can be defined in terms of $$\mathbf{r}$$ (params have been eliminated).
* This suggests an alternative strategy to find $$\hat{\theta}$$:
* Find optimal residual estimates $$\hat{\mathbf{r}}$$ by solving:
$$
\min\_{\mathbf{r}} \frac{1}{2}\mathbf{r}^T\mathbf{K}\mathbf{K}\mathbf{r}-\mathbf{r}^T\mathbf{K}\mathbf{y}+\frac{1}{2}\mathbf{y}^T\mathbf{y}+\frac{\lambda}{2}\mathbf{r}^T\mathbf{K}\mathbf{r}
$$
* Recover $$\hat{\theta}=\boldsymbol{\Phi}^T\hat{\mathbf{r}}$$.
* It turns out, however, that $$\hat{\theta}$$ are *not needed at all* to make predictions.
* Specifically, optimal model prediction is:
$$
\begin{aligned}
\hat{\mathbf{y}}&= \boldsymbol{\Phi}\hat{\theta}\\
&=\boldsymbol{\Phi}\boldsymbol{\Phi}^T\hat{\mathbf{r}}\\
&=\mathbf{K}(\mathbf{K}+\lambda I)^{-1}\mathbf{y}.
\end{aligned}
$$
* Consequently, prediction only depends on input data $$\mathbf{K}$$ and output data $$\mathbf{y}$$.
* Params $$\hat{\theta}$$ are *intermediary* variables.
* Kernel matrix $$\mathbf{K}$$ captures interactions of input variables $$X$$.
* Specifically, given $$x\_\omega,; \omega \in \mathcal{S}$$, kernel matrix is:
$$
\mathbf{K}\_{i,j}=k(x\_i,x\_j),; i,j\in \mathcal{S}
$$
where $$k(x\_i,x\_j)$$ is known as *kernel function*.
* Kernel generalizes strategies discussed; e.g., linear estimation corresponds to:
$$
k(x\_i,x\_j)=x\_i^Tx\_j,; i,j\in \mathcal{S}
$$
* Case of mixtures of basis functions corresponds to:
$$
k(x\_i,x\_j)=\phi(x\_i)^T\phi(x\_j),; i,j\in \mathcal{S}
$$
* Kernel function captures nonlinear interactions in input data.
* A widely used kernel is radial basis function (RBF):
$$
k(x\_i,x\_j)=\textrm{exp}(-\gamma|x\_i-x\_j|^2),; i,j\in \mathcal{S}
$$
where $$\gamma\in \mathbb{R}\_+$$ is a param; this is also known as Gaussian kernel.
* Gaussian kernel captures nonlinear interactions in input data.
* To see this, consider scalar $$x\_\omega$$ and notice that:
$$
\begin{aligned}
&\exp(-\gamma(x\_i-x\_j)^2)\\
&=\exp(-\gamma x\_i^2+2\gamma x\_ix\_j-\gamma x\_j^2)\\
&=\exp(-\gamma x\_i^2-\gamma x\_j^2)\left(1+\frac{2\gamma x\_ix\_j}{1, !}+\frac{(2\gamma x\_ix\_j)^2}{2,!}+...\right)\\
&=\exp(-\gamma x\_i^2)\exp(-\gamma x\_j^2)
\left(1\cdot1 +\left(x\_i\sqrt{2\gamma/1,!}\right)\cdot\left(x\_j\sqrt{2\gamma/1,!}\right)\right.\\
&\left.\qquad\qquad+\left(x\_i^2\sqrt{(2\gamma)^2/2,!}\right)\cdot \left(x\_j^2\sqrt{(2\gamma)^2/2,!}\right)+...\right)
\end{aligned}
$$
* Gaussian kernel is a separable function and can be written as:
$$
\begin{aligned}
\exp(-\gamma(x\_i-x\_j)^2)&=\phi(x\_i)^T\phi(x\_j)
\end{aligned}
$$
* Where $$\phi(x)$$ is an *infinite* collection of polynomial basis functions:
$$
\phi(x)=e^{-\gamma x^2}(1,; x\sqrt{2\gamma/1,!},; x^2\sqrt{(2\gamma)^2/2,!},; x^3\sqrt{(2\gamma)^3/3,!},...).
$$
* Here, a *single* param $$\gamma$$ controls coefficients of basis functions.
* In kernel methods, we do not need to postulate input-output model (e.g., $$Y=\theta^T\phi(X)$$) and estimate its params $$\theta$$).
* Instead, we postulate kernel model and estimate its params $$\gamma$$.
* Kernel can be interpreted as a generating function of basis functions and requires few params (typically a handful).
* This is key advantage over standard estimation because overparameterization is avoided.
* Kernels can capture highly complex nonlinear behavior.
## Example: Kernel Estimation for Hougen-Watson
* Consider Hougen-Watson model:
$$
Y=\frac{(\theta\_0X\_2-X\_3/\theta\_4)}{(1+\theta\_1X\_1+\theta\_2X\_2+\theta\_3X\_3)}
$$
* This model might not fit data if mechanism proposed is not correct.
* Instead of using this model, we use data $$(x\_\omega,y\_\omega)$$ to construct a kernel model.
* We use $$k(x\_\omega,x\_{\omega'})=\exp(-\gamma|x\_\omega-x\_{\omega'}|^2)$$ with $$\gamma=0.1,0.001,0.0001$$.
#### Python implementation
{% hint style="success" %}
Notebook [`ch6_hougen_watson_kernel.ipynb`](https://colab.research.google.com/github/vzavala/StatsBook/blob/master/Python/ch6_hougen_watson_kernel.ipynb) — open it in Colab and run the cells alongside this section.
{% endhint %}
Kernel model for the Hougen-Watson data (Victor Z and Ugochukwu Ikegwu, UW-Madison, 2024).
| Slides | Python | Meaning |
| -------------------------------------------------------------------------------- | -------------------- | ------------------------------------ |
| $$x\_\omega$$, $$y\_\omega$$ | `X`, `y` | partial pressures and reaction rates |
| $$k(x\_\omega,x\_{\omega'})=\exp(-\gamma\lVert x\_\omega-x\_{\omega'}\rVert^2)$$ | `rbf(xi, xj, gamma)` | RBF kernel |
| $$\mathbf{K}$$, $$\lambda$$ | `K`, `lam` | kernel matrix and regularization |
| $$\hat{\mathbf{y}}=\mathbf{K}(\mathbf{K}+\lambda I)^{-1}\mathbf{y}$$ | `y_hat` | kernel prediction |
**Note:** This notebook reproduces Example 6.23 (Section 6.5) of *Statistics for Chemical Engineers*: a kernel model of the Hougen-Watson reaction rate for different kernel parameters gamma.
```python
import numpy as np
import matplotlib.pyplot as plt
data = np.loadtxt(path+'Data/hougenwatson.dat')
X, y = data[:, :3], data[:, -1]
S = len(y)
lam = 0.01
def rbf(xi, xj, gamma):
return np.exp(-gamma * np.sum((xi - xj)**2))
fig, axs = plt.subplots(1, 3, figsize=(15, 4.5))
for ax, gamma in zip(axs, (0.1, 0.001, 0.0001)):
K = np.array([[rbf(X[i], X[j], gamma) for j in range(S)] for i in range(S)])
y_hat = K @ np.linalg.solve(K + lam * np.eye(S), y)
ax.plot(y_hat, y, 'ko', markerfacecolor='w')
ax.plot([0, 15], [0, 15], 'k--')
ax.set_title(rf'$\gamma={gamma}$: $\|y-\hat{{y}}\|$ = {np.linalg.norm(y - y_hat):.2f}')
ax.set_xlabel(r'$\hat{y}$'); ax.set_ylabel('$y$')
ax.grid(True)
plt.tight_layout()
plt.show()
```
**Note:** As in the book, the fit is rather insensitive to gamma, but it degrades for gamma = 0.0001: the kernel parameter must be tuned.
