Gaussian process - Why adding data points cannot increase the predictive bias? I've seen this question here:
How to increase variance in Gaussian Process regression?
And trying to complete the proof. 
I'm looking at this book: 
Rasmussen & Williams 2006: Gaussian Processes for Machine Learning
Page 31 starts the relevant subject, and question 2.9.4 is just the question I want to ask:

Let $\operatorname{var}_{n}\left(f\left(\mathbf{x}_{*}\right)\right)$ be the predictive variance of a Gaussian process regression model at $x∗$ given a dataset of size $n$. The corresponding predictive
  variance using a dataset of only the first $n − 1$ training points is denoted $\operatorname{var}_{n-1}\left(f\left(\mathbf{x}_{*}\right)\right)$. Show that $\operatorname{var}_{n}\left(f\left(\mathbf{x}_{*}\right)\right) \leq \operatorname{var}_{n-1}\left(f\left(\mathbf{x}_{*}\right)\right)$, i.e. that
  the predictive variance at x∗ cannot increase as more training data is obtained. One way to approach this problem is to use the partitioned matrix
  equations given in section A.3 to decompose 
  $\operatorname{var}_{n}\left(f\left(\mathbf{x}_{*}\right)\right)=k\left(\mathbf{x}_{*}, \mathbf{x}_{*}\right)-\mathbf{k}_{*}^{\top}\left(K+\sigma_{n}^{2} I\right)^{-1} \mathbf{k}_{*}$.

I want to represent $K_{n}+\sigma^{2} I_{n}$ as a function of $K_{n-1}$, and then to decompose $K_{n}+\sigma^{2} I_{n}$ by the following lemma from the book : (A.11 page 219)
Let the invertible $n \times n$ matrix $A$ and its inverse $A^{-1}$ be partitioned into
$$ {A=\left( \begin{array}{cc}{P} & {Q} \\ {R} & {S}\end{array}\right), \quad A^{-1}=\left( \begin{array}{cc}{\tilde{P}} & {\tilde{Q}} \\ {\tilde{R}} & {\tilde{S}}\end{array}\right)} $$
where $P$ and $\tilde{P}$ are $n1 \times n1$ matrices and S and $\tilde{S}$ are $n2 \times n2$ matrices with
$n = n1 + n2$.
The submatrices of $A^{-1}$ are given in Press et al. [1992, p. 77] as:
$$\tilde{P}=P^{-1}+P^{-1} Q M R P^{-1}$$
$$\tilde{Q}=-P^{-1} Q M$$
$$\tilde{R}=-M R P^{-1}$$
$$\tilde{S}=M$$, Where $M=\left(S-R P^{-1} Q\right)^{-1}$ 
But then I don't know how to reach $k_{*}^{\top}\left(K_{n}+\sigma^{2} I_{n}\right)^{-1} k_{*} \geq k_{n-1}\left(x^{*}\right)^{\top}\left(K_{n-1}+\sigma^{2} I_{n-1}\right)^{-1} k_{n-1}\left(x^{*}\right)$ in order to complete the proof. 
It will finish the proof because of $\operatorname{var}_{n}\left(f\left(\mathbf{x}_{*}\right)\right) \leq \operatorname{var}_{n-1}\left(f\left(\mathbf{x}_{*}\right)\right)$ and the fact that $\operatorname{var}_{n}\left(f\left(\mathbf{x}_{*}\right)\right)=k\left(\mathbf{x}_{*}, \mathbf{x}_{*}\right)-\mathbf{k}_{*}^{\top}\left(K_{n}+\sigma^{2} I_{n}\right)^{-1} \mathbf{k}_{*}$...
How can follow the clue there? Thanks
 A: Recall the expression for the posterior variance of a Gaussian Process:  $\operatorname{var}_{n}\left(f\left(\mathbf{x}_{*}\right)\right)=k\left(\mathbf{x}_{*}, \mathbf{x}_{*}\right)-\mathbf{k}_{*}^{\top}\left(K_n+\sigma^{2} I_n\right)^{-1} \mathbf{k}_{*}$. We want to show that $\text{var}_{n-1}(f(x_*)) \geq \text{var}_n(f(x_*))$. The first term is the same for both $\text{var}_n$ and $\text{var}_{n-1}$, so we just need to break down the second term in a way such that we can compare between $\text{var}_n$ and $\text{var}_{n-1}$. 
