Hamiltonian systems provide a canonical description of purely conservative dynamics and serve as idealized reference models for real physical systems. In engineering applications involving large-scale Hamiltonian systems, it is often necessary to reduce the dimensionality of the problem to enable efficient analysis and design. Data-driven reduction methods are particularly appealing in this context, as they allow for the construction of reduced-order models (ROMs) without explicit knowledge of the underlying governing equations. However, standard approaches may fail to preserve fundamental structural properties of the original system: in the case of Hamiltonian dynamics, this typically results in ROMs that do not conserve energy and lack symplectic structure. To address these limitations, the work of [Gruber and Tezaur, SIAM J. Appl. Dyn. Syst., 2025] introduces a model reduction framework that yields variationally-consistent Hamiltonian ROMs by construction. While the advances in [Gruber and Tezaur, SIAM J. Appl. Dyn. Syst., 2025] of choices that satisfy the desired structural constraints. It is therefore natural to ask whether it is possible to find the optimal/most accurate among these. In this work, we address this question using the recently-introduced NiTROM formulation [Padovan et al., SIAM J. Appl. Dyn. Syst., 2024], which is designed to identify optimal ROMs from high-fidelity trajectory data by solving a non-intrusive optimization problem over a carefully chosen matrix manifold. In particular, we show that the desired variationally-consistent structure can be strongly enforced with minor modifications to the NiTROM framework, and the optimization problem can be solved straightforwardly using existing descent algorithms on matrix manifolds with closed-form gradient computations. We test our formulation on engineering applications of increasing complexity, demonstrating both energy preservation and superior predictive accuracy of the reduced models.