The numerical unitarity approach, which allows for the complete reduction of phenomenologically-relevant scattering amplitudes to linear combinations of master integrals, has recently been extended to two loops. In these proceedings we discuss the method and a proof-of-principle calculation of the four-gluon amplitudes in an arbitrary helicity configuration. Furthermore, we combine the method with functional reconstruction techniques to extract analytic formulae for the amplitudes.