Kernels in Gaussian Processes I am trying to understand intuitively how a kernel works in a Gaussian process. I know that GP are distributions over functions, in short you have the model $y = f(x)+\epsilon$ and the $f(x)$ follows a Gaussian process
$GP(m(x),K(x,x^{\prime}))$.   
After deriving the conditional distribution we arrive at the predictive equations (with noise):  
$\hat{f}=K(x,x^{\prime})(K(x,x)+\sigma^{2}I)^{-1}y$ for the mean  
and  
$cov(f)=K(x^{\prime},x^{\prime})-K(x^{\prime},x)(K(x,x)+\sigma^{2}I)^{-1}K(x,x^{\prime})$ for the covariance.  
The covariance kernel is the SE (Squared Exponential):  
$k_{f}(x_{i},x_{j}) = \sigma^{2}exp(-\frac{1}{2\ell^{2}}\sum_{j=1}^{q}(x_{i,j}-x_{k,j})^{2})$
I know that a kernel, when we have a vector with inputs $x$ and size $(n)$, has dimensions $(n\times n)$. But how it works in multidimensional environment?
For example, the covariance above $cov(f)$ with one dimensional inputs $(n)$ and for a single test point $(p = 1)$ will have the form:
$(p \times p) - (p \times n)(n \times n)(n \times p)$ which gives a scalar (the variance). 
I can't comprehend how the above equation will work when we accept as input a matrix $(n \times d)$ size instead. 
 A: Notation / setting
We are considering a GP regression model:
\begin{equation}
y_i = f(x_i) + \epsilon_i
\end{equation}
where $y_i\in \mathbb{R}$,$x_i \in \mathbb{R}^d$, $f$ a Gaussian process (whose realizations are functions $f:\mathbb{R}^d\rightarrow \mathbb{R}$), 
\begin{equation}
f \sim \mathrm{GP}(m(x_i), \kappa(x_i,x_j)).
\end{equation}
$n$ datapoints $(y_1,x_1), (y_2,x_2), (y_3,x_3),\ldots, (y_n,x_n)$ are given. (I  use $\kappa$ to distinguish the function from the matrices $K(\cdot,\cdot)$ that contain values of $\kappa$ evaluated at certain points. The question denotes both by $K$)
How to handle $d$-dimensional inputs
The question covers computing the posterior predictive distribution for a test point (or $p$ test points) in the case $d=1$ and asks how to extend to the general $d=2,3,\ldots$.
Answer: nothing changes and the formulas from the one-dimensional case work as well in this case. Note that $m$ is then a function from $\mathbb{R}^d$ to $\mathbb{R}$ and $\kappa$ a function from $\mathbb{R}^d \times \mathbb{R}^d$ to $\mathbb{R}$.
So, for example the matrix denoted by $K(x,x)$ in the question is a $n\times n$ matrix for which $K(x,x)_{i,j} = \kappa(x_i, x_j)$  ($x_i$ and $x_j$ are $d$-dimensional but since $\kappa$ maps two $d$-dimensional vectors to a scalar, $\kappa(x_i, x_j)$ is a scalar. Similarly for $K(x,x')$ and $K(x',x')$ where $x'$ are the test points.
Thus, the dimensions of the matrices in the predictive covariance equation are $(p\times p) - (p \times n)\,(n \times n)\,(n \times p)$ independent of whether the elements of the matrices are obtained by evaluating a function $\kappa(\cdot,\cdot)$ whose arguments are $1$-dimensional of a function $\kappa(\cdot, \cdot)$ whose arguments are $d$-dimensional. In fact, the inputs could even be in some space other than $\mathbb{R}^d$ (such as if we have a categorical predictor) as long as a positive-definite covariance function can be defined.
An extra remark about the SE kernel appearing in the question
The question mentions the SE kernel 
\begin{equation}
k_{f}(x_{i},x_{j}) = \sigma^{2}\exp\!\Big(-\frac{1}{2\ell^{2}}\sum_{j=1}^{q}(x_{i,j}-x_{k,j})^{2}\Big)
\end{equation}
Note that this is already a function from $\mathbb{R}^q \times \mathbb{R}^q$ to $\mathbb{R}$ (with scalar inputs there would be no "$x_{i,j}$ and $x_{k,j}$" for different values of $j$. And $q$ should be $d$ if $d$ is the dimension of inputs.
Optionally, the length scale $\ell$ could be made different for each input dimension as $\ell_{j}$, such that the term $\frac{1}{2\ell_{j}^{2}}$ is instead placed inside the summation
