Problems involving large-deformation hyperelastic materials in contact remain a central challenge in computational mechanics due to strong nonlinearities and non-smooth constraints. Traditional approaches rely on finite element discretization combined with Newton-type solvers, but convergence difficulties often arise from bifurcations, severe nonlinearity, and the non-differentiable nature of contact conditions. Path-following techniques such as predictor-corrector methods, which are widely used in smooth analysis, are potential candidates to mitigate these issues by introducing an additional constraint and treating the load parameter as an unknown. However, their direct application to non-smooth problems is not well documented, and a rigorous theoretical foundation that incorporates notions of subdifferentiability is therefore required. In this work, we extend the Moore-Penrose continuation method to contact mechanics problems in hyperelasticity, building on prior applications to large-deformation hyperelasticity [1]. We present a detailed convergence analysis of the method for contact problems, leveraging tools from nonsmooth analysis [2,3], which includes subdifferentials, to characterize behavior near non-smooth regions. The proposed framework is compatible with standard contact enforcement strategies, such as penalty and augmented Lagrangian formulations, both widely used in practice for frictionless and frictional contact. Numerical experiments in two and three dimensions demonstrate the effectiveness of the proposed approach in maintaining convergence under challenging conditions. The results indicate that Moore-Penrose continuation is a robust addition for nonlinear contact problems and opens avenues for development of reliable solution strategies and bifurcation detection in contact mechanics.