Magnetic particle images are currently most often reconstructed using classical Tikhonov regularization (i.e. an ℓ2 regularization term) combined with Kaczmarz method. Quality enhancing choices like sparsity promoting ℓ1-regularization or TV regularization lead to problems that cannot be solved by standard Kaczmarz method. We propose to use stochastic primal-dual hybrid gradient method to gain more flexibility concerning the choice of data fitting term and regularization, respectively, and still obtain an algorithm which is at least as fast as Kaczmarz method. The proposed algorithm performs comparably to the current state-of-the-art method in terms of run time. The quality of reconstructions can be significantly improved as different regularization terms can be easily integrated. Moreover, in order to achieve further speed up of the method, we propose two new step size rules which lead to fast convergence and make the algorithm very easy to handle. We improve the performance of the algorithm further by applying a data-driven splitting scheme leading to a significant speed-up during the first iterations.