In many applications, it is desired to obtain extreme eigenvalues and eigenvectors of large Hermitian matrices by efficient and compact algorithms. In particular, orthogonalization-free methods are preferred for large-scale problems for finding eigenspaces of extreme eigenvalues without explicitly computing orthogonal vectors in each iteration. For the top p eigenvalues, the simplest orthogonalization-free method is to find the best rank-p approximation to a positive semi-definite Hermitian matrix by algorithms solving the unconstrained Burer–Monteiro formulation. We show that the nonlinear conjugate gradient method for the unconstrained Burer–Monteiro formulation is equivalent to a Riemannian conjugate gradient method on a quotient manifold with the Bures–Wasserstein metric, thus its global convergence to a stationary point can be proven. Numerical tests suggest that it is efficient for computing the largest k eigenvalues for large-scale matrices if the largest k eigenvalues are nearly distributed uniformly.