In this paper, we propose a high-order graph matching formulation to address non-rigid surface matching. The singleton terms capture the geometric and appearance similarities (e.g., curvature and texture) while the high-order terms model the intrinsic embedding energy. The novelty of this paper includes: 1) casting 3D surface registration into a graph matching problem that combines both geometric and appearance similarities and intrinsic embedding information, 2) the first implementation of high-order graph matching algorithm that solves a non-convex optimization problem, and 3) an efficient two-stage optimization approach to constrain the search space for dense surface registration. Our method is validated through a series of experiments demonstrating its accuracy and efficiency, notably in challenging cases of large and/or non-isometric deformations, or meshes that are partially occluded.