Abstract Estimates of the generalization error are proved for a residual neural network with $L$ random Fourier features layers $\bar z_{\ell +1}=\bar z_\ell + \textrm {Re}\sum _{k=1}^K\bar b_{\ell k}\,e^{\textrm {i}\omega _{\ell k}\bar z_\ell }+ \textrm {Re}\sum _{k=1}^K\bar c_{\ell k}\,e^{\textrm {i}\omega ^{\prime}_{\ell k}\cdot x}$. An optimal distribution for the frequencies $(\omega _{\ell k},\omega ^{\prime}_{\ell k})$ of the random Fourier features $e^{\textrm {i}\omega _{\ell k}\bar z_\ell }$ and $e^{\textrm {i}\omega ^{\prime}_{\ell k}\cdot x}$ is derived. This derivation is based on the corresponding generalization error for the approximation of the function values $f(x)$. The generalization error turns out to be smaller than the estimate ${\|\hat f\|^2_{L^1({\mathbb {R}}^d)}}/{(KL)}$ of the generalization error for random Fourier features, with one hidden layer and the same total number of nodes $KL$, in the case of the $L^\infty $-norm of $f$ is much less than the $L^1$-norm of its Fourier transform $\hat f$. This understanding of an optimal distribution for random features is used to construct a new training method for a deep residual network. Promising performance of the proposed new algorithm is demonstrated in computational experiments.