Purpose: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.Methods:The PPG‐OS algorithm is proposed to address two issues encountered in the optimization of PWLS function with edge‐preserving regularizers. First, the second derivative of this function (Hessian matrix) is shift‐variant and ill‐conditioned due to the weighting matrix (which includes emission data, attenuation, and normalization correction factors) and the regularizer. As a result, the paraboloidal surrogate functions (used in the optimization transfer techniques) end up with high curvatures and gradient‐based algorithms take smaller step‐sizes toward the solution, leading to a slow convergence. In addition, preconditioners used to improve the condition number of the problem, and thus to speed up the convergence, would poorly act on the resulting ill‐conditioned Hessian matrix. Second, the PWLS function with a nondifferentiable penalty such as TV is not amenable to optimization using gradient‐based algorithms. To deal with these issues and also to enhance edge‐preservation of the TV and Huber regularizers by incorporating adaptively or anatomically derived boundary side information, the authors followed a proximal splitting method. Thereby, the optimization of the PWLS function is split into a gradient descent step (upgraded by preconditioning, step size optimization, and ordered subsets) and a proximal mapping associated with boundary weighted TV and Huber regularizers. The proximal mapping is then iteratively solved by dual formulation of the regularizers.Results:The convergence performance of the proposed algorithm was studied with three different diagonal preconditioners and compared with the state‐of‐the‐art separable paraboloidal surrogates accelerated with ordered‐subsets (SPS‐OS) algorithm. In simulation studies using a realistic numerical phantom, it was shown that the proposed algorithm depicts a considerably improved convergence rate over the SPS‐OS algorithm. Furthermore, the results of bias‐variance and signal‐to‐noise evaluations showed that the proposed algorithm with anatomical edge information depicts an improved performance over conventional regularization. Finally, the proposed PPG‐OS algorithm is used for image reconstruction of a clinical study with adaptively derived boundary edge information, demonstrating the potential of the algorithm for fast and edge‐preserving PET image reconstruction.Conclusions:The proposed PPG‐OS algorithm shows an improved convergence rate with the ability of incorporating additional boundary information in regularized PET image reconstruction.