  We consider generalized quantum Ising models, including those which could
describe disordered materials or quantum annealers, and we prove that for all
temperatures above a system-size independent threshold the path integral Monte
Carlo method based on worldline heat-bath updates always mixes to stationarity
in time $\mathcal{O}(n \log n)$ for an $n$ qubit system, and therefore provides
a fully polynomial-time approximation scheme for the partition function. This
result holds whenever the temperature is greater than four plus twice the
maximum interaction degree (valence) over all qubits, measured in units of the
local coupling strength. For example, this implies that the classical
simulation of the thermal state of a superconducting device modeling a
frustrated quantum Ising model with maximum valence of 6 and coupling strengths
of 1 GHz is always possible at temperatures above 800 mK. Despite the quantum
system being at high temperature, the classical spin system resulting from the
quantum-to-classical mapping contains strong couplings which cause the
single-site Glauber dynamics to mix slowly, therefore this result depends on
the use of worldline updates (which are a form of cluster updates that can be
implemented efficiently). This result places definite constraints on the
temperatures required for a quantum advantage in analog quantum simulation with
various NISQ devices based on equilibrium states of quantum Ising models.
