Searches for gravitational wave events use models, or templates, for the signals of interest. The templates used in current searches in the LIGO-Virgo-Kagra (LVK) data model the dominant quadrupole mode $(\ell,|m|)=(2,2)$ of the signals, and omit sub-dominant higher-order modes (HM) such as $(\ell,|m|)=(3,3)$, $(4,4)$, which are predicted by general relativity. This omission reduces search sensitivity to black hole mergers in interesting parts of parameter space, such as systems with high masses and asymmetric mass-ratios. We develop a new strategy to include HM in template banks: instead of making templates containing a combination of different modes, we separately store normalized templates corresponding to $(2,2)$, $(3,3)$ and $(4,4)$ modes. To model aligned-spin $(3,3)$, $(4,4)$ waveforms corresponding to a given $(2,2)$ waveform, we use a combination of post-Newtonian formulae and machine learning tools. In the matched filtering stage, one can filter each mode separately with the data and collect the timeseries of signal-to-noise ratios (SNR). This leads to a HM template bank whose matched-filtering cost is just $\approx 3\times$ that of a quadrupole-only search (as opposed to $\approx\! 100 \times$ in previously proposed HM search methods). Our method is effectual and generally applicable for template banks constructed with either stochastic or geometric placement techniques. New GW candidate events that we detect using our HM banks and details for combining the different SNR mode timeseries are presented in accompanying papers: Wadekar et al. [1] and [2] respectively. Additionally, we discuss non-linear compression of $(2,2)$-only geometric-placement template banks using machine learning algorithms.