Modeling a high-dimensional Hamiltonian system in reduced dimensions with respect to coarse-grained (CG) variables can greatly reduce computational cost and enable efficient bottom-up prediction of main features of the system for many applications. However, it usually experiences significantly altered dynamics due to loss of degrees of freedom upon coarse-graining. To establish CG models that can faithfully preserve dynamics, previous efforts mainly focused on equilibrium systems. In contrast, various soft matter systems are known out of equilibrium. Therefore, the present work concerns non-equilibrium systems and enables accurate and efficient CG modeling that preserves non-equilibrium dynamics and is generally applicable to any non-equilibrium process and any observable of interest. To this end, the dynamic equation of a CG variable is built in the form of the non-stationary generalized Langevin equation (nsGLE) to account for the dependence of non-equilibrium processes on the initial conditions, where the two-time memory kernel is determined from the data of the two-time auto-correlation function of the non-equilibrium trajectory-averaged observable of interest. By embedding the non-stationary non-Markovian process in an extended stochastic framework, an explicit form of the non-stationary random noise in the nsGLE is introduced, and the cost is significantly reduced for solving the nsGLE to predict the non-equilibrium dynamics of the CG variable. To prove and exploit the equivalence of the nsGLE and extended dynamics, the memory kernel is parameterized in a two-time exponential expansion. A data-driven hybrid optimization process is proposed for the parameterization, a non-convex and high-dimensional optimization problem.