Generative models
The hidden units are just structural support and we don't care about what those hidden vectors really are.
Generative modeling is concerned about $P(X)$, to be able to compute it, we use representation learning (aka deep learning) to identify and disentangle the underlying explanatory factors hidden in the observed data,
because if we could separate them out like this
$$P(X) = P(X \mid Z) P(Z),$$
where $Z$ could encode not just the digit identity but the angle that the digit is drawn, the stroke width, and also abstract stylistic properties.
Any classifier for some task will decide which of those $z$ elements could help him to make a better decision, i.e. just the part that holds the digit identity and ignore the rest, this kind of modeling is called distributed representation, where each element is a simple indicator about some feature.
No conditional dependency within a layer
(my understanding for this reconstructed visible vector is this vector is a vector encoded in the defined RBM in the first place, we are not really construct something new, but we just happen to sample this vector from the defined RBM)
RBM is an undirected graphical model and no conditional dependency within a layer. Thus, the conditional distribution over the hidden units $h$ given the input image $v$ factorizes, that is
$$p(h \mid v) = \prod_i p(h_i \mid v) $$
And the conditional distribution over the visible units $v$ given the hidden units $h$ also factorize.
$$p(v \mid h) = \prod_j p(v_j \mid h) $$
Inference is tractable
Inference in RBM is tractable and we could compute $p(z|x)$, but the tractability of the RBM does not extend to its partition function.
To minimize the negative log-likelihood of the data, consider the following equation.
$$-\frac{\partial \log p(x)}{\partial \theta} = \frac{\partial \mathcal{F}(x)}{\partial \theta}-\sum_{\tilde{x}} p(\tilde{x}) \frac{\partial \mathcal{F}(\tilde{x})}{\partial \theta}$$
The first term $\frac{\partial \mathcal{F}(x)}{\partial \theta}$ increases the probability of training data. The second term $\sum_{\tilde{x}} p(\tilde{x}) \frac{\partial \mathcal{F}(\tilde{x})}{\partial \theta}$ decreases the probability of samples generated by the model.
The second term is nothing less than an expectation over all possible configurations of the input $x$ (under the distribution $P$ formed by the model).
Monte Carlo-based algorithms
So, we use Monte Carlo based algorithms to replace the expectation with an average sum. Samples of $p(x)$ from the model can be obtained by running a Markov chain to convergence, using Gibbs sampling as the transition operator.
For RBMs, the set of visible and hidden units, since they are conditionally independent, one can perform block Gibbs sampling. In this setting, visible units are sampled simultaneously given fixed values of the hidden units. Similarly, hidden units are sampled simultaneously given the visible units.
A step in the Markov chain is thus taken as follows:
- $v$: input values
- $h$: hidden features, latent space, a.k.a $Z$
- $W$: the weight parameters; note that I am using the transposed $W$,
since it's an undirected graphical model.
So, we have the formulas.
$$h^{(n+1)} \sim \operatorname{sigm}\left(W^{\prime} v^{(n)}+c\right)$$
and
$$v^{(n+1)} \sim \operatorname{sigm}\left(W h^{(n+1)}+b\right)$$
where $h^{n+1}$ refers to the set of all hidden units at the $n$th step of the Markov chain. What it means is that, for example, $h^(n+1)_i$ is randomly chosen to be $1$ (versus $0$) with the probability of the sigmoid function and similarly for $v^(n+1)_j$.
In theory, each parameter update in the learning process would require running one such chain to convergence. It is needless to say that doing so would be prohibitively expensive. As such, several algorithms have been devised for RBMs, see Contrastive Divergence.
We initialize the Markov chain with a training example (i.e., from a distribution that is expected to be close to $p$, so that the chain will be already close to having converged to its final distribution $p$).
CD does not wait for the chain to converge. Samples are obtained after only $k$-steps of Gibbs sampling. In practice, $k=1$ has been shown to work surprisingly well.
RBM for dimensionality reduction
RBM can be used to perform dimensionality reduction, and those hidden vectors are some abstract representations of the raw inputs:
According to the paper Reducing the Dimensionality of Data with Neural Networks
In 2006, a breakthrough in feature learning and deep learning took place (Hinton et al., 2006; Bengio et al., 2007; Ranzato et al., 2007), A central idea, referred to as greedy layerwise unsupervised pre-training, was to learn a hierarchy of features one level at a time, using unsupervised feature learning to learn a new transformation at each level
to be composed with the previously learned transformations; essentially, each iteration of unsupervised feature learning adds one layer of weights to a deep neural network. Finally, the set of layers could be combined to initialize a deep supervised predictor, such as a neural network classifier.

Summary
Train single layer RBM to reconstruct to input and then keep adding layer after layer until we reach some good low-level representation as depicted in the figure below, W1, W2, W3 are the matrices that have some lower dimensionality of the input space.
Transpose those matrices $W_1, W_2, W_3$ to be the decoder part.
Fine-tune the whole encoder-decoder stack to reconstruct the input. image

RBM vs VAE and GAN
Side note: RBM is now obsolete, see VAE, GAN they are way better at modeling the manifold in which the data is concentrated, yet you could sample from the manifold directly with a single inference step.