In the era of 4th generation synchrotrons, high photon flux and an increase in spatial resolution leads to faster acquisition of large data sets that demands accuracy and computational speed of imaging reconstruction techniques. However, the current computational complexity of most algorithms struggles to meet the demands of big data sizes, requiring robust hardware and advanced numerical techniques. In this manuscript, we present a novel and fast GPU implementation of the Expectation Maximization method, aiming to address this challenge head-on. Our approach capitalizes on frequency domain adaptations of the forward projection and backprojection operators to enhance computational efficiency in tomographic data processing. Moreover, we improve efficiency through parallelization across multiple GPUs. Numerical simulations of large data sets are presented to show the robustness and efficiency of our approach.