The calculation of thermal averages plays a central role in equilibrium statistical mechanics. For interacting quantum systems, these averages are challenging to compute due to the high dimensionality and the need to capture quantum effects. Path integral methods offer a route to transform these quantum averages into equivalent classical averages by introducing an extended coordinate space. However, the computational cost grows rapidly with the number of path integral beads. This talk will present two complementary strategies to enable more efficient sampling of quantum thermal averages. First, we investigate the continuum limit as the number of beads approaches infinity, providing algorithms that have exponential ergodicity with respect to the number of beads. Second, we introduce a random batch method that greatly reduces the per-timestep complexity for simulating multiple interacting particles. For both approaches, we establish rigorous error estimates demonstrating improved accuracy and efficiency.