  Efficient sampling of complex high-dimensional probability distributions is a
central task in computational science. Machine learning methods like
autoregressive neural networks, used with Markov chain Monte Carlo sampling,
provide good approximations to such distributions, but suffer from either
intrinsic bias or high variance. In this Letter, we propose a way to make this
approximation unbiased and with low variance. Our method uses physical
symmetries and variable-size cluster updates which utilize the structure of
autoregressive factorization. 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.
