  We present the first surrogate model for gravitational waveforms from the
coalescence of precessing binary black holes. We call this surrogate model
NRSur4d2s. Our methodology significantly extends recently introduced
reduced-order and surrogate modeling techniques, and is capable of directly
modeling numerical relativity waveforms without introducing phenomenological
assumptions or approximations to general relativity. Motivated by GW150914,
LIGO's first detection of gravitational waves from merging black holes, the
model is built from a set of $276$ numerical relativity (NR) simulations with
mass ratios $q \leq 2$, dimensionless spin magnitudes up to $0.8$, and the
restriction that the initial spin of the smaller black hole lies along the axis
of orbital angular momentum. It produces waveforms which begin $\sim 30$
gravitational wave cycles before merger and continue through ringdown, and
which contain the effects of precession as well as all $\ell \in \{2, 3\}$
spin-weighted spherical-harmonic modes. We perform cross-validation studies to
compare the model to NR waveforms \emph{not} used to build the model, and find
a better agreement within the parameter range of the model than other,
state-of-the-art precessing waveform models, with typical mismatches of
$10^{-3}$. We also construct a frequency domain surrogate model (called
NRSur4d2s_FDROM) which can be evaluated in $50\, \mathrm{ms}$ and is suitable
for performing parameter estimation studies on gravitational wave detections
similar to GW150914.
