In this paper we propose a novel variance reduction approach for additive functionals of Markov chains based on minimization of an estimate for the asymptotic variance of these functionals over suitable classes of control variates. A distinctive feature of the proposed approach is its ability to significantly reduce the overall finite sample variance. This feature is theoretically demonstrated by means of a deep non asymptotic analysis of a variance reduced functional as well as by a thorough simulation study. In particular we apply our method to various MCMC Bayesian estimation problems where it favourably compares to the existing variance reduction approaches.Assume that P has the unique invariant distribution π, that is, X π(dx)P (x, dy) = π(dy). Under appropriate conditions, the Markov kernel P may be shown to converge to the stationary distribution π, that is, for any x ∈ X, lim n→∞ P n (x, ·) − π TV = 0,