Using the equivalent inclusion method (a method strongly related to the Hashin–Shtrikman variational principle) as a surrogate model, we propose a variance reduction strategy for the numerical homogenization of random composites made of “inclusions” (or rather inhomogeneities) embedded in a homogeneous matrix. The efficiency of this strategy is demonstrated within the framework of two-dimensional, linear conductivity. Significant computational gains vs full-field simulations are obtained even for high contrast values. We also show that our strategy allows to investigate the influence of parameters of the microstructure on the macroscopic response. A specific attention is given to the computational cost of the surrogate model. In particular, an inexpensive approximation of the so-called influence tensors (that are used to compute the surrogate model) is proposed.