We describe ordered subsets (OS) algorithms applied to regularized expectation-maximization (EM) algorithms for emission tomography. Our reconstruction algorithms are based on a maximum a posteriori approach, which allows us to incorporate a priori information in the form of a regularizer to stabilize the unstable EM algorithm. In this work, we use two-dimensional smoothing splines as regularizers. Our motivation for using such regularizers stems from the fact that, by relaxing the requirement of imposing significant spatial discontinuities and using instead quadratic smoothing splines, solutions are easier to compute and hyperparameter calculation becomes less of a problem. To optimize our objective function, we use the method of iterated conditional modes, which is useful for obtaining convenient closed-form solutions. In this case, step sizes or line-search algorithms necessary for gradient-based descent methods are also avoided. We finally accelerate the resulting algorithm using the OS principle and propose a principled way of scaling smoothing parameters to retain the strength of smoothing for different subset numbers. Our experimental results indicate that our new methods provide quantitatively robust results as well as a considerable acceleration.