  We perform statistically rigorous uncertainty quantification (UQ) for chiral
effective field theory ($\chi$EFT) applied to infinite nuclear matter up to
twice nuclear saturation density. The equation of state (EOS) is based on
high-order many-body perturbation theory calculations with nucleon-nucleon and
three-nucleon interactions up to fourth order in the $\chi$EFT expansion. From
these calculations our newly developed Bayesian machine-learning approach
extracts the size and smoothness properties of the correlated EFT truncation
error. We then propose a novel extension that uses multitask machine learning
to reveal correlations between the EOS at different proton fractions. The
inferred in-medium $\chi$EFT breakdown scale in pure neutron matter and
symmetric nuclear matter is consistent with that from free-space
nucleon-nucleon scattering. These significant advances allow us to provide
posterior distributions for the nuclear saturation point and propagate
theoretical uncertainties to derived quantities: the pressure and
incompressibility of symmetric nuclear matter, the nuclear symmetry energy, and
its derivative. Our results, which are validated by statistical diagnostics,
demonstrate that an understanding of truncation-error correlations between
different densities and different observables is crucial for reliable UQ. The
methods developed here are publicly available as annotated Jupyter notebooks.
