  We investigate joint modeling of longevity trends using the spatial
statistical framework of Gaussian Process regression. Our analysis is motivated
by the Human Mortality Database (HMD) that provides unified raw mortality
tables for nearly 40 countries. Yet few stochastic models exist for handling
more than two populations at a time. To bridge this gap, we leverage a spatial
covariance framework from machine learning that treats populations as distinct
levels of a factor covariate, explicitly capturing the cross-population
dependence. The proposed multi-output Gaussian Process models straightforwardly
scale up to a dozen populations and moreover intrinsically generate coherent
joint longevity scenarios. In our numerous case studies we investigate
predictive gains from aggregating mortality experience across nations and
genders, including by borrowing the most recently available "foreign" data. We
show that in our approach, information fusion leads to more precise (and
statistically more credible) forecasts. We implement our models in \texttt{R},
as well as a Bayesian version in \texttt{Stan} that provides further
uncertainty quantification regarding the estimated mortality covariance
structure. All examples utilize public HMD datasets.
