  We have computed the even-$N$ moments $N \leq 20$ of the gluon-to-quark
splitting function $P_{\rm qg}$ at the fourth order of perturbative QCD via the
renormalization of off-shell operator matrix elements. Our results, derived
analytically for a general gauge group, agree with all results obtained for
this function so far, in particular with the lowest five moments obtained via
physical cross sections. Using our new moments and all available endpoint
constraints, we construct approximations for the four-loop $P_{\rm qg}(x)$ that
should be sufficient for a wide range of collider-physics applications. The
N$^3$LO corrections resulting from these and the corresponding quark-quark
splitting functions lead to a marked improvement of the perturbative accuracy
for the scale derivative of the singlet quark distribution, with effects of 1%
or less at $x \gtrsim 10^{\,-4}$ at a standard reference scale with $\alpha_s =
0.2$.
