  When two black holes merge, a tremendous amount of energy is released in the
form of gravitational radiation in a short span of time, making such events
among the most luminous phenomenon in the universe. Models that predict the
peak luminosity of black hole mergers are of interest to the gravitational wave
community, with potential applications in tests of general relativity. We
present a surrogate model for the peak luminosity that is directly trained on
numerical relativity simulations of precessing binary black holes. Using
Gaussian process regression, we interpolate the peak luminosity in the
7-dimensional parameter space of precessing binaries with mass ratios $q\leq4$,
and spin magnitudes $\chi_1,\chi_2\leq0.8$. We demonstrate that our errors in
estimating the peak luminosity are lower than those of existing fitting
formulae by about an order of magnitude. In addition, we construct a model for
the peak luminosity of aligned-spin binaries with mass ratios $q\leq8$, and
spin magnitudes $|\chi_{1z}|,|\chi_{2z}|\leq0.8$. We apply our precessing model
to infer the peak luminosity of the GW event GW190521, and find the results to
be consistent with previous predictions.
