In this paper we address the construction of efficient algorithms for the estimation of gradients of general performance measures of Gaussian systems. Exploiting a clever coupling between the normal and the Maxwell distribution, we present a new gradient estimator, and we show that it outperforms both the single-run based infinitesimal perturbation analysis (IPA) estimator and the score function (SF) estimator, in the one-dimensional case, for a dense class of test functions. Next, we present an example of the multi-dimensional case with a system from the area of stochastic activity networks. Our numerical experiments show that this new estimator also has significantly smaller sample variance than IPA and SF. To increase efficiency, in addition to variance reduction, we present an optimized method for generating the Maxwell distribution, which minimizes the CPU time. © 2007 Elsevier B.V. All rights reserved.