We present a data-driven surrogate modeling for parameterized electromagnetic simulation. This method extracts a set of reduced basis (RB) functions from full-order solutions through a two-step proper orthogonal decomposition (POD) method. A mapping from the time/parameter to the principal components of the projection coefficients, extracted by canonical polyadic decomposition (CPD) , is approximated by a cubic spline interpolation (CSI) approach. The reduced-order model (ROM) is trained in the offline phase, while the RB solution of a new time/parameter value is recovered fast during the online phase. We evaluate the performance of the proposed method with numerical tests for the scattering of a plane wave by a 2-D multi-layer dielectric disk and a 3-D multi-layer dielectric sphere.