How to estimate normalization constant during optimization of complex parameterized PDF using MLE? If I have some data in $\mathbb{R}^N$ space, with $N$ beiling large and I want to estimate the density function of this data, $\mathbf{P}_{data}(x)$, using the fantastic property of neural networks that allows approximation of very complex functions...I will then start with some parameterized formulation $\mathbf{P}_{model}(x|{\theta})$, with ${\theta}$ being the weights in the neural network.
It seems of course I could use the traditional MLE formulation to estimate model parameters, ${\theta}^* = \mathrm{argmax}_{\theta}\prod_{i=1}^m \left[\mathbf{P}_{model}(x|{\theta})\right]$.
The problem I see though is that MLE traditionally relies on a known PDF being used....satisfying the property that $\int_x\mathbf{P}_{model}=1$.  Thus I'm not sure how to do an optimization where the normalizing constant changes throughout the optimization, and/or the integral of $\int_x\mathbf{P}_{model}$ is very hard.  I also can't quite see how the normalizing constant could be pulled in to the backprop equation for weight updates when it is seems like a separate constraint (unless this can be formulated as a constrained optimization problem)?
Doing some reading, I've seen some hints that importance sampling (or some kind of MCMC method) could be used here, but I've not been able to figure out how to implement it, or if its even tractable on a typical deep learning machine (>= 8 cores, 1 high end GPU, >100 GB RAM, etc).
How can I solve this problem...is it even possible?
 A: There are many algorithms for training "unnormalized" models (e.g., score matching, minimum probability flow, noise contrastive estimation, contrastive divergence). In the following, I explain the basic idea behind persistent contrastive divergence (PCD).
Assuming our distribution is positive, we express it in terms of an energy function $E$ and normalization constant $Z$:
$$\textbf{P}_\text{model}(x \mid \theta) = \frac{e^{-E_\theta(x)}}{Z_\theta}$$
The maximum likelihood gradient is
\begin{align}
\nabla \mathbb{E}_\text{data}[\log \textbf{P}_\text{model}(x \mid \theta)]
&= \mathbb{E}_\text{data}[-\nabla E_\theta(x)] - \nabla \log Z_\theta \\
&= \mathbb{E}_\text{model}[\nabla E_\theta(x)] - \mathbb{E}_\text{data}[\nabla E_\theta(x)]
\end{align}
since
$$
\nabla \log Z_\theta 
= \frac{1}{Z_\theta} \nabla \int e^{-E_\theta(x)} \, dx 
= - \int \frac{e^{-E_\theta(x)}}{Z_\theta} \nabla E_\theta(x) \, dx 
= -\mathbb{E}_\text{model}[\nabla E_\theta(x)].
$$
The problem with this gradient is that we cannot easily estimate an expectation over the model distribution, since it is difficult to sample from an unnormalized model. We can try to run an MCMC sampler to sample from the model, but this can take a long time and we have to do it for every gradient step. So in PCD we initialize our sampler with the samples from the last gradient step, which if the model hasn't changed much should already be close to the desired distribution. We then only apply one or a few MCMC updates before using the samples. A commonly used sampler is Hamiltonian Monte Carlo.
One example paper which used deeper nets to represent the energy function and which used PCD is Learning Deep Energy Models (Ngiam et al., ICML, 2011).
A: This is just supplementary information (that won't fit in a comment) after studying the topic much further based on the accepted answer.
I've seen this particular question come up extensively under the area of Restricted Boltzmann Machines (RBMs) with physics analogies like here.
There is still an issue with tracking progress in such a scenario because the partition function $Z_\theta$ still cannot be estimated within reason.  Directly from the link:

Tracking Progress
RBMs are particularly tricky to train. Because of
the partition function Z of Eq. (1), we cannot estimate the
log-likelihood $\log(P(x))$ during training. We therefore have no direct
useful metric for choosing the optimal hyperparameters.
Several options are available to the user.
Inspection of Negative Samples
Negative samples obtained during training can be visualized. As
training progresses, we know that the model defined by the RBM becomes
closer to the true underlying distribution, $p_{train}(x)$. Negative
samples should thus look like samples from the training set. Obviously
bad hyperparameters can be discarded in this fashion.
Visual Inspection of Filters
The filters learnt by the model can be visualized. This amounts to
plotting the weights of each unit as a gray-scale image (after
reshaping to a square matrix). Filters should pick out strong features
in the data. While it is not clear for an arbitrary dataset, what
these features should look like, training on MNIST usually results in
filters which act as stroke detectors, while training on natural
images lead to Gabor like filters if trained in conjunction with a
sparsity criteria.
Proxies to Likelihood
Other, more tractable functions can be used as a proxy to the
likelihood. When training an RBM with PCD, one can use
pseudo-likelihood as the proxy. Pseudo-likelihood (PL) is much less
expensive to compute, as it assumes that all bits are independent.

Others have used annealed importance sampling (AIS) to estimate the partition function directly, although I'm not exactly sure how well it scales to high dimensions it is claimed to be invented for such cases.
DBNs
AIS
