  The interpretation of redshift surveys requires modeling the relationship
between large-scale fluctuations in the observed number density of tracers,
$\delta_\mathrm{h}$, and the underlying matter density, $\delta$. Bias models
often express $\delta_\mathrm{h}$ as a truncated series of integro-differential
operators acting on $\delta$, each weighted by a bias parameter. Due to the
presence of `composite operators' (obtained by multiplying fields evaluated at
the same spatial location), the linear bias parameter measured from clustering
statistics does not coincide with that appearing in the bias expansion. This
issue can be cured by re-writing the expansion in terms of `renormalised'
operators. After providing a pedagogical and comprehensive review of bias
renormalisation in perturbation theory, we generalize the concept to
non-perturbative dynamics and successfully apply it to dark-matter haloes
extracted from a large suite of N-body simulations. When comparing numerical
and perturbative results, we highlight the effect of the window function
employed to smooth the random fields. We then measure the bias parameters as a
function of halo mass by fitting a non-perturbative bias model (both before and
after applying renormalisation) to the cross spectrum
$P_{\delta_\mathrm{h}\delta}(k)$. Finally, we employ Bayesian model selection
to determine the optimal operator set to describe
$P_{\delta_\mathrm{h}\delta}(k)$ for $k<0.2\,h$ Mpc$^{-1}$ at redshift $z=0$.
We find that it includes $\delta, \nabla^2\delta, \delta^2$ and the square of
the traceless tidal tensor, $s^2$. Considering higher-order terms (in $\delta$)
leads to overfitting as they cannot be precisely constrained by our data. We
also notice that next-to-leading-order perturbative solutions are inaccurate
for $k\gtrsim 0.1\,h$ Mpc$^{-1}$.
