Does gaussian process regression model interactions? Suppose there are two feature variables and one target variable. We want to predict the target variable based on a Gaussian Process Regression model. Does this model (e.g., the GPR model of the type provided in scikit-learn) automatically include interactions between features? Or new features have to be created for this?
Clarification: I am not talking about whether the sklearn API would add new interaction features in the preprocessing step, but whether GPR models can model interactions automatically like RandomForests and Deep Neural Nets.
 A: Before discussing GPs I want to be more precise about what we're looking for and what it means for a model to be able to capture interactions. For the notation I'll consider the inputs as $x \in \mathcal X\subseteq \mathbb R^p$ and the target is $y \in \mathbb R$.
Suppose we are estimating $y$ with a function $\hat f(x)$. If there are no interactions then this means that the effects on $\hat y$ for varying one coordinate of $x$ should not depend on the other coordinates. I'll formalize this by saying that a model $f$ does not capture interactions if and only if it can be written in additive form as
$$
\hat f(x) = \hat\beta_0 + \sum_{i=1}^p \hat\beta_i h_i(x_i).
$$
A consequence of this is that, if $\hat f$ is differentiable, then I can check for additivity by checking if
$$
\nabla \hat f(x) = (h_1(x_1), \dots, h_p(x_p))^T
$$
for some functions $h_1,\dots,h_p$ and all $x \in \mathcal X$.
I'll check this with three examples: (1) linear regression with no interaction terms; (2) linear regression with an interaction; (3) a decision tree. I'll then apply these insights to GPs.
1. Linear regression with no interactions
We have $\hat f(x) = x^T\hat\beta$ so this is additive with $h_i = \text{id}$ and
$$
\nabla \hat f(x) = \hat\beta.
$$
2. Linear regression with interactions
I'll take $p=2$ and consider $\hat f(x) = \hat \beta_0 + \hat \beta_1 x_1 + \hat \beta_2 x_2 + \hat \beta_{12}x_1x_2$. Now
$$
\nabla \hat f(x) = \begin{bmatrix}\hat\beta_1 + \hat\beta_{12}x_2 \\ \hat\beta_2 + \hat\beta_{12}x_1\end{bmatrix}
$$
so this does model interactions, because now the effect of an increase in one variable on the target depends on the value of the other variable.
3. a single decision tree
I'll again take $p=2$. A decision tree is of the form
$$
\hat f(x) = \sum_{k=1}^m \hat \beta_k I(x \in C_k)
$$
where the $C_k$ are axis-parallel rectangles. Whether or not this is additive depends on the $C_k$.
For this example I'll suppose we have a cut at $x_1 = 1$ and another at $x_2 = -2$ so the four regions are
$$
C_1 = C_{++} = [1, \infty)\times [-2, \infty) \\
C_2 = C_{+-} = [1, \infty) \times (-\infty, -2) \\
C_3 = C_{-+} = (-\infty, 1)\times [-2, \infty) \\
C_4 = C_{--} = (-\infty, 1)\times(-\infty, -2).
$$
Suppose initially $x_2 = 2.5$ so if $x_1 < 1$ we're in $C_{+-}$ while if $x_1 \geq 1$ then $x \in C_{++}$. Near the line $\{x_1 = 1\}$ we'll find that a small increase in $x_1$ can lead to the predictions changing from $\hat\beta_{-+}$ to $\hat\beta_{++}$. But if instead $x_2 = -2.5$, say, then we'll change from $\hat\beta_{--}$ to $\hat\beta_{+-}$. This shows that the change in predictions as we vary one input depends on the value of the other input, so we are capturing some interactions even though we can't use the gradient test.
If instead I had $C_1 = [1, \infty)\times\mathbb R$ and $C_2 = \mathbb R\times [-2, \infty)$ then the model could be written in additive form so it does not capture interactions.
Applying to GPs
$\newcommand{\e}{\varepsilon}\newcommand{\y}{\mathbf y}\newcommand{\0}{\mathbf 0}\newcommand{\kl}{k_{\text L}}\newcommand{\kse}{k_{\text{SE}}}\newcommand{\fl}{f_{\text L}}\newcommand{\fse}{f_{\text{SE}}}$I'll now do this with GPs. As @user20160 said in their comment, it depends on the kernel/covariance function we use.
We have $y = f(x) + \e(x)$ with $f \sim \mathcal {GP}(0, k)$ and $\e(x)\stackrel{\text{iid}}\sim \mathcal N(0, \alpha)$ for $\alpha > 0$. I'll use $\kl(x,x') = x^Tx'$ as the linear kernel and $\kse(x,x') = \exp(-\gamma\|x-x'\|^2)$ as the squared exponential kernel. The corresponding posterior mean functions will be $\hat \fl$ and $\hat \fse$ respectively.
We'll have
$$
\hat \fl(x) = x^TX^T(XX^T + \alpha I)^{-1}\y
$$
and
$$
\hat\fse(x) = \begin{bmatrix}\exp(-\gamma\|x-x_1\|^2) \\ \vdots \\ \exp(-\gamma\|x-x_n\|^2)\end{bmatrix}^T (K + \alpha I)^{-1}\y.
$$
For $\hat\fl$ we can use the fact that $X^T(XX^T + \alpha I_n)^{-1} = (X^TX + \alpha I_p)^{-1}X^T$ to rewrite this as
$$
\hat \fl (x) = x^T(X^TX + \alpha I_p)^{-1}X^T\y = x^T\hat\beta_\alpha
$$
where $\hat\beta_\alpha$ is the result of a ridge regression of $\y$ on $X$. This means $\nabla \hat \fl(x) = \hat\beta_\alpha$ and this matches the requirement for $\hat\fl$ to not capture interactions, and indeed we knew that already at this point since it's just a ridge regression.
But for $\hat\fse$ you can check that the gradient is not so simple and the effect of varying one input on the output depends on where we are in $\mathcal X$, so interactions are captured.
