We present a formulation of third-order density-functional perturbation theory which is manifestly invariant with respect to unitary transfomations within the occupied-states manifold and is particularly suitable for a practical implementation of the so called2n+ 1'theorem. Our implementation is demonstrated with the calculation of the third-order anharmonic coupling coefficients for some high-simmetry phonons in Silicon.