Efficient sampling of complex high-dimensional probability densities is a central task in computational science. Machine learning techniques based on autoregressive neural networks provide good approximations to probability distributions of interest in physics. In this work, we propose a systematic way to make this approximation unbiased by using it as an automatic generator of Markov chain Monte Carlo cluster updates. Symmetry enforcing and variable-size cluster updates are found to be essential to the success of this technique. We test our method for first- and second-order phase transitions of classical spin systems, showing its viability for critical systems and in the presence of metastable states.