It is well known that the Saint Venant-Kirchhoff model and the incompressible Mooney-Rivlin model result in nonlinear Finite Element (F.E.) models featuring, exclusively, quadratic and cubic nonlinearities with respect to the displacement vector[1, 2]. In this contribution, these nonlinearities are represented, respectively, by third and forth order symmetric sparse tensors which are never explicitly computed nor stored in memory in order to ease the computational and spatial (memory) complexity of the models. By using the harmonic balance method (HBM)[3], semi-analytical equations can be derived for the free and forced vibration of geometrically nonlinear structures. The HBM however, multiplies the effective size of the problems to be solved by the number of harmonics considered, which quickly becomes prohibitive for large scale F.E. models. We therefore propose combining the semi-analytical HBM equations with an implicit moment-matching[4] technique to tackle this issue. The resulting methodology is applied to speed up the modal analysis and the frequency response computation of geometrically nonlinear structures.