Abstract: We develop an optimization-based maximum likelihood approach to analyze the cross-correlation of the Cosmic Microwave Background (CMB) and the Large-Scale Structure (LSS) induced by the kinetic Sunyaev-Zeldovich effect. Our main goal is to reconstruct the radial velocity field of the universe. While the existing quadratic estimator (QE) is statistically optimal for current and near-term experiments, the likelihood can extract more signal-to-noise in the future. Our likelihood formulation has further advantages over the QE, such as the possibility of jointly fitting cosmological and astrophysical parameters. We implement an auto-differentiable likelihood pipeline in JAX, which is computationally tractable for a realistic survey resolution, and evaluate it on the Agora simulation.