  We introduce the first complete non-parametric model for the astrophysical
distribution of the binary black hole (BBH) population. Constructed from basis
splines, we use these models to conduct the most comprehensive data-driven
investigation of the BBH population to date, simultaneously fitting
non-parametric models for the BBH mass ratio, spin magnitude and misalignment,
and redshift distributions. With GWTC-3, we report the same features previously
recovered with similarly flexible models of the mass distribution, most notably
the peaks in merger rates at primary masses of ${\sim}10\,M_\odot$ and
${\sim}35\,M_\odot$. Our model reports a suppressed merger rate at low primary
masses and a mass ratio distribution consistent with a power law. We infer a
distribution for primary spin misalignments that peaks away from alignment,
supporting conclusions of recent work. We find broad agreement with the
previous inferences of the spin magnitude distribution: the majority of BBH
spins are small ($a<0.5$), the distribution peaks at $a\sim0.2$, and there is
mild support for a non-spinning subpopulation, which may be resolved with
larger catalogs. With a modulated power law describing the BBH merger rate's
evolution in redshift, we see hints of the rate evolution either flattening or
decreasing at $z\sim0.2-0.5$, but the full distribution remains entirely
consistent with a monotonically increasing power law. We conclude with a
discussion of the astrophysical context of our new findings and how
non-parametric methods in gravitational-wave population inference are uniquely
poised to complement to the parametric approach as we enter the data-rich era
of gravitational-wave astronomy.
