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