Equation error, also known as one-step-ahead prediction error, is a common quality-of-fit metric in dynamical system identification and learning. In this letter, we use Lagrangian relaxation to construct a convex upper bound on equation error that can be optimized over a convex set of nonlinear models that are guaranteed to be contracting, a strong form of nonlinear stability. We provide theoretical results on the tightness of the relaxation, and show that the method compares favorably to established methods on a variety of case studies.