  We introduce a new method for calculating the mixing matrix for non-singlet
quark operators including total derivatives, based solely on their
renormalization structure in the chiral limit. As input, the method requires
the well-known forward anomalous dimensions, which determine the evolution of
parton distribution functions, and a calculation of the matrix elements of
operators without total derivatives. Assuming a large number of quark flavors
$n_f$, we are able to calculate the mixing matrix to fifth order in the strong
coupling $\alpha_s$ in the $\overline{\mbox{MS}}$-scheme.
