We present here a non-perturbative cumulant expansion method for computing the grand partition function of quantum systems. It embeds the classical component of grand partition function exactly and treats the quantum contributions by a systematic cluster expansion of the Boltzmann operator. The cluster expansion Ansatz exploits our recently developed notion of thermal normal ordering and a Wick-like expansion formula, which makes evaluation of the thermal trace particularly easy. The thermal normal ordering also confers a finite expansion structure to the equations for the cluster amplitudes. Our formulation provides manifestly extensive free energy, and works very well over a wide range of temperatures and coupling strengths. As an illustrative application, we have computed the grand partition function of an anharmonic oscillator with a pure quartic perturbation.