  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. Hence,
these searches could lose 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 that exploits
the natural connection between the modes. We use a combination of
post-Newtonian formulae and machine learning tools to model aligned-spin
$(3,3)$, $(4,4)$ waveforms corresponding to a given $(2,2)$ waveform. Each of
these modes can be individually filtered against the data to yield separate
timeseries of signal-to-noise ratios (SNR), which can be combined in a
relatively inexpensive way to marginalize over extrinsic parameters of the
signals. This leads to a HM search pipeline whose matched-filtering cost is
just $\approx 3\times$ that of a quadrupole-only search (in contrast to being
$\approx\! 100 \times$, as in previously proposed HM search methods). Our
method is effectual and is generally applicable for template banks constructed
with either stochastic or geometric placement techniques. Additionally, we
discuss compression of $(2,2)$-only geometric-placement template banks using
machine learning algorithms.
