Lagrange multiplier-based (LM) coupling approaches have been shown to be stable, accurate, and applicable to multiphysics problems coupled over an interface. For explicitly time-integrated coupled problems (IVR), the resulting Schur complement that is used to solve for the coupling traction forces is tractable but requires access to operators generally embedded deeply within the software modeling the subdomain. For implicitly time-integrated problems (IFR), the Schur complement involves inversion of stiffness matrices and is not computationally feasible. Through a data-driven discovery of coupling operators, we are able to loosen the requirement for access to multiple complex data structures and the resulting operator is more compact and less computationally expensive to evaluate than the full-order coupling operator. With the rise in popularity of data-driven modeling techniques applied to subdomain problems, we make the assumption of the availability of full-resolution subdomain solution data on at least one subdomain and then leverage an idea by Carey et. al. from which to generate high-order approximations of traction forces for snapshots. We use the traction force snapshots to generate a reduced basis and then perform ordinary least squares operator regression to approximate the Schur complement system. The approach is minimally burdensome, requiring only Gramian matrices capturing a surface integral of solution trial and LM test functions over the interface. We demonstrate the effectiveness of the technique with several numerical experiments and make comparison with respect to speed against Schwarz-based approaches and a full-order, Schur complement based approach.