# Quantifying the effect of noise on the posterior variance in Gaussian processes / multivariate Gaussian vectors

Consider a real-valued Gaussian process $f$ on some compact domain $\mathcal{X}$ with mean zero and covariance function $k(x,x') \in [0,1]$ (also known as the kernel function). This question concerns a finite collection of points in $\mathcal{X}$, namely, a set of sampled points $\mathbf{x} = [x_1,\dotsc,x_n]^T$ and a query point $x$. Writing $\mathbf{f}(\mathbf{x}) = [f(x_1),\dotsc,f(x_n)]^T$, we have the joint distribution $$\left[ \begin{array}{c} \mathbf{f}(\mathbf{x}) \\ f(x) \end{array}\right] \sim N\left( \left[ \begin{array}{c} \mathbf{0} \\ 0 \end{array}\right], \left[ \begin{array}{cc} \mathbf{K} & \mathbf{k}(x) \\ \mathbf{k}(x)^T & k(x,x) \end{array}\right] \right),$$ where $\mathbf{k}(x)$ is an $n\times 1$ vector with $i$-th entry $k(x,x_i)$, and $\mathbf{K}$ is an $n \times n$ matrix with $(i,j)$-th entry $k(x_i,x_j)$. The samples corresponding to $\mathbf{x} = [x_1,\dotsc,x_n]^T$ are denoted by $\mathbf{y} = [y_1,\dotsc,y_n]^T$, and take the form $$y_i = f(x_i) + z_i,$$ where $z_i \sim N(0,\sigma^2)$ is additive Gaussian noise (independent for each sample) with $\sigma^2 \le 1$.

It is well known that the posterior distribution of $f(x)$ given $\mathbf{y}$ (with $\mathbf{x}$ assumed fixed and known) is Gaussian, with the posterior mean and variance taking the form $$\mu_n(x) = \mathbf{k}(x)^T(\mathbf{K} + \sigma^2 \mathbf{I})^{-1}\mathbf{y}$$ $$\sigma_n^2(x) = k(x,x) - \mathbf{k}(x)^T(\mathbf{K} + \sigma^2 \mathbf{I})^{-1}\mathbf{k}(x).$$ My question is as follows: If we let $\widetilde{\sigma}_n^2(x) = k(x,x) - \mathbf{k}(x)^T \mathbf{K}^{-1}\mathbf{k}(x)$ be the posterior variance we would get under noiseless samples, then is it true that $$\sigma_n^2(x) \le \widetilde{\sigma}_n^2(x) + C\sigma^2$$ for some universal constant $C$? Intuitively, if our samples are each corrupted by noise of variance $\sigma^2$, we shouldn't expect to incur more than $O(\sigma^2)$ additional uncertainty on the unknown function value $f(x)$ that we are trying to predict.

Notes: A potential starting point is to use the Woodbury matrix identity to write $$\mathbf{k}(x)^T(\mathbf{K} + \sigma^2 \mathbf{I})^{-1}\mathbf{k}(x) = \mathbf{k}(x)^T \mathbf{K}^{-1}\mathbf{k}(x) + \sigma^2\mathbf{k}(x)^T \Big(\mathbf{K}^{-1} \big(\mathbf{I} + \sigma^2\mathbf{K}^{-1} \big)^{-1} \mathbf{K}^{-1}\Big)\mathbf{k}(x).$$ By a matrix Taylor expansion, the final term should behave as $O(\sigma^2)$ as $\sigma^2 \to 0$, which appears to yield the desired result. However, this approach leads to a constant factor depending on $\mathbf{x}$ and $n$, whereas I would like to show the above result with an absolute constant $C$.

Having said that, if it makes things easier, I would be happy for $C$ to depend on the covariance function $k$, and/or on the input domain $\mathcal{X}$ (e.g., even the simple choices $\mathcal{X} = [0,1]$ and $k(x,x') = e^{-c\cdot(x-x')^2}$ would be of interest).

No, you would generally expect it to depend on n. Consider $k(x,y) = 1$, which you may want to rule out. Here you are observing a single value with error, and you average and get a variance like $\frac 1 n$. If you don't like that, consider a quite smooth Gaussian where you observe m sets of n quite close together points, say near $x_1,..., x_m$. In the error free case the extra observations do you no good, and the error is essentially that of observing the process at $x_1,...,x_m$. In the case with error, average each of the m sets, and you again have the value of the process at $x_1, ..., x_m$, with much smaller error. Cleary (hope you don't mind that) as $n \rightarrow \infty$ you will do as well with the observations with error.

• I believe that in this example, the property I am seeking does hold. If k(x,y)=1 then we get sigma_noiseless^2 = 0, and sigma_noisy^2 = sigma^2 / n, so it is certainly true that sigma_noisy^2 <= sigma_noiseless^2 + C*sigma^2. – jmscarlett Jan 29 '18 at 23:59
• Sorry, you weren't asking the question I thought you were asking. This is ridge regression, try en.wikipedia.org/wiki/Tikhonov_regularization. I think the decomposition of RSS = RSS_0 + ... gives a constant C that is like $||Y||^2$ – user83457 Jan 30 '18 at 11:32

Consider the joint Gaussian distribution of $(Y, Z, f(x))$. Observe that knowing both $Y$ and $Z$ together is equivalent to knowing $f(\mathbf{x})$ (the noiseless version of $Y$). Then we can compute $\mbox{Var}(f(x) \mid f(\mathbf{x}))$ as $\mbox{Var}(f(x) \mid Y, Z)$. Furthermore, we can compute $\mbox{Var}(f(x) \mid Y)$ via the law of total variance:

$$\mbox{Var}(f(x) \mid Y) = \mbox{E}(\mbox{Var}(f(x) \mid Y, Z)) + \mbox{Var}(\mbox{E}(f(x) \mid Y, Z)).$$

$$\sigma^2_n(x) = \tilde{\sigma}^2_n(x) + \mbox{Var}(\mathbf{k}^t(x)\mathbf{K}^{-1}(Y - Z))$$
where the second variance is with respect to the random variable $Z$. Simplifying further yields
$$\sigma^2_n(x) = \tilde{\sigma}^2_n(x) + \sigma^2 \mathbf{k}^t(x)\mathbf{K}^{-2}\mathbf{k}(x).$$
Accordingly, it seems that $C = \mathbf{k}^t(x)\mathbf{K}^{-2}\mathbf{k}(x)$.