  Loop diagram calculations typically rely on reduction to a finite set of
master integrals in $4 - 2\epsilon$ dimensions. It has been shown that for any
problem, the masters can be chosen so that their coefficients are finite as
$\epsilon \rightarrow 0$. I propose a definition of renormalized
$\epsilon$-finite master integrals, which incorporate ultraviolet divergence
subtractions in a specific way. A key advantage of this choice is that in
expressions for physical observables, expansions to positive powers in
$\epsilon$ are never needed. As an example, I provide the subtractions for
general three-loop self-energy integrals. The differential equations method is
used to compute numerically the renormalized $\epsilon$-finite master integrals
for arbitrary external momentum invariant, in special cases with internal
masses equal to a single scale or zero. These include the ones necessary for
the three-loop QCD corrections to the self-energies of the W, Z, and Higgs
bosons. In principle, the same method should provide for numerical computation
of general three-loop self energies with any masses.
