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).
