We speed up thermal simulations of quantum many-body systems in both one- (1D) and two-dimensional (2D) models in an exponential way by iteratively projecting the thermal density matrix $\hat\rho=e^{-\beta \hat{H}}$ onto itself. We refer to this scheme of doubling $\beta$ in each step of the imaginary time evolution as the exponential tensor renormalization group (XTRG). This approach is in stark contrast to conventional Trotter-Suzuki-type methods which evolve $\hat\rho$ on a linear quasi-continuous grid in inverse temperature $\beta \equiv 1/T$. In general, XTRG can reach low temperatures exponentially fast, and thus not only saves computational time but also merits better accuracy due to significantly fewer truncation steps. We work in an (effective) 1D setting exploiting matrix product operators (MPOs) which allows us to fully and uniquely implement non-Abelian and Abelian symmetries to greatly enhance numerical performance. We use our XTRG machinery to explore the thermal properties of Heisenberg models on 1D chains and 2D square and triangular lattices down to low temperatures approaching ground state properties. The entanglement properties, as well as the renormalization group flow of entanglement spectra in MPOs, are discussed, where logarithmic entropies (approximately $\ln\beta$) are shown in both spin chains and square lattice models with gapless towers of states. We also reveal that XTRG can be employed to accurately simulate the Heisenberg XXZ model on the square lattice which undergoes a thermal phase transition. We determine its critical temperature based on thermal physical observables, as well as entanglement measures. Overall, we demonstrate that XTRG provides an elegant, versatile, and highly competitive approach to explore thermal properties in both 1D and 2D quantum lattice models.