  Feynman amplitudes at higher orders in perturbation theory generically have
complex singular structures. Notwithstanding the emergence of many powerful new
methods, the presence of infrared divergences poses significant challenges for
their evaluation. In this article, we develop a systematic method for the
removal of the infrared singularities, by adding appropriate counterterms that
approximate and cancel divergent limits point-by-point at the level of the
integrand. We provide a proof of concept for our method by applying it to
master-integrals that are found in scattering amplitudes for representative
two-to-two scattering processes of massless particles. We demonstrate that,
after the introduction of counterterms, the remainder is finite in four
dimensions. In addition, we find in these cases that the complete singular
dependence of the integrals can be obtained simply by analytically integrating
the counterterms. Finally, we observe that our subtraction method can be also
useful in order to extract in a simple way the asymptotic behavior of Feynman
amplitudes in the limit of small mass parameters.
