  The LIGO Scientific, Virgo and KAGRA Collaborations recently released the
third gravitational wave transient catalog or GWTC-3, significantly expanding
the number of gravitational wave (GW) signals. To address the -- still
uncertain -- formation channels of the source compact binaries, their
population properties must be characterized. The computational cost of the
Bayesian hierarchical methods employed thus far scales with the size of the
event catalogs, and such methods have until recently assumed fixed functional
forms for the source distribution. Here we propose a fast and flexible method
to reconstruct the population of LIGO--Virgo merging black hole (BH) binaries
without such assumptions. For sufficiently high event statistics and
sufficiently low individual event measurement error (relative to the scale of
population features) a kernel density estimator (KDE) reconstruction of the
event distribution will be accurate. We improve the accuracy and flexibility of
KDE for finite event statistics using an adaptive bandwidth KDE (awKDE). We
apply awKDE to publicly released parameter estimates for 44 significant (69) BH
binary mergers in GWTC-2 (GWTC-3), in combination with a fast polynomial fit of
search sensitivity, to obtain a non-parametric estimate of the mass
distribution, and compare to Bayesian hierarchical methods. We also demonstrate
a robust peak detection algorithm based on awKDE and use it to calculate the
significance of the apparent peak in the BH mass distribution around $35\,
M_\odot$. We find such a peak is very unlikely to have occurred if the true
distribution is a featureless power-law (significance of $3.6\sigma$ for
confident GWTC-2 BBH events, $3.0\sigma$ for confident GWTC-3 BBH events).
