In iterative positron emission tomography (PET) image reconstruction, the statistical variability of the PET data precorrected for random coincidences or acquired in sufficiently high count rates can be properly approximated by a Gaussian distribution, which can lead to a penalized weighted least-squares (PWLS) cost function. In this study, the authors propose a proximal preconditioned gradient algorithm accelerated with ordered subsets (PPG-OS) for the optimization of the PWLS cost function and develop a framework to incorporate boundary side information into edge-preserving total variation (TV) and Huber regularizations.