We consider sample covariance matrices $S_N=\frac{1}{p}Σ_N^{1/2}X_NX_N^* Σ_N^{1/2}$ where $X_N$ is a $N \times p$ real or complex matrix with i.i.d. entries with finite $12^{\rm th}$ moment and $Σ_N$ is a $N \times N$ positive definite matrix. In addition we assume that the spectral measure of $Σ_N$ almost surely converges to some limiting probability distribution as $N \to \infty$ and $p/N \to γ>0.$ We quantify the relationship between sample and population eigenvectors by studying the asymptotics of functionals of the type $\frac{1}{N} \text{Tr} (g(Σ_N) (S_N-zI)^{-1})),$ where $I$ is the identity matrix, $g$ is a bounded function and $z$ is a complex number. This is then used to compute the asymptotically optimal bias correction for sample eigenvalues, paving the way for a new generation of improved estimators of the covariance matrix and its inverse.