This work introduces a novel hybrid framework that combines operator learning with finite difference discretization to efficiently and accurately solve elliptic interface problems. These problems are characterized by discontinuous coefficients and large jump conditions across complex interfaces. Reformulating the interior-domain flux as an augmented variable denoted by B, we construct an operator \(\mathcal {G}\) that maps B to the solution u of the elliptic problem. This operator is then embedded in a convex optimization problem constrained by the first interface condition (solution jump). Theoretical analysis demonstrates the convexity and uniqueness of the reformulated problem, with error estimates derived under the assumption of Lipschitz continuity. Numerical experiments on five benchmark problems, including smooth/non-smooth coefficients, complex geometries, and high-contrast material properties, show that our method achieves exponential convergence rates, requiring 2-3 orders of magnitude fewer degrees of freedom than traditional mesh-based methods. The proposed approach is novel, as it eliminates the need for interface-fitted meshes while maintaining spectral accuracy, and demonstrates significant potential for multi-physics simulations.