Following the hint to look at section A.3 in the book, we see that for the partitioned matrix $A$ and its inverse $A^{-1}$, 
\begin{align}
A = \begin{pmatrix} P & Q \\ R & S \end{pmatrix}, \qquad A^{-1} = \begin{pmatrix} \tilde P & \tilde Q \\ \tilde R & \tilde S \end{pmatrix},
\end{align}
\begin{align}
\tilde P & = P^{-1} + P^{-1}QMRP^{-1},\\
\tilde Q & = -P^{-1}QM,\\
\tilde R & = -MRP^{-1},\\
\tilde S & = M,\\
\text{where}\  M & =  (S - RP^{-1}Q)^{-1}.\\
\end{align}
We can decompose the Gram matrix with added noise $K_n + \sigma^2I_n$ as follows:
\begin{align}
K_n + \sigma^2 I_n = 
\begin{pmatrix}
K_{n-1} + \sigma^2 I_{n-1} & k_{n-1}(x')\\k_{n-1}(x')^\top & k(x',x') + \sigma^2
\end{pmatrix}
\end{align}
where $x'$ is the $n$th point sampled. 
Writing out the decomposed inverse, we get 
\begin{align}
(K_n + \sigma^2 I_n)^{-1} = \begin{pmatrix}\kappa + \kappa k_{n-1}(x')Mk_{n-1}(x')^\top\kappa & -\kappa k_{n-1}(x')M \\ -Mk_{n-1}(x')^\top \kappa & M \end{pmatrix},
\end{align}
with 
\begin{align}
M & = (k(x', x') + \sigma^2 - k_{n-1}(x')^\top(K_{n-1} + \sigma^2I_{n-1})^{-1}k_{n-1}(x'))^{-1}, \\
\kappa & = (K_{n-1} + \sigma^2I_{n-1})^{-1}
\end{align}
Finally we can actually compute the term $k_*^\top(K_n + \sigma^2I_n)^{-1}k_*$ by multiplying through:
\begin{align}
k_*^\top(K_n + \sigma^2I_n)^{-1}k_* = & \begin{pmatrix}k_{n-1}(x^*)\\k'(x^*)\end{pmatrix}^\top\begin{pmatrix}\kappa + \kappa k_{n-1}(x')Mk_{n-1}(x')^\top\kappa & -\kappa k_{n-1}(x')M \\ -Mk_{n-1}(x')^\top\kappa  & M \end{pmatrix}\begin{pmatrix}k_{n-1}(x^*)\\k'(x^*)\end{pmatrix}\\
& = \begin{pmatrix} k_{n-1}^\top(x^*)\kappa + k_{n-1}^\top(x^*)\kappa k_{n-1}(x')Mk_{n-1}(x')^\top\kappa - k'(x^*)Mk_{n-1}(x')^\top \kappa \\ -k_{n-1}^\top(x^*)\kappa k_{n-1}(x')M + k'(x^*)M \end{pmatrix} ^\top
\begin{pmatrix}k_{n-1}(x^*)\\k'(x^*)\end{pmatrix}
\end{align}
where $\kappa = (K_{n-1} + \sigma^2I_{n-1})^{-1}$ and $k'(x^*) = k(x', x^*)$.
We end up with 
\begin{align}
k_*^\top(K_n + \sigma^2I_n)^{-1}k_* & = k_{n-1}^\top(x^*)\kappa k_{n-1}(x^*) + k_{n-1}^\top(x^*)\kappa k_{n-1}(x')Mk_{n-1}(x')^\top\kappa k_{n-1}(x^*)\\
& - k'(x^*)Mk_{n-1}(x')^\top\kappa k_{n-1}(x^*) - k_{n-1}(x^*)^\top\kappa k_{n-1}(x')Mk'(x^*) + k'(x^*)Mk'(x^*).
\end{align}
Note that the first term is simply the relevant term from $\text{var}_{n-1}$, so we just need to show that the subsequent terms are non-negative in order to show that the term is larger for $n$ compared to $n-1$ (and so the variance is smaller). Note that $M$ is the reciprocal of the variance at $x'$ plus $\sigma^2$, so it's positive and we can factor it out. 
We are left with 
\begin{align}
\alpha^2 - 2k'(x^*)\alpha + k'(x^*)^2,
\end{align}
where $\alpha = k_{n-1}^\top(x^*)\kappa k_{n-1}(x')$ (a scalar). So we have 
\begin{align}
k_*^\top(K_n + \sigma^2I_n)^{-1}k_* & = k_{n-1}^\top(x^*)\kappa k_{n-1}(x^*) + (\alpha - k'(x^*))^2\\
& = k_{n-1}(x^*)^\top(K_{n-1} + \sigma^2I_{n-1})^{-1}k_{n-1}(x^*) + \tfrac{1}{M}(\alpha - k'(x^*))^2\\
& \geq k_{n-1}(x^*)^\top(K_{n-1} + \sigma^2I_{n-1})^{-1}k_{n-1}(x^*).
\end{align}
So the variance after $n$ points is smaller than the variance after $n-1$ points, with equality achieved when the quantity $\tfrac{1}{M}(k_{n-1}^\top(x^*)(K_{n-1} + \sigma^2I_{n-1})^{-1} k_{n-1}(x') - k(x^*, x'))$ is zero. 
