# Gaussian process with boundaries on unobserved variables

I'm trying to use Gaussian process regression of which I only have basic knowledge and the problem I have to deal with has a bit of a twist relative to the natural set-up, so I was wondering if there are any results known for the type of problem I face.

I'm interested in the posterior distribution of a set of unobserved $y_i$ associated to some $x_i$ for $i$ in some set $I$. The data that I have available is of the form:

1. $(x_i,lb_i)$ for $i \in L \subset I$ where $lb_i$ is a lower bound for the unobserved $y_i$ to $x_i$.
2. $(x_i,ub_i)$ for $i \in U \subset I$ where $ub_i$ is an upper bound for the unobserved $y_i$ to $x_i$.

It can be assumed that there won't be any conflicting conditions, and that $L \neq U$, but also that $L \cup U \neq \emptyset$ and that $E[y_i] = 0$

• What do you mean with "I only have unobserved points and for some of them I have upper and lower bounds"? I am not sure if fully understand which information you have. Anyhow, it seems you might want to look at my answer here, I think (if I understand correctly your question) it covers a similar case. May 12, 2016 at 19:14
• I have edited the question to be (hopefully!) more clear. The link provided is somewhat useful (in the sense that it seems similar, I know to little about GP to directly understand the implications of assocaited observation noise. The difference with the answered question is that in some cases I have intervals available, but in some cases only one side of the interval is finite, so the technique of starting with the average point doesn't work
– Bram
May 12, 2016 at 19:31
• For those points which have finite lower and upper bound, can you say anything about the distribution of true value within those bounds? Your original version seemed to say there might be some small probability of truth being outside the bounds. May 12, 2016 at 19:41
• Much better now. Yes, for your problem the trick I suggested for the other question would not work. For the open-ended intervals I don't see an obvious way to link to a Gaussian likelihood, without making other assumptions (e.g., the typical width of the interval for an unobserved bound). If you are not willing to make further assumptions, you may have to deal with a pretty nasty non-Gaussian and observation-dependent likelihood, but that seems it's going to be quite complicated. May 12, 2016 at 19:41
• Thanks for the pointers both, very helpful! can I mark the questions answered based on those?
– Bram
May 13, 2016 at 5:33

Recapping and expanding on my discussion in the comments.

If you had $(lb_i, ub_i)$ for every $i$, one approach would be the one I recommended to this other question, that is to have a single GP with observations $\tilde{y}_i = \frac{lb_i + ub_i}{2}$ and observation-dependent noise proportional to $ub_i - lb_i$.

Since in your data $lb_i$ or $ub_i$ are not available for all data points, the alternative I can think of is to fit two GPs, $G_l$ and $G_u$, respectively to $(x_i,lb_i)$ and $(x_i,ub_i)$. This will allow you to generate predictions for the bounds at unobserved points. We are not explicitly requiring $G_l < G_u$ in training but we can add the constraint in prediction.

In particular, for a given point $x$: $$p(y|x) \propto \int_{-\infty}^{\infty} \mathcal{N}(l| \mu_l(x),\sigma^2_l(x)) \left\{\int_{l}^\infty \mathcal{N}(u| \mu_u(x),\sigma^2_u(x)) \frac{\left[l \le y \le u \right]}{u - l} du \right\} dl$$ where $[\cdot]$ denotes Iverson's bracket, and $\mu_l, \sigma^2_l$ ($\mu_u, \sigma^2_u$) are the latent mean and variance of $G_l(x)$ (resp. $G_u(x)$). I doubt that you can evaluate the above integral analytically but you can compute it numerically pretty easily.

For a better model, you probably want to couple $G_l$ and $G_u$ via a coupling kernel as per Bonilla, Chai and Williams (2008), Multi-task Gaussian Process Prediction. For multi-task training you need the inputs to be well-defined for all GPs, which is not your case. However, in first approximation (if you are okay with a bit of double dipping) you can fix that by using the uncoupled $G_l$ and $G_u$ to fill in the unseen values. Ideally, you should fill in as uncertain observations of the bounds, with observation-dependent uncertainty as per $\sigma_l$ and $\sigma_u$.

You can actually compute the posterior mean using the following remark if you know the kernel $k$ of the GP. Let $f\sim \mathcal{GP}(0,k)$ from $\mathcal{X}$ to $\mathbb{R}$ (such as $\mathcal{X} \subseteq \mathbb{R}^d$), and $\{x_j\}_{j\in L}$ the set of points for which you have a lower bound $l_j$ and $\{x_j\}_{j\in U}$ the set of points for which you have an upper bound $u_j$. For any point $x\in\mathcal{X}$ you can compute the prediction given by the posterior mean: $$\mathbb{E}\Big[f(x) ~\big|~ \{f(x_j)>l_j\}_{j\in L}, \{f(x_j)<u_j\}_{j\in U}\Big] \\ = \mathbb{E}\Big[\mathbb{E}\big[f(x)\mid \{f(x_j)\}_{j\in L\cup U}\big] \mid \{x_j>l_j\}_{j\in L}, \{x_j<u_j\}_{j\in U}\Big] \\ = \mathbb{E}\Big[\mathrm{k}_x\mathrm{K}_{LU}^{-1}\mathrm{f}_{LU} \mid \{x_j>l_j\}_{j\in L}, \{x_j<u_j\}_{j\in U}\Big] \\ = \mathrm{k}_x\mathrm{K}_{LU}^{-1} \mathbb{E}\Big[\mathrm{f}_{LU} \mid \{x_j>l_j\}_{j\in L}, \{x_j<u_j\}_{j\in U}\Big]\,,$$ where $\mathrm{k}_x$ is the vector of covariances $k(x,x_j)$ for all $x_j\in L\cup U$, and similarly $\mathrm{K}_{LU}$ is the matrix of the $k(x_j,x_{j'})$ and $\mathrm{f}_{LU}$ is a shorthand for the column vector of $f(x_j)$. Then the last expectation is the mean of a truncated multivariate normal. It is hard to find exact formulae, but you have efficient algorithms to sample from this distribution, available in Matlab, R, etc.

Here is an example of the output, where the prediction is in blue, the lower bounds in green and the upper bounds in red. The kernel $k$ is the Squared Exponential with unit length scale.

Without using approximations, the computations costs is driven by the inversion of the covariance matrix $\mathrm{K}_{LU}$, that is $\mathcal{O}(n^3)$ for $n=|L\cup U|$. Cholesky decomposition improves the constant and the stability, but is still $\mathcal{O}(n^3)$. Note that this decomposition is used both in the formula for the posterior mean and for the sampling of the truncated normal. The computational cost does not depend on the dimension $d$, only on the number of training points $n$.

• Nice remark. As far as I understand, this assumes that you already have the GP over the unobserved values. How would you fit the GP hyperparameters? I think that's an important point for what the OP was asking. May 13, 2016 at 18:12
• Hum that's true. I edited the answer to make it clearer. I don't know how the classical methods to estimate the GP hyper-parameters would perform here (maximisation of (pseudo-)likelihood). I can add experiments later if that helps. May 13, 2016 at 18:20
• Problem is, the marginal likelihood is non-analytical due to the non-Gaussian observation likelihood, so you'll have to go with approximate inference (I would try MCMC or variational Bayes here, due to the high uncertainty in the observations). This is the reason why I proposed a solution based on GPs on the bounds as opposed to a GP on the underlying function, since, with the former, one can use the standard analytical machinery. Anyhow, it is interesting to see how different approaches would work. May 13, 2016 at 18:36
• I'll have to think on this answer and the comments here to be honest, but that's due to that I'm rusty and not that the answer in itself isn't clear. To be precise, I'm struggilng with what the $x$ in your $k_x$ is. A variable, such meaning that your final answer is a function of the new $x$ values I would want to evaluate at? Also, I'm not sure how expensive sampling from a truncated multivariate normal is when the number of dimensions is in the order of 50 to 100? At the risk of being a cheapskate, would you be able to somehow share/link to the code used to produce the graph?
– Bram
May 14, 2016 at 15:12
• No problem, I added a short paragraph on computational cost, which does not directly depend on $d$. I could add Matlab code if you think it may clarify the procedure. May 16, 2016 at 10:42