Quantum embedding theories provide a feasible route for obtaining quantitative descriptions of correlated materials. However, a critical challenge is solving an effective impurity model of correlated orbitals embedded in an electron bath. Many advanced impurity solvers require the approximation of a bath continuum using a finite number of bath levels, producing a highly nonconvex, ill-conditioned inverse problem. To address this drawback, this study proposes an efficient fitting algorithm for matrix-valued hybridization functions based on a data-science approach, sparse modeling, and a compact representation of Matsubara Greens functions. The efficiency of the proposed method is demonstrated by fitting random hybridization functions with large off-diagonal elements as well as those of a 20-orbital impurity model for a high-Tc compound, LaAsFeO, at low temperatures (T). The results set quantitative goals for the future development of impurity solvers toward quantum embedding simulations of complex correlated materials.