  We demonstrate that in non-Abelian ${\cal N}=1$ supersymmetric gauge theories
the NSVZ relation is valid for terms quartic in the Yukawa couplings
independently of the subtraction scheme if the renormalization group functions
are defined in terms of the bare couplings and the theory is regularized by
higher covariant derivatives. The terms quartic in the Yukawa couplings appear
in the three-loop $\beta$-function and in the two-loop anomalous dimension of
the matter superfields. We have obtained that the three-loop contribution to
the $\beta$-function quartic in the Yukawa couplings is given by an integral of
double total derivatives. Consequently, one of the loop integrals can be taken
and the three-loop contribution to the $\beta$-function is reduced to the
two-loop contribution to the anomalous dimension. The remaining loop integrals
have been calculated for the simplest form of the higher derivative
regularizing term. Then we construct the renormalization group functions
defined in terms of the renormalized couplings. In the considered approximation
they do not satisfy the NSVZ relation for a general renormalization
prescription. However, we verify that the recently proposed boundary conditions
defining the NSVZ scheme in the non-Abelian case really lead to the NSVZ
relation between the terms of the considered structure.
