Thermalization algorithms for classical gauge theories

I propose a method, based on a set of Langevin equations, for bringing classical gauge theories to thermal equilibrium while respecting the set of Gauss' constraints exactly. The algorithm is described in detail for the SU(2) gauge theory with or without the Higgs doublet. As an example of application, canonical average of the maximal Lyapunov exponent is computed for the SU(2) Yang-Mills theory.