In this talk, we discuss an ongoing calculation of the four-loop form factors in massless QCD. We begin by discussing our novel approach to the calculation in detail. Of particular interest are a new polynomial-time integration by parts reduction algorithm and a new method to algebraically resolve the IR and UV singularities of dimensionally-regulated bare perturbative scattering amplitudes. Although not all integral topologies are linearly reducible for the more non-trivial color structures, it is nevertheless feasible to obtain accurate numerical results for the finite parts of the complete four-loop form factors using publicly available sector decomposition programs and bases of finite integrals. Finally, we present first results for the four-loop gluon form factor Feynman diagrams which contain three closed fermion loops.