Template fits to observed galaxy fluxes allow calculation of K-corrections and conversions among observations of galaxies at various wavelengths. We present a method for creating model-based template sets given a set of heterogeneous photometric and spectroscopic galaxy data. Our technique, nonnegative matrix factorization, is akin to principal component analysis (PCA), except that it is constrained to produce nonnegative templates, it can use a basis set of models (rather than the delta-function basis of PCA), and it naturally handles uncertainties, missing data, and heterogeneous data (including broadband fluxes at various redshifts). The particular implementation we present here is suitable for ultraviolet, optical, and near-infrared observations in the redshift range 0 < z < 1.5. Since we base our templates on stellar population synthesis models, the results are interpretable in terms of approximate stellar masses and star formation histories. We present templates fitted with this method to data from Galaxy Evolution Explorer, Sloan Digital Sky Survey spectroscopy and photometry, the Two Micron All Sky Survey, the Deep Extragalactic Evolutionary Probe, and the Great Observatories Origins Deep Survey. In addition, we present software for using such data to estimate K-corrections.