Is estimation partition functions without resorting to markov chain monte carlo still an open question? I was told estimation of partition functions without resorting to MCMC was still an open question in physics about a year and a half ago.  An example is that say you have some physical model that obeys a Boltzmann distribution: $p(x) \sim \frac {e^{f(x)}}{Z(x)}$.  The objective is to calculate $Z(x)$, without resorting to some form of rejection sampling.  I don't think partition function estimation was this person's field so he may have not been fully informed, but he was a professor of physics.  I'm curious because this paper is one of many that have dealt with and may have solved this issue in machine learning if I'm understanding this right from a physics perspective (which I may not be). 
In order to help me understand this question, it may be helpful to explain what issues there are with this method for estimating partition functions that are needed in physics if it still is an open question.  If it is not an open question, perhaps it may be interesting to point to the papers that have solved this if they were proposed recently.   
PS That paper also owes a debt of gratitude for physicists for introducing Russian roulette estimator.  
 A: I am going to discuss my knowledge of sampling methods and their application to partition function estimation from a physical point of view. I am not a statistician, so I can't provide any meaningful input on the paper in your link: maybe some statisiticians could help with that?
For a reduced potential energy function $u(\vec{x})$ the configurational partition function over a configurational volume $V$ is (barring any proportionality constants):
$$Z_1 = \int_{V} e^{-u(\vec{x})} d\vec{x}$$
Now, as you have already seen, this is a $3N$-dimensional integral which can only be estimated by probabilistic methods (quadrature results in a combinatorial explosion). The problem with the exponential term is that an overwhelming number of sample points will have a negligible Boltzmann factor, while the terms that contribute the most (the low energy configurations) are extremely rare. Therefore, the variance of any integral estimate is going to be extremely high and not practically useful. That's why your best bet is to convert this into a sampling problem.
The way to do this is to start with a distribution whose partition function you know (usually a uniform distribution, in which case $Z_0=V^{N}$) and then use any method which can estimate ratios of normalisation constants and/or free energy differences $\Delta f_{01}$, since $\frac{Z_1}{Z_0} = e^{-\Delta f_{01}}$. There are many methods which can estimate free energy differences, but all of them rely on intermediate connecting distributions $Z_{\lambda},\lambda\in (0,1)$. However, in order to sample from these you need global sampling, which with regular sampling methods is extremely inefficient (the efficiency decreases exponentially with the kinetic barrier heights, i.e. roughly with system size).
Alternatively, you can recast the partition function as a one-dimensional integral:
$$Z_1 = \int_{-\infty}^{\infty} e^{-u}\Omega(u)du$$
for a density of states $\Omega(u)$:
$$\Omega(u) = \int_{V} \delta(u(\vec{x})-u) e^{-(u(\vec{x})-u)} d\vec{x}$$
which means that you "only" need to estimate the density of states and you are done! Of course, this is the difficult bit, but there is a method which claims to do that even for high-dimensional systems: nested sampling (NS). I am not going to go into detail, you can read about it if you are interested, but suffice it to say that this method results in this impressive paper, where they calculate a whole phase diagram based on estimating the density of states.
Another method that provides the level of sampling needed to estimate the density of states is sequential Monte Carlo (SMC). In fact, statisticians much prefer SMC to NS and it has been suggested that SMC is indeed a generalisation of NS. SMC is a path-integral-like method which also traverses intermediate distributions from the uniform distribution to the distribution of interest, but it uses a survival of the fittest mechanism, which only results in good trajectories through parameter space.
In all cases, estimating partition functions is still very difficult and the scaling is highly polynomial at best (but practically still likely exponential), so I would say this is probably still an open problem, but it is entirely conceivable that it is impossible to come up with an algorithm with a lower computational complexity than what we currently have - only time will tell.
