We consider sample covariance matrices SN=p1ΣN1/2XNXN∗ΣN1/2 where X N is a N × p real or complex matrix with i.i.d. entries with finite 12th moment and ΣN is a N × N positive definite matrix. In addition we assume that the spectral measure of ΣN almost surely converges to some limiting probability distribution as N → ∞ and p/N → γ >0. We quantify the relationship between sample and population eigenvectors by studying the asymptotics of functionals of the type N1Tr(g(ΣN)(SN−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 invers