In this paper, we propose and analyze a Riemannian Gauss–Newton method for the low-rank matrix completion problem. Different from the existing second-order Riemannian optimization methods for this problem, we adopt the approximate Riemannian Hessian of Gauss–Newton methods instead of the true Riemannian Hessian. This approximate Riemannian Hessian has a computationally efficient and symmetric positive semidefinite structure. Added with a regularization term, the corresponding Riemannian Gauss–Newton equation can be solved efficiently by the linear conjugate gradient method. Global and local convergence of the proposed method are established under mild assumptions. Preliminary numerical results on synthetic and image recovery problems with various dimensionalities and ranks are reported to demonstrate the efficiency of the proposed method.