  We present the largest and most homogeneous collection of near-infrared (NIR)
spectra of Type Ia supernovae (SNe Ia): 339 spectra of 98 individual SNe
obtained as part of the Carnegie Supernova Project-II. These spectra, obtained
with the FIRE spectrograph on the 6.5 m Magellan Baade telescope, have a
spectral range of 0.8--2.5 $\mu$m. Using this sample, we explore the NIR
spectral diversity of SNe Ia and construct a template of spectral time series
as a function of the light-curve-shape parameter, color stretch $s_{BV}$.
Principal component analysis is applied to characterize the diversity of the
spectral features and reduce data dimensionality to a smaller subspace.
Gaussian process regression is then used to model the subspace dependence on
phase and light-curve shape and the associated uncertainty. Our template is
able to predict spectral variations that are correlated with $s_{BV}$, such as
the hallmark NIR features: Mg II at early times and the $H$-band break after
peak. Using this template reduces the systematic uncertainties in K-corrections
by ~90% compared to those from the Hsiao template. These uncertainties, defined
as the mean K-correction differences computed with the color-matched template
and observed spectra, are on the level of $4\times10^{-4}$ mag on average. This
template can serve as the baseline spectral energy distribution for light-curve
fitters and can identify peculiar spectral features that might point to
compelling physics. The results presented here will substantially improve
future SN~Ia cosmological experiments, for both nearby and distant samples.
