We propose a control variate method with multiple controls to effectively reduce variances of Monte Carlo simulations for pricing European options under multifactor stochastic volatility models. Based on an application of Ito’s formula, the option price is decomposed by its discounted payoff and stochastic integrals with the appearance of gradients of the unknown option price with respect to state variables (the risky asset price and driving volatility levels). Taking advantage of the closed-form option price approximations obtained by Fouque et al. (SIAM Journal on Multiscale Modeling and Simulation 2(1), 2003), we are able to build controls by substituting approximate option prices into the stochastic integrals. This setup leads to an unbiased control variate and naturally suggests estimates for control parameters. Several numerical experiments are provided to demonstrate the performance of variance reductions by this control variate method. In comparison with variance reduction ratios obtained from importance sampling, we generally find that the control variate method is numerically more stable and efficient than importance sampling for the European option pricing problems under stochastic volatility models.