Understanding how complex systems transition between states requires mapping the energy landscape that governs these changes. Local transition-state networks reveal the barrier architecture that explains observed behaviour and enables mechanism-based prediction across computational chemistry, biology, and physics, yet current practice either prescribes endpoints or randomly samples only a few saddles around an initial guess. We present a general optimization framework that systematically expands local coverage by coupling a multi-objective explorer with a bilayer minimum-mode kernel. The inner layer uses Hessian-vector products to recover the lowest-curvature subspace (smallest k eigenpairs), the outer layer optimizes on a reflected force to reach index-1 saddles, then a two-sided descent certifies connectivity. The GPU-based pipeline is portable across autodiff backends and eigensolvers and, on large atomistic-spin tests, matches explicit-Hessian accuracy while cutting peak memory and wall time by orders of magnitude. Applied to a DFT-parameterized Néel-type skyrmionic model, it recovers known routes and reveals previously unreported mechanisms, including meron-antimeron-mediated Néel-type skyrmionic duplication, annihilation, and chiral-droplet formation, enabling up to 32 pathways between biskyrmion (Q=2) and biantiskyrmion (Q=-2). The same core transfers to Cartesian atoms, automatically mapping canonical rearrangements of a Ni(111) heptamer, underscoring the framework's generality.