In this letter, a novel optimization method for the data-driven local coordinates approach to a prediction error method, which is one of system identification methods, is developed. In the proposed method, parameters of system matrices to be estimated are kept as their original forms of matrices whereas they are converted into a vector in the literature. Optimization with respect to the matrix variables can improve the computational efficiency. Furthermore, it is effective to equate input–output equivalent systems with each other since they attain the same value of the objective function. To this end, a Riemannian quotient manifold is introduced. The geometry of the proposed manifold is investigated and a novel Riemannian conjugate gradient method on the manifold is provided. Numerical experiments show that the difference of the proposed Riemannian conjugate gradient method from Euclidean one greatly improves the efficiency of the algorithm.