The phase field method has become a widely adopted computational tool for modeling microstructure evolution during martensitic transformations (MTs), owing to its ability to handle complex morphologies without explicit interface tracking [1]. In recent years, machine learning has shown great potential in accelerating phase-field simulations and extracting physically meaningful models from data, as highlighted in reviews on phase-field modeling advancements and data-driven acceleration strategies [2,3]. However, its application to large-scale systems in MTs remains constrained by the necessity of extremely fine mesh resolution at diffuse interfaces, which leads to prohibitive computational costs. While strategies such as manually constructed gradient-based correction functions have been proposed to decouple interfacial energy from numerical thickness [4], these approaches often lack adaptability and precision across varying scales and microstructural configurations. In this work, we introduce a machine learning (ML) framework to autonomously derive an optimized, data-driven gradient correction function for the order parameter field. The model is trained on high-fidelity, small-scale phase field simulations, enabling it to accurately represent the inhomogeneous distribution of order parameter gradients across evolving interfaces. By adaptively adjusting the gradient energy coefficient and bulk free energy density in a thermodynamically consistent manner, the ML-augmented model maintains the correct interfacial energy density while effectively increasing the numerical interface width. This allows for significantly coarser spatial discretization without sacrificing physical accuracy, thereby enabling efficient simulations of MTs in large, representative volume elements.