An effective and flexible numerical scheme is proposed to simulate the dissipative quantum dynamics of a linearized system-bath Hamiltonian. Based on the observation that the Feynman path integrals for a Gaussian bath have a quadratic functional form, the bath average can be performed by directly sampling paths of the discretized harmonic modes and then propagating the system under the influence of quantum Gaussian force. The algorithm is amenable to all known quantum propagation methods and can thus be flexibly applied to study quantum dissipation in the condensed phase. Nontrivial numerical examples based on the spin-boson and damped quantum oscillator models are presented to demonstrate the application of the new algorithm.