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 similar to 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 x 10(-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.