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 seven-dimensional parameter space of precessing binaries with mass ratios q ≤ 4 and spin magnitudes χ₁, χ₂ ≤ 0.8. We demonstrate that our errors in estimating the peak luminosity are lower than those of existing fitting formulas by about an order of magnitude. In addition, we construct a model for the peak luminosity of aligned-spin binaries with mass ratios q ≤ 8 and spin magnitudes |χ₁_z|,|χ₂_z| ≤ 0.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
