  The all-loop resummation of SU$(N)$ gauge theory amplitudes is known to
factorize into an IR-divergent (soft and collinear) factor and a finite (hard)
piece. The divergent factor is universal, whereas the hard function is a
process-dependent quantity. We prove that this factorization persists for the
corresponding celestial amplitudes. Moreover, the soft/collinear factor becomes
a scalar correlator of the product of renormalized Wilson lines defined in
terms of celestial data. Their effect on the hard amplitude is a shift in the
scaling dimensions by an infinite amount, proportional to the cusp anomalous
dimension. This leads us to conclude that the celestial-IR-safe gluon amplitude
corresponds to a expectation value of operators dressed with Wilson line
primaries. These results hold for finite $N$. In the large $N$ limit, we show
that the soft/collinear correlator can be described in terms of vertex
operators in a Coulomb gas of colored scalar primaries with nearest neighbor
interactions. In the particular cases of four and five gluons in planar
$\mathcal{N}=4$ SYM theory, where the hard factor is known to exponentiate, we
establish that the Mellin transform converges in the UV thanks to the fact that
the cusp anomalous dimension is a positive quantity. In other words, the very
existence of the full celestial amplitude is owed to the positivity of the cusp
anomalous dimension.
