In this paper, we are concerned with the parameterized least squares inverse eigenvalue problems for the case that the number of parameters to be constructed is less than the number of prescribed realizable eigenvalues. Through equivalent transformation, the original problem becomes a nonlinear least squares problem associated with a specific over-determined mapping defined between a Riemannian manifold and a Euclidean space. We propose the Riemannian two-step perturbed Gauss–Newton method combined with a specific second-order nonmonotone backtracking line search technique for solving general nonlinear least squares problem on Riemannian manifold. Global convergence of this algorithm is discussed under some mild assumptions. Meanwhile, a cubical convergence rate is obtained under injectivity of the differential of the underlying map and zero residue of this map at an accumulation point. To apply the proposed method to solving the parameterized least squares inverse eigenvalue problems, exact solution of the perturbed Riemannian Gauss–Newton equation is constructed. Finally, numerical experiments show the efficiency of the proposed method.