  In order to extract cosmological information from observations of the
millimeter and submillimeter sky, foreground components must first be removed
to produce an estimate of the cosmic microwave background (CMB). We developed a
machine-learning approach for doing so for full-sky temperature maps of the
millimeter and submillimeter sky. We constructed a Bayesian spherical
convolutional neural network architecture to produce a model that captures both
spectral and morphological aspects of the foregrounds. Additionally, the model
outputs a per-pixel error estimate that incorporates both statistical and model
uncertainties. The model was then trained using simulations that incorporated
knowledge of these foreground components that was available at the time of the
launch of the Planck satellite. On simulated maps, the CMB is recovered with a
mean absolute difference of $<4\mu$K over the full sky after masking map pixels
with a predicted standard error of $>50\mu$K; the angular power spectrum is
also accurately recovered. Once validated with the simulations, this model was
applied to Planck temperature observations from its 70GHz through 857GHz
channels to produce a foreground-cleaned CMB map at a Healpix map resolution of
NSIDE=512. Furthermore, we demonstrate the utility of the technique for
evaluating how well different simulations match observations, particularly in
regard to the modeling of thermal dust.
