Gaussian process smoothing with eigen-decomposition From the book 'Gaussian Process for Machine learning', the author talked about the smoothing w.r.t eigen-decomposition. The prediction $f$ at training points can be written as
$$f=K(K+\sigma_n^2I)^{-1}y$$
where $K$ is the covariance matrix and $\sigma_n$ is the noise variance. Note that this is the predictions for training points.
Then the author use eigen-decomposition on $K$ such that $K=\sum^n_{i=1}\lambda_iu_iu_i^T$ where $n$ is the number of training points. Then, let the responses linear combined by $u_i$ such that $y=\sum^n_{i=1}\gamma_iu_i$. Next, the author gives the formula
$$f=\sum^n_{i=1}\frac{\gamma_i\lambda_i}{\lambda_i+\sigma_n^2}$$
I don't get how can the outer product $u_iu_i^T$ cancel out. I tested the result and it is indeed correct. Any one knows the detail?
 A: The version that I found has a different formula for $f$. Equation (2.36) in Section 2.6 (Smoothing, Weight Functions and Equivalent Kernels) states the following:
\begin{equation}
f = \sum_{i=1}^{n}\frac{\gamma_i\lambda_i}{\lambda_i + \sigma_n^2}u_i.
\end{equation}
This makes sense, as $K(K + \sigma_n^2I)^{-1}$ is a matrix and $y$ is a vector, so $f$ should be a vector. I will show how to get this formula.

Since $K$ is a real covariance matrix, it is a positive-definite real matrix, which guarantees that it is orthogonally diagonalizable. This means that the $u_i$ in the expression
\begin{equation}
K = \sum_{i=1}^{n}\lambda_i u_i u_i^{\mathsf{T}}
\end{equation}
form an orthonormal basis of $\mathbb{R}^n$.
\begin{eqnarray}
f &=& K(K + \sigma_n^2I)^{-1}y\\
&=& \left(\sum_{i=1}^{n}\lambda_i u_i u_i^{\mathsf{T}}\right)\left(\sum_{j=1}^{n}\lambda_j u_j u_j^{\mathsf{T}} + \sigma_n^2 I\right)^{-1}\left(\sum_{\ell=1}^{n}\gamma_{\ell} u_{\ell}\right)\\
&=& \left(\sum_{i=1}^{n}\lambda_i u_i u_i^{\mathsf{T}}\right)\sigma_n^{-2}\left( \sigma_n^{-2}\sum_{j=1}^{n}\lambda_j u_j u_j^{\mathsf{T}} +I\right)^{-1}\left(\sum_{\ell=1}^{n}\gamma_{\ell} u_{\ell}\right)
\end{eqnarray}
Recall the Taylor-MacLaurin series for $1 - x$:
\begin{equation}
1 - x = \sum_{n=0}^{\infty}x^n~~\textrm{for}~~|x| < 1.
\end{equation}
Then
\begin{equation}
1+ x = \sum_{n=0}^{\infty}(-x)^n~~\textrm{for}~~|x| < 1.
\end{equation}
An analogous formula can be used for the square-matrix analogue:
\begin{equation}
I + A = \sum_{m=0}^{\infty}(-A)^m~~\textrm{for}~~\left\Arrowvert A\right\Arrowvert < 1
\end{equation}
or if $A^N = O$, the zero matrix, for some $N$. This latter condition cannot hold for a positive-definite matrix.
Let's assume that this will work for $(\sigma_n^{-2}K + I)^{-1}$:
\begin{equation}
\left(I + \sigma_n^{-2}K\right)^{-1} = \sum_{m=0}^{\infty}(-\sigma_n^{-2})^mK^m.
\end{equation}
A valuable thing to notice about orthogonally diagonalizable matrices such as $K$ is that
\begin{equation}
K^m = \left(\sum_{i=1}^{n}\lambda_i u_iu_i^{\mathsf{T}}\right)^m = \sum_{i=1}^{n}\lambda_i^m u_iu_i^{\mathsf{T}}.
\end{equation}
I'll let you verify that. Use the fact that the $u_i$ form an orthonormal basis.
Now we can write $f$ as
\begin{eqnarray}
f &=& \sigma_n^{-2} \color{red}{K}\left[\sum_{m=0}^{\infty}(-\sigma_n^{-2})^m\color{red}{K^m}\right]y\\
&=& \sigma_n^{-2}\left[\sum_{m=0}^{\infty}(-\sigma_n^{-2})^m\color{red}{K^{m+1}}\right]y\\
&=& \sigma_n^{-2}\left[\sum_{m=0}^{\infty}(-\sigma_n^{-2})^m\color{red}{\sum_{i=1}^{n}\lambda_i^{m+1}u_i u_i^{\mathsf{T}}}\right]\sum_{\ell=1}^{n}\gamma_{\ell} u_{\ell}\\
&=& \sigma_n^{-2}\sum_{m=0}^{\infty}(-\sigma_n^{-2})^m\left[\color{red}{\sum_{i=1}^{n}\lambda_i^{m+1}u_i}\sum_{\ell=1}^{n}\gamma_{\ell}\color{red}{u_i^{\mathsf{T}}}u_{\ell}\right]
\end{eqnarray}
Since the $u_i$ are orthonormal,
\begin{equation}
u_i^{\mathsf{T}}u_{\ell} = \delta_{i,\ell},
\end{equation}
the Kronecker delta function. When we sum over $\ell$, the only term that "survives" will be the term for which $\ell = i$.
\begin{eqnarray}
f &=& \sigma_n^{-2}\sum_{m=0}^{\infty}(-\sigma_n^{-2})^m\left[\sum_{i=1}^{n}\lambda_i^{m+1}u_i\sum_{\ell=1}^{n}\gamma_{\ell}\delta_{i,\ell}\right]\\
&=& \sigma_n^{-2}\sum_{m=0}^{\infty}(-\sigma_n^{-2})^m\left[\sum_{i=1}^{n}\gamma_i\lambda_i^{m+1}u_i\right]\\
&=& \sigma_n^{-2}\sum_{i=1}^{n}\left[\sum_{m=0}^{\infty}(-\sigma_n^{-2})^m\lambda_i^{m}\right]\lambda_i\gamma_i u_i\\
&=& \sigma_n^{-2}\sum_{i=1}^{n}\left[\sum_{m=0}^{\infty}(-\sigma_n^{-2}\lambda_i)^m\right]\lambda_i\gamma_i u_i\\
&=& \sigma_n^{-2}\sum_{i=1}^{n}\left[(1 + \sigma_n^{-2}\lambda_i)^{-1}\right]\lambda_i\gamma_i u_i\\
&=& \sigma_n^{-2}\sum_{i=1}^{n}\frac{\gamma_i\lambda_i}{1 + \sigma_n^{-2}\lambda_i}u_i\\
&=& \sum_{i=1}^{n}\frac{\gamma_i\lambda_i}{\sigma_n^{2} + \lambda_i}u_i
\end{eqnarray}

Upon questioning from @yiming-che, I realized an easier approach. This approach is easy to see because $K$ is orthogonally diagonalizable. First, note that since the $u_i$ form an orthogonormal basis for $\mathbb{R}^n$, $I = \sum_{i=1}^{n}u_iu_i^{\mathsf{T}}$. This gives us a simpler expression for $I + \sigma_n^{-2}K$.
\begin{eqnarray}
\left(I + \sigma_n^{-2}K\right)^{-1} &=& \left(\sum_{i=1}^{n}u_iu_i^{\mathsf{T}} + \sum_{i=1}^{n}\sigma_n^{-2}\lambda_i u_iu_i^{\mathsf{T}}\right)^{-1}\\
&=& \left(\sum_{i=1}^{n}(1 + \sigma_n^{-2}\lambda_i)u_iu_i^{\mathsf{T}}\right)^{-1}
\end{eqnarray}
I claim that the last line can be re-written as
\begin{equation}
\left(\sum_{i=1}^{n}\left(1 + \sigma_n^{-2}\lambda_i\right)u_iu_i^{\mathsf{T}}\right)^{-1} = \sum_{i=1}^{n}\left(1 + \sigma_n^{-2}\lambda_i\right)^{-1}u_iu_i^{\mathsf{T}}.
\end{equation}
Let's prove it by multiplying the original matrix and the suspected inverse.
\begin{eqnarray}
&&\left[\sum_{j=1}^{n}\left(1 + \sigma_n^{-2}\lambda_j\right)u_ju_j^{\mathsf{T}}\right]\left[\sum_{i=1}^{n}\left(1 + \sigma_n^{-2}\lambda_i\right)^{-1}u_iu_i^{\mathsf{T}}\right]\\
&=& \sum_{i=1}^{n}\sum_{j=1}^{n}\left(1 + \sigma_n^{-2}\lambda_j\right)\left(1 + \sigma_n^{-2}\lambda_i\right)^{-1}u_j~\underbrace{u_j^{\mathsf{T}}u_i}_{\delta_{i,j}}~u_i^{\mathsf{T}}\\
&=& \sum_{i=1}^{n}\sum_{j=1}^{n}\left(1 + \sigma_n^{-2}\lambda_j\right)\left(1 + \sigma_n^{-2}\lambda_i\right)^{-1}u_j\delta_{i,j}u_i^{\mathsf{T}}\\
&=& \sum_{i=1}^{n}\left(1 + \sigma_n^{-2}\lambda_i\right)\left(1 + \sigma_n^{-2}\lambda_i\right)^{-1}u_iu_i^{\mathsf{T}}\\
&=& \sum_{i=1}^{n}u_iu_i^{\mathsf{T}}\\
&=& 
\end{eqnarray}
This works as long as $1 + \sigma_n^{-2}\lambda_i \neq 0$ for each $i$.
