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: 


*

*$(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$.

*$(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$
 A: 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$.
A: 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$.
