We develop a gauge covariant neural network for four dimensional non-abelian gauge theory, which realizes a map between rank-2 tensor valued vector fields.We find that the conventional smearing procedure and gradient flow for gauge fields can be regarded as known neural networks, residual networks and neural ordinal differential equations for rank-2 tensors with fixed parameters.In terms of machine learning context, projection or normalization functions in the smearing schemes correspond to an activation function in neural networks.Using the locality of the activation function, we derive the backpropagation for the gauge covariant neural network. Consequently, the smeared force in hybrid Monte Carlo (HMC) is naturally derived with the backpropagation. As a demonstration, we develop the self-learning HMC (SLHMC) with covariant neural network approximated action for non-abelian gauge theory with dynamical fermions, and we observe SLHMC reproduces results from HMC.