  The computation of two-loop amplitudes for the production of multiple Higgs
and electroweak gauge bosons via gluon fusion with exact dependence on quark
masses relies primarily on numerical methods. We propose a framework that
enables their numerical evaluation in momentum space. The method is inspired by
the factorization of infrared divergences in QCD scattering amplitudes. It
extends techniques introduced for electroweak gauge boson production from
quark-antiquark annihilation to processes with external gluons. By combining
diagrammatic integrands, we make use of local cancellations between diagrams
that automatically eliminate most non-factoring infrared singularities. With a
limited number of counterterms, we then derive two-loop integrands for which
all soft and collinear singularities factorize locally. We hope that the local
subtraction techniques presented in this article will play a useful role in
extending the local factorization formalism to two-loop amplitudes for
arbitrary processes.
