The Boltzmann Machine
How a 1985 probabilistic network became the foundation of every modern generative model
Introduction
In 1985, Ackley, Hinton, and Sejnowski introduced the Boltzmann machine, a network of binary units that switch between active and inactive states according to probabilities rather than fixed rules. The model learns a probability distribution over observed data by adjusting an energy function defined over network configurations. The name comes from statistical physics and the Boltzmann distribution, the probability law that determines the states of physical systems at thermal equilibrium, where low-energy states are more probable than high-energy ones. Sejnowski recounts the development of the model and its reception across neuroscience, physics, and computer science in The Deep Learning Revolution (Sejnowski, 2018), where it appears as one of the foundational ideas that shaped everything that followed.
The Boltzmann machine’s direct descendants, the restricted Boltzmann machine, the deep belief network, the variational autoencoder, and the score-based diffusion model, each resolved a specific limitation while preserving that objective. The lineage is examined here through a sequence of core questions. Every major generative model in use today traces a direct mathematical lineage to the 1985 paper.
What does the Boltzmann machine’s energy function represent?
The Boltzmann machine is a generative model. Its goal is to learn the probability distribution of the training data well enough to produce new samples that look like they came from the same source. The mathematical structure it uses to represent that distribution is an energy function defined over network configurations. The network consists of binary units, each taking a value of 0 or 1. Some units are visible, corresponding to observed data such as the pixels of an image. Others are hidden, capturing structure in the data that is not directly observed. Every unit is connected to every other unit, with a weight on each connection and a bias on each unit that controls its baseline probability of being active.
The energy function assigns a single number to every possible configuration of the network, where a configuration is a particular assignment of 0s and 1s to all units. For a joint configuration of visible units v and hidden units h, the energy is:
E(v, h) = −Σᵢⱼ wᵢⱼ vᵢ hⱼ − Σ_{i<k} wᵢₖ vᵢ vₖ − Σ_{j<l} wⱼₗ hⱼ hₗ − Σᵢ bᵢ vᵢ − Σⱼ cⱼ hⱼ.
Here wᵢⱼ are weights between visible unit i and hidden unit j, wᵢₖ are weights between visible units i and k, wⱼₗ are weights between hidden units j and l, and bᵢ and cⱼ are biases. Configurations where connected units tend to be simultaneously active produce strongly negative sums and low energy. Configurations where connected units tend to be inactive together, or active in opposing patterns, produce weakly negative or positive sums and high energy.
The Boltzmann distribution is what turns the energy function into a probability distribution over configurations. It assigns a probability to each configuration proportional to the exponential of its negative energy:
p(v, h) = exp(−E(v, h)) / Z.
Here Z is the partition function, a normalizing constant equal to the sum of exp(−E(v, h)) over all possible configurations, whose role is to ensure that all probabilities sum to one. Low energy means high probability. The network learns by adjusting weights so that configurations consistent with the training data receive low energy, and all other configurations receive higher energy. Training succeeds when the probability distribution defined by the energy function closely matches the distribution of the observed data.
A concrete example shows how the energy is computed. Two visible units v₁, v₂ and one hidden unit h have weights w₁ₕ = 1.5, w₂ₕ = 1.5, with visible-visible weights and all biases set to zero. For the configuration
(v₁ = 1, v₂ = 1, h = 1), the energy is −(1.5)(1)(1) − (1.5)(1)(1) = −3. For the configuration (v₁ = 1, v₂ = 0, h = 0), the energy is 0. Applying exp(−E) gives exp(3) ≈ 20.1 for the first configuration and exp(0) = 1 for the second. After normalization, the first configuration has much higher probability. When the training data always shows v₁ = v₂ = 1, this is the distribution the network should learn.
Hidden units give the model its generative power, and understanding their role explains why the architecture works at all. Visible units alone can represent correlations between observed variables, but many real datasets contain dependencies that pairwise statistics of the observed data cannot capture. A hidden unit that activates in the presence of a particular pattern across multiple visible units captures a higher-order dependency compactly. Adding hidden units increases the expressive power of the model, allowing it to capture distributions that no energy function over visible units alone could express. The tradeoff is that hidden units are never observed, so their contribution must be inferred rather than measured directly.
What happens during the wake phase and the sleep phase?
Learning in the Boltzmann machine proceeds through two alternating phases. In the wake phase, the network observes real data and measures how units respond to it. In the sleep phase, the network runs freely without any data input and measures how units behave according to its current weights. The weight update rule compares the two measurements and uses the difference to steer the model distribution toward the data distribution.
The wake phase fixes the visible units to a training example and measures how the hidden units respond. The hidden units then update stochastically, meaning each unit flips between 0 and 1 according to probabilities determined by the weighted inputs from connected units. The updates continue until the network reaches thermal equilibrium, a condition where the probability distribution over hidden unit configurations stops changing. At equilibrium, the expected correlation between each pair of connected units i and j is measured. This quantity is written ⟨vᵢ hⱼ⟩_data and represents the average co-activation of units i and j when the visible units are fixed to training data.
The sleep phase releases all units and measures correlations under the model’s own distribution. Starting from a random configuration, visible and hidden units update freely according to their joint probability under the current model, with no data constraining the visible units. The network runs again until thermal equilibrium is reached. The expected correlation under these conditions is written ⟨vᵢ hⱼ⟩_model.
The weight update rule follows from the difference between the two phases, and this difference is what drives learning. The update is:
ΔWᵢⱼ = ε (⟨vᵢ hⱼ⟩_data − ⟨vᵢ hⱼ⟩_model).
Here ε is the learning rate, a small positive number that controls the size of each update step. When a pair of units is more correlated under the data than under the model, their connection weight increases, lowering the energy of configurations where both are active and raising their probability. When the model produces stronger correlations than the data, the weight decreases. Applied repeatedly across all weights in the network, the rule pushes the model distribution toward the data distribution.
A concrete example shows how the update works. When units v₁ and h₁ co-activate in 80 percent of training examples, ⟨v₁h₁⟩_data = 0.8. When the current model generates this co-activation only 40 percent of the time in the sleep phase, ⟨v₁h₁⟩_model = 0.4. The update is ΔW₁₁ = ε(0.8 − 0.4) = 0.4ε, which increases the weight, lowers the energy of configurations where both units are active, and raises the probability of those configurations under the model. Across many training steps, all weights converge to values where the model correlations match the data correlations, at which point the model distribution closely approximates the data distribution.
The fundamental cost of this procedure is reaching true thermal equilibrium in both phases. In the wake phase, the visible units are fixed, which constrains the hidden unit dynamics and makes equilibrium faster to reach. In the sleep phase, all units are free, and the network must explore its full state space starting from a random configuration. For networks with many units and complex weight structures, this exploration requires a large number of sequential stochastic update steps. A network that runs both phases to true equilibrium and applies the weight update rule repeatedly is guaranteed to converge to the maximum likelihood solution, meaning the set of weights that assigns the highest possible probability to the training data. For networks large enough to model real data, reaching equilibrium in the sleep phase takes prohibitively long, making the original training algorithm impractical for most applications.
What did the restricted Boltzmann machine change about the architecture?
The restricted Boltzmann machine (RBM) introduces a single structural change to the original architecture that transforms the computational difficulty of the problem. All connections within the visible layer and within the hidden layer are removed, leaving connections only between visible units and hidden units. This constraint is what the word “restricted” refers to.
The practical consequence of this restriction is that each layer can be sampled in a single parallel step given the other layer. In the full Boltzmann machine, units within a layer influence each other, creating dependencies that can only be resolved through iterative stochastic updates running to equilibrium. In the RBM, removing within-layer connections creates a structure where visible and hidden units form two separate groups. Given the visible units, all hidden units become conditionally independent of each other. Given the hidden units, all visible units become conditionally independent.
The conditional distributions that result from this structure take a closed form, meaning they can be written as a direct equation and computed without iteration. For hidden unit j given the visible units:
p(hⱼ = 1 | v) = σ(cⱼ + Σᵢ wᵢⱼ vᵢ),
and for visible unit i given the hidden units:
p(vᵢ = 1 | h) = σ(bᵢ + Σⱼ wᵢⱼ hⱼ),
where σ is the logistic sigmoid σ(x) = 1 / (1 + exp(−x)), mapping any real input to a probability between 0 and 1. Both conditionals require only a weighted sum of the opposite layer’s activations followed by a sigmoid.
A concrete example shows how the conditional probabilities are computed. A network with two visible units v₁, v₂ and two hidden units h₁, h₂ has weights w₁₁ = 1.0, w₁₂ = 0.5, w₂₁ = 0.5, w₂₂ = 1.0 and all biases zero. Given the input v = (1, 1), the activation probability for h₁ is σ(1.0 + 0.5) = σ(1.5) ≈ 0.82 and for h₂ is σ(0.5 + 1.0) = σ(1.5) ≈ 0.82. Both hidden units are sampled in one parallel step. Given sampled states h = (1, 1), the reconstruction probabilities for v₁ and v₂ are σ(1.0 + 0.5) ≈ 0.82 and σ(0.5 + 1.0) ≈ 0.82, again computed simultaneously. A complete forward and backward pass takes a fixed number of operations regardless of network size.
The energy function of the RBM simplifies accordingly, and this simplification has consequences for how the model is trained. Without within-layer connections, the terms for visible-visible and hidden-hidden interactions disappear:
E(v, h) = −Σᵢ Σⱼ wᵢⱼ vᵢ hⱼ − Σᵢ bᵢ vᵢ − Σⱼ cⱼ hⱼ.
The probability of a visible configuration v is obtained by summing exp(−E(v, h)) over all hidden configurations. For a hidden layer with n units this sum has 2ⁿ terms, which is intractable for large hidden layers, meaning the computation grows exponentially with the number of hidden units and cannot be completed in a reasonable amount of time. The conditional independence structure of the RBM makes approximate sampling fast and reliable, and it is this property that makes contrastive divergence, the training approximation, both possible and effective.
Why was the original training procedure too slow, and how was it fixed?
Get Eleventh Hour Enthusiast’s stories in your inbox
Join Medium for free to get updates from this writer.
Training the Boltzmann machine requires computing the exact gradient of the log-likelihood, which has a two-phase structure. The positive phase, corresponding to the wake phase, measures correlations under the data distribution and is tractable in the RBM because one forward pass gives the hidden activations given the data. The negative phase, corresponding to the sleep phase, measures correlations under the model distribution. Computing these correlations requires running a Markov chain from random noise all the way to thermal equilibrium, a process requiring many sequential random steps to complete. For large networks this process takes prohibitively long, making exact training impractical. Contrastive divergence is the approximation that solved this problem.
Contrastive divergence replaces the requirement for full equilibrium with a short approximation that starts from the data rather than from random noise. Hinton’s key observation was that the Markov chain does not need to run to equilibrium. Starting from a training example, which is already a point of high probability under the data distribution, and stopping after only k steps produces a sample close enough to the original to capture the gradient direction reliably. Contrastive divergence with k steps is written CD-k, and CD-1 is the most widely used setting in practice.
The CD-1 procedure follows a fixed sequence of steps. Gibbs sampling is the mechanism used, meaning the network alternates between sampling hidden units given visible units and sampling visible units given hidden units. Starting from a training example v⁰, the procedure samples h⁰ from p(h | v⁰), then samples v¹ from p(v | h⁰), then computes h¹ as the expected activation p(h | v¹) without sampling. The weight update becomes:
ΔWᵢⱼ = ε (vᵢ⁰ hⱼ⁰ − vᵢ¹ hⱼ¹),
where h⁰ and h¹ are the expected hidden activations given v⁰ and v¹ respectively. Each update involves two forward passes and one backward pass, and the cost is fixed regardless of network size.
A concrete example shows how the CD-1 steps proceed. A network with input v⁰ = (1, 0) has weights w₁₁ = 0.8, w₁₂ = −0.3, w₂₁ = −0.2, w₂₂ = 0.9 and all biases zero. The hidden activations are h₁⁰ = σ(0.8) ≈ 0.69 and h₂⁰ = σ(−0.2) ≈ 0.45. The reconstruction probabilities are
v₁¹ = σ(0.8 × 0.69 + (−0.2) × 0.45) = σ(0.46) ≈ 0.61 and
v₂¹ = σ((−0.3) × 0.69 + 0.9 × 0.45) = σ(0.20) ≈ 0.55. The reconstructed hidden activations are h₁¹ = σ(0.8 × 0.61 + (−0.3) × 0.55) = σ(0.32) ≈ 0.58 and
h₂¹ = σ((−0.2) × 0.61 + 0.9 × 0.55) = σ(0.37) ≈ 0.59. The update for w₁₁ is
ε(1.0 × 0.69 − 0.61 × 0.58) = ε(0.69 − 0.35) = 0.34ε. The weight increases because v₁ is more active in the data than in the reconstruction, meaning the model is not assigning enough probability to configurations where v₁ and h₁ are simultaneously active.
CD-k does not compute the exact gradient of the log-likelihood, and understanding this limitation explains where the approximation is strong and where it is weak. Early in training, when the model is far from the data distribution and the gradient is large, the approximation is rough but the correct direction is easy to identify. As the model moves closer to the data distribution and the gradient shrinks, the approximation becomes more accurate. The result is a training procedure that is fast, stable, and produces models that assign high probability to the training data.
Why did stacking RBMs revive deep learning, and why was that approach abandoned?
The 2006 paper by Hinton, Osindero, and Teh showed that deep networks could be trained by stacking RBMs one layer at a time, and the procedure had a formal justification that made it more than a heuristic. The approach was greedy and sequential, meaning each layer was trained independently before the next was added. An RBM was first trained on the raw data using contrastive divergence, and its hidden activations were then treated as a new dataset for the next RBM. This process continued for as many layers as desired, with each layer learning a representation of the layer below. The resulting model was called a deep belief network, where the lower layers retain the RBM structure and the top two layers form an undirected associative memory.
The formal justification for greedy training is what separated this approach from earlier attempts at deep networks. Hinton et al. showed that training a new RBM on the hidden representations of the previous one improves a variational lower bound on the log-likelihood of the data, where a variational lower bound is a quantity that is always less than or equal to the true log-likelihood and can be computed without knowing the partition function. The argument is that the lower bound increases each time a new layer is added and trained, because the new layer learns a better prior over the representations produced by the layer below. Adding depth always improved the bound, making the full deep model a better generative model of the training data than any single RBM trained in isolation.
The practical consequence was a solution to the vanishing gradient problem that had made deep networks untrainable with backpropagation alone. In 2006, gradients propagated backward through many sigmoid layers diminished exponentially with depth, leaving early layers unchanged during training. Greedy pretraining avoided this by never requiring gradient propagation across multiple layers at once. Each RBM trained independently with contrastive divergence produced weights that already captured meaningful structure in the data. Fine-tuning the full network with backpropagation from this starting configuration was stable because the pretrained weights placed the network in a region of the loss surface, meaning a region where the training objective responds usefully to small changes in the weights, so that gradients remained informative throughout all layers.
The approach was eventually replaced as better tools for direct end-to-end training became standard, and each tool addressed a specific weakness that greedy pretraining had worked around. Rectified linear units replaced sigmoid activations and eliminated saturation, the condition where a sigmoid unit’s output is very close to 0 or 1, making its gradient nearly zero and blocking learning in earlier layers. Dropout regularization allowed training on limited data without overfitting. Batch normalization stabilized the distribution of activations across layers, making deeper networks trainable from random initialization. Initialization schemes such as Xavier and He initialization gave starting weights scaled according to the number of units in each layer, keeping activations and gradients from growing or shrinking uncontrollably as signals passed through the network. With these advances in place, end-to-end backpropagation from random initialization matched or exceeded greedy pretraining in most settings, and the additional complexity of the layer-by-layer procedure became unnecessary.
The historical importance of the 2006 paper was nonetheless substantial. It demonstrated that deep networks could be trained successfully, and that the obstacle was the training procedure rather than the architecture. That demonstration shifted the field’s attention back to deep networks after a period of limited progress and opened the research directions that led to the tools that eventually made greedy pretraining unnecessary. The idea that useful representations could be learned layer by layer from unlabeled data also outlasted the specific technique and continues to inform approaches to representation learning in settings where labeled data is scarce.
Did wake-sleep training have a formal objective, and what did it become?
The Boltzmann machine’s learning objective is maximum likelihood, meaning finding the weights that assign the highest possible probability to the training data. Computing the log-likelihood
log p(v) = log Σ_h exp(−E(v, h)) − log Z requires the partition function Z, which sums exp(−E(v, h)) over all visible and hidden configurations. For a network with n binary units, this sum has 2ⁿ terms, and since 2ⁿ grows exponentially with n, direct computation becomes impossible for any network large enough to model real data. Wake-sleep training avoids computing Z directly by working with correlations instead, using Gibbs sampling to approximate the gradient. The approximation is practical but is not derived from a formally stated objective function.
The variational autoencoder (VAE), introduced by Kingma and Welling in 2013, arrives at the same generative objective from a different direction and provides the mathematical framework that wake-sleep training lacks. Rather than approximating the intractable gradient of log p(x) directly, the VAE maximizes a lower bound on log p(x) called the evidence lower bound (ELBO). Here x is the observed variable and z is the latent variable, an unobserved internal representation the model uses to describe the structure underlying the data. The approximate posterior q(z | x) is a learned distribution over latent variables given an observation, and p(z) is the prior, the assumed distribution over latent variables before any data is seen. The ELBO is:
ELBO = E_q[log p(x | z)] − KL(q(z | x) || p(z)).
The first term measures reconstruction quality and the second is the KL divergence, a penalty for how far the learned latent distribution departs from the prior. Maximizing the ELBO is equivalent to maximizing a lower bound on log p(x), and the bound is tight when q(z | x) matches the true posterior p(z | x).
Stating both objectives together reveals the connection between wake-sleep and the ELBO. The wake phase fixes visible units to data and infers hidden unit activations. This is the encoder, the network that maps an observation to a distribution over latent variables, computing q(z | x). The sleep phase runs the network freely and generates visible data from sampled latent variables. This is the decoder, the network that maps a latent variable back to the data space, computing p(x | z). The weight update
ΔW = ε(⟨v h⟩_data − ⟨v h⟩_model) approximates the gradient of the same bound the ELBO makes explicit. Wake-sleep optimizes a lower bound on the log-likelihood without a formal framework to identify it as such. The VAE provides that framework.
Three structural changes separate the VAE from the Boltzmann machine, each resolving a specific limitation of the original model. Discrete binary hidden units become a continuous latent space, with the posterior q(z | x) parameterized as a Gaussian with learned mean and variance. The partition function problem disappears because the KL divergence between two Gaussians has a closed form. The reparameterization trick replaces stochastic sampling with a differentiable operation. It expresses z as
z = μ + σ ⊙ ε, where ε is drawn from a standard Gaussian, separating the randomness from the parameters so that gradients can flow through both encoder and decoder during training. Thermal equilibrium is replaced by a closed-form divergence and a differentiable sampling procedure.
The generative adversarial network (GAN), introduced by Goodfellow et al. in 2014, arrives at the same generative goal through a structurally different path. A generator and a discriminator are trained jointly until the generator’s distribution matches the data distribution, using neither an energy function nor a variational bound. GANs produced the first photorealistic generated images and dominated generative modeling for nearly a decade, demonstrating that the problem of learning a data distribution could be solved through competition between two networks rather than through explicit likelihood maximization. The diversity of methods converging on this same objective by 2014 reflects how central the problem of learning the data distribution is.
How does the Boltzmann energy function live on in modern generative models?
Every energy-based model carries a mathematical property that makes the partition function disappear, and this property is the direct link between the 1985 Boltzmann machine and modern diffusion models. Here x refers to a continuous data point rather than a binary configuration, generalizing the energy function beyond the discrete setting of the original Boltzmann machine. For any energy-based model, the probability of a data point x is p(x) = exp(−E(x)) / Z. Taking the gradient of log p(x) with respect to x gives:
∇x log p(x) = −∇x E(x) − ∇x log Z = −∇x E(x).
The partition function Z does not depend on x, so its gradient vanishes. The gradient of the log probability is the negative gradient of the energy. This quantity is called the score function of the distribution, and it points in the direction of steepest increase in log probability from any starting point in the data space. Computing the score requires only the energy function, with no need to compute or approximate Z.
Score-based generative modeling, introduced by Song and Ermon in 2019, takes this gradient as the primary object of learning, and it is the mathematical foundation of what are now widely known as diffusion models. A neural network s_θ(x) is trained to estimate the score function
∇x log p_data(x) of the data distribution. Denoising score matching resolves the problem that the true score requires knowledge of the data distribution by adding Gaussian noise of various magnitudes to training examples and training the network to predict the noise that was added. Predicting the added noise is mathematically equivalent to estimating the score of the noised distribution, and as the noise level decreases toward zero the score converges to that of the original data distribution.
Generation uses the learned score function through Langevin dynamics, a stochastic process that follows the score while adding small controlled noise at each step:
xₜ₊₁ = xₜ + (α/2) sθ(xₜ) + √α εₜ.
Starting from random noise and following the score repeatedly moves the state toward regions of high probability under the data distribution. The Boltzmann machine generated samples by running the network to thermal equilibrium through random unit flips. Diffusion models generate samples by following the score function through a continuous state space, step by step, from noise to data. The mechanism differs. The underlying mathematical object, a surface over the data space whose gradient points toward high probability, is the same one the Boltzmann machine defined in 1985.
Conclusion
The Boltzmann machine defined the problem of learning a probability distribution over observed data and proposed a mathematical structure for solving it. The structure was an energy function whose gradient, when followed, leads from any starting point toward regions of high probability. Training adjusted the energy function so that configurations matching the data received low energy and high probability.
The architecture changed substantially across the four decades that followed. Discrete binary units became continuous latent spaces. Thermal equilibrium was replaced by closed-form divergences and differentiable sampling. Layer-by-layer greedy training gave way to end-to-end backpropagation. The intractable partition function was first approximated, then avoided through variational bounds, and finally made irrelevant by working directly with its gradient. The core objective, learning the distribution well enough to generate from it, did not change. Every generative model in the lineage, the RBM, the deep belief network, the VAE, the GAN, the diffusion model, is each a more scalable and tractable approach to the same problem Ackley, Hinton, and Sejnowski wrote down in 1985.
References:
Ackley, D.H., Hinton, G.E., and Sejnowski, T.J. (1985). A learning algorithm for Boltzmann machines. Cognitive Science, 9(1), 147–169. https://doi.org/10.1207/s15516709cog0901_7
Hinton, G.E., Osindero, S., and Teh, Y.W. (2006). A fast learning algorithm for deep belief nets. Neural Computation, 18(7), 1527–1554. https://doi.org/10.1162/neco.2006.18.7.1527
Kingma, D.P. and Welling, M. (2013). Auto-encoding variational bayes. arXiv e-print, arXiv:1312.6114. https://doi.org/10.48550/arXiv.1312.6114
Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y. (2014). Generative adversarial nets. Advances in Neural Information Processing Systems, 27.
https://doi.org/10.48550/arXiv.1406.2661
Song, Y. and Ermon, S. (2019). Generative modeling by estimating gradients of the data distribution. Advances in Neural Information Processing Systems, 32. https://doi.org/10.48550/arXiv.1907.05600
Sejnowski, T. (2018). The Deep Learning Revolution. MIT Press. https://mitpress.mit.edu/9780262038034/the-deep-learning-revolution/
