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).
 A: 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.  
A: 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)).$$
Or, using your notation,
$$\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)$.
