Despite decades of developments, combining multiple ground-based astronomy exposures into high signal-to-noise coadds while also improving their spatial resolution is still an outstanding challenge today. Here, we present an expectation-maximization framework for multiframe image denoising and deconvolution, which can be readily extended for performing (1) blind deconvolution, where we self-consistently solve for the point-spread function of each exposure, (2) super-resolution, and (3) improved background subtraction. Our TensorFlow implementation is computationally efficient, scalable, and flexible, as it benefits from advanced algorithmic solutions and leverages Graphical Processing Unit (GPU) acceleration. The testbed for our method is a set of 4K Hyper Suprime-Cam exposures, which are closest to the quality of imaging data from the upcoming Rubin Observatory. Our analysis of the data reveals that performing denoising and deconvolution using traditional gradient descent-based methods tends to result in physically uninterpretable coadds, as the optimizers converge to bad local minima due to their greedy nature. We developed an extension that not only consistently enforces desired constraints, such as non-negativity of pixel values in the coadd, but also proceeds in a way that yields no ``usual'' artifacts. The preliminary results are extremely promising: unprecedented details such as the shape of the spiral arms of galaxies are recovered, while we also manage to deconvolve stars perfectly into essentially single pixels. Statistical tests of the extracted source catalogs are ongoing.