  We review recent progress on operator mixing in the light of the theory of
canonical forms for linear systems of differential equations and, in
particular, of the Poincar\'e-Dulac theorem. We show that the matrix $A(g) =
-\frac{\gamma(g)}{\beta(g)} =\frac{\gamma_0}{\beta_0}\frac{1}{g} + \cdots $
determines which different cases of operator mixing can occur, and we review
their classification. We derive a sufficient condition for $A(g)$ to be set in
the one-loop exact form $A(g) = \frac{\gamma_0}{\beta_0}\frac{1}{g}$. Finally,
we discuss the consequences of the unitarity requirement in massless QCD-like
theories, and we demonstrate that $\gamma_0$ is always diagonalizable if the
theory is conformal invariant and unitary in its free limit at $g =0$.
