Suppose data sample is $D = (X, \mathbf{y}) = \{\mathbf{x}_i, y_i = y(x_i)\}_{i = 1}^N$.
Also suppose, that we have a covariance function $k(\mathbf{x}_1, \mathbf{x}_2)$ and zero mean specified for a Gussian process. Distribution for a new point $\mathbf{x}$ will be Gaussian with mean $$m(\mathbf{x}) = \mathbf{k} K^{-1} \mathbf{y}$$ and variance $$V(\mathbf{x}) = k(\mathbf{x}, \mathbf{x}) - \mathbf{k} K^{-1} \mathbf{k}^T.$$ Vector $\mathbf{k} = \{k(\mathbf{x}, \mathbf{x}_1), \ldots, k(\mathbf{x}, \mathbf{x}_N)\}$ is a vector of covariances, matrix $K = \{k(\mathbf{x}_i, \mathbf{x}_j)\}_{i, j = 1}^N$ is a matrix of sample covariances. In case we make prediction using mean value of posterior distribution for sample interpolation property holds. Really,
$$m(X) = K K^{-1} \mathbf{y} = \mathbf{y}.$$
But, it isn't the case if we use regularization i.e. incorporate white noise term. in this case covariance matrix for sample has form $K + \sigma I$, but for covariances with real function values we have covariance matrix $K$, and posterior mean is
$$
m(X) = K (K + \sigma I)^{-1} \mathbf{y} \neq \mathbf{y}.
$$
In addition, regularization makes problem more computationally stable.
Choosing noise variance $\sigma$ we can select if we want interpolation ($\sigma = 0$) or we want to handle noisy observations ($\sigma$ is big).
Also, the Gaussian processes regression is local method because variance of predictions grows with distance to learning sample, but we can select appropriate covariance function $k$ and handle more complex problems, than with RBF. Another nice property is small number of parameters. Usually it equals $O(n)$, where $n$ is data dimension.
