We consider the numerical approximation of the inertial Landau-Lifshitz-Gilbert (iLLG) equation, which describes the dynamics of the magnetization in ferromagnetic materials at subpicosecond time scales. We propose and analyze two fully discrete numerical schemes based on two different approaches: The first method is based on a reformulation of the problem as a linear constrained variational formulation for the time derivative of the magnetization. The second method exploits a reformulation of the problem as a first order system in time for the magnetization and the angular momentum. Both schemes are implicit, based on first-order finite elements, and the constructed numerical approximations satisfy the inherent unit-length constraint of iLLG at the vertices of the underlying mesh. For both schemes, we establish a discrete energy law and prove convergence of the approximations towards a weak solution of the problem. Numerical experiments validate the theoretical results and show the applicability of the methods for the simulation of ultrafast magnetic processes.