Recently developed permanent magnets, featuring specially engineered microstructures of inhomogeneous magnetic phases, are being considered as cost-effective alternatives to homogeneous single-main-phase hard magnets composed of Nd2Fe14B, without compromising performance. In this study, we conducted a comprehensive examination of a core–shell sphere cluster model of Ce-substituted inhomogeneous Nd2-δCeδFe14B phases versus homogeneous magnetic phases, utilizing finite-element micromagnetic simulation and machine learning methods. This involved a meticulous, sphere-by-sphere analysis of individual demagnetization curves calculated from the cluster model. The grain-by-grain analyses unveiled that these individual demagnetization curves can elucidate the overall magnetization reversal in terms of the nucleation and coercive fields for each sphere. Furthermore, it was observed that Nd-rich spheres exhibited much broader ranges of nucleation and coercive field distributions, while Nd-lean spheres showed relatively narrower ranges. To identify the key parameter responsible for the notable differences in the nucleation fields, we constructed a machine learning regression model. The model utilized numerous hyperparameter sets, optimized through the very fast simulated annealing algorithm, to ensure reliable training. Using the kernel SHapley Additive eXplanation (SHAP) technique, we inferred that stray fields among the 11 parameters were closely related to coercivity. We further substantiated the machine learning models’ inference by establishing an analytical model based on the eigenvalue problem in classical micromagnetic theory. Our grain-by-grain interpretation can guide the optimal design of granular hard magnets from Nd2Fe14B and other abundant rare earth transition elements, focusing on extraordinary performance through the careful adjustment of microstructures and elemental compositions.