# Learning a Gaussian Process from function observations (not GP regression)

Suppose we have a set of observations, where each observation represents a function. For example, our set is $$\{f_1, f_2, ..., f_n\}$$ where each $$f_i = \{(x_1, y_1), (x_2, y_2), ..., (x_{p_i}, y_{p_i})\}$$ representing evaluations of that function.

Now I want to learn a distribution over functions from these function observations, for example a Gaussian Process. What methods can you use for this?

Please note this is not GP regression, where you get a posterior distribution over functions representing your uncertainty over a single true function. I'm instead interested in maintaining a distribution over functions from which you sample genuinely different functions.

As a concrete example, each function observation $$f_i$$ could represent temperature over the course of a year, I have records from many years, and I'm interested in learning about what temperature profiles generally look like (i.e. warmer in summer, colder in winter), as well as what variations between years looks like.

• As I understand you are interested in a GP with both a continuous input $x$ and a categorical input, say $u$. Your observations $Y(x_i,\,u_i)$ are made at a finite set of design points. For instance $x$ can be an time within the year (in practice a day number between $1$ and $365$) while the year is categorical or ordinal. In this setting you can not extrapolate at an unobserved level of $u$. – Yves Mar 9 '19 at 8:03
• @Yves Thanks for the reply. I think my question is slightly different from what you suggest. A GP is distribution over functions, so we can take draws from this distribution, $f_i \sim GP(\mu, k)$, where $f_i$ is a finite-dimensional representation of the function (e.g. evaluations at ${(x_1, y_i), ..., (x_p, y_p)}$. For example, each draw could be an year, "i.i.d." from the GP. Essentially we should be able to learn about the GP's mean and covariance functions this way. The difference with your interpretation is that each year is just an identical draw - there is no additional variable $u$. – BayesCruncher Mar 10 '19 at 1:09
• So you have several paths of a $\text{GP}(\mu,\,k)$; you want to infer on $\mu$ and $k$ or predict/smooth within a specific path. This enters the previous framework with $u$ corresponding to a path. Indeed, a special covariance structure for this setting is $\text{Cov}([x,\,u],\,[x',\,u']) = k_{\text{cont}}(x,\,x') k_{\text{cat}}(u,\,u')$ where $k_{\text{cont}}(x, \, x')$ is a continuous kernel and $k_{\text{cat}}(u,\,u')$ is a categorical kernel, which essentially is a covariance matrix $[\Gamma_{u,u'}]_{u,u'}$. For i.i.d. paths take an identity matrix. – Yves Mar 10 '19 at 9:14
• @Yves Yes, that's exactly my framework. As I'm making the i.i.d. assumption $k_{cat}(u, u')$ is the identity matrix as you say, so it does not need to be inferred. In this case do you know how I may go about inferring the mean function $\mu(x)$ and covariance function $k_{cont}(x, x')$? – BayesCruncher Mar 10 '19 at 22:01
• You need use a software that can cope with a product kernel as described and makes maximum-likelihood or Bayesian inference. The computational task depends on the design, the kernel $k_{\text{cont}}$ and the dimension $d$ of the continuous input $x$. The simple case is $d = 1$, a tensor design (same $x_i$ for all paths) and a semi-Markov kernel such as Matérn with shape $1/2$, $3/2$, $5/2$, ... But in general, the computational cost is $O(N^3)$ where $N$ is the number of sample points. – Yves Mar 11 '19 at 6:42

First off, we need to specify a parametric form for both the mean function $$\mu$$ (e.g. linear function $$\alpha + \beta^Tx$$) and the covariance $$k$$ (e.g. RBF or Matérn). These parameters $$\alpha$$, $$\beta$$, $$\sigma^2$$ etc. are often known as hyperparameters in GP regression.
In the general case when the function paths are not i.i.d., we can index each with a variable $$u$$. Then the solution is to simply perform GP regression with the mean and covariance functions taking the forms $$\mu(x, u)$$ and $$k((x, u), (x', u'))$$, and then optimize the hyperparameters using e.g. maximum marginal likelihood or Bayesian inference. However, the complexity for the regression is $$O(N^3)$$, where N is the total number of points $$(x, u)$$, i.e. $$N = \sum_{i=1}^n p_i$$.
However, in the i.i.d. case, we can simply run a GP regression on each of the $$n$$ datasets $$f_i$$ individually, optimizing or Bayesian updating the same hyperparameters for all of these. In this case the complexity is (assuming the number of datapoints per dataset $$p_i = p$$ is fixed) $$O(np^3)$$ instead of $$O(n^3p^3)$$.