### Helpful resources
1. http://www.cs.toronto.edu/~rgrosse/courses/csc321_2017/slides/lec19.pdf
2. http://willwolf.io/2018/10/20/thorough-introduction-to-boltzmann-machines/

### Probabilistic model
A Boltzmann machine is a (fully observable) graphical model whose nodes are binary random variables, i.e., $x_i = \{0, 1\}$ for all $i=1, \dots, n$. A bias $b_i$ determines how likely it is that $x_i=1$ and a weight $w_{ij}$ determines how likely it is that $x_i$ and $x_j$ take the same value. 
#### Model assumption
A Boltzmann machine assumes that the joint probability distribution $p(x) = p(x_1, \dots, x_n)$ can be modeled as 
\begin{align}
p(x) = \frac{e^{H(x)}}{Z},
\end{align}
where $H$ is the happiness function defined as 
\begin{align}
H(x) = \sum_{\substack{i=1 \\ j < i}}^n w_{ij} x_i x_j + \sum_i b_i,
\end{align}
and $Z$ is a normalization constant.

### Learning a Boltzmann machine
Given a dataset $\left\{X^{(j)}\right\}_{j \in [m]}$, we can learn the weights and biases of the model by the principle of maximum likelihood. In contrast to the naive Bayes model, there is no closed-form solution for the parameters, however, we can update the parameters using gradient descent.

#### Log-likelihood function
The log-likelihood function $l$ is given as 
\begin{align}
l = \left[ \frac{1}{m} \sum_{k=1}^m H(X^{(k)} \right] - \log Z,
\end{align}
and it turns out that 
\begin{align}
    \frac{\partial l}{\partial w_{rs}} &= \frac{1}{m} \sum_{k=1}^m X^{(k)}_r X^{(k)}_s - \sum_{x} p(x) x_r x_s, \\
    \frac{\partial l}{\partial b_{r}} &= \frac{1}{m} \sum_{k=1}^m X^{(k)}_r- \sum_{x} p(x) x_r,
\end{align}
where $\sum_x$ is the sum of all possible configurations of the Boltzmann machine.

#### Learning algorithm
We first initialize the weights and biases randomly and can then update them using gradient descent, i.e.,
\begin{align}
    w^{n+1}_{rs} &\gets w^{n}_{rs} + \frac{\partial l}{\partial w^{n}_{rs}}, \\
    b^{n+1}_{r} &\gets b^{n}_{r} + \frac{\partial l}{\partial b^{n}_{r}}.
\end{align}
Note that the sum over all configurations of the Boltzmann machine has $2^n$ terms making exact gradient descent steps infeasible. Hence, we approximate $\sum_{x} p(x) x_r x_s = \mathbb{E}_{\text{model}} [x_r x_s]$ and $\sum_{x} p(x) x_r = \mathbb{E}_{\text{model}} [x_r]$ using a Markov chain Monte Carlo method. For big datasets, we could approximate $\frac{1}{m} \sum_{k=1}^m X^{(k)}_r X^{(k)}_s = \mathbb{E}_{\text{data}} [x_r x_s]$ as well as $\frac{1}{m} \sum_{k=1}^m X^{(k)}_r = \mathbb{E}_{\text{data}} [x_r]$ using mini-batches.

### Sampling
For sampling we use a Markov chain Monte Carlo method called Gibbs sampling. The reason for doing so is that Gibbs sampling is based on the conditional probability distributions
\begin{align}
\text{Pr}\left(x_i = 1 \mid x_{-i}\right),
\end{align}
which can be easily computed with our model assumption; note that $x_{-i} = \{x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_n\}$. It turns out that 
\begin{align}
\text{Pr}\left(x_i = 1 \mid x_{-i}\right) = \sigma \left( \sum_{j \neq i} w_{ij} x_j + b_i \right),
\end{align}
where $\sigma$ is the logistic function
\begin{align}
\sigma(x) = \frac{1}{1 + \exp ( -x )}.
\end{align}
As usual with Markov chain Monte Carlo methods, Gibbs sampling generates a Markov chain of samples. Generally, samples from the beginining of the chain might not represent the true probability distribution well, and therefore we discard the first $n_{\text{burn-in}}$ samples.


```python
import numpy as np
```


```python
class BoltzmannMachine:
    def __init__(self, init_weights, init_biases, configs):
        self.weights = init_weights
        self.biases = init_biases
        self.configs = configs
        
    @staticmethod
    def _logistic_function(x):
        return 1 / (1 + np.exp(-x))
        
    def _happiness_function(self, x):
        happiness = 0
        
        for i in range(self.weights.shape[0]):
            for j in range(i):
                happiness += self.weights[i, j] * x[i] * x[j]
        
        for i, _ in enumerate(self.biases):
            happiness += self.biases[i] * x[i]
                
        return happiness
    
    def likelihood(self, x, compute_log_likelihood=True):
        exponential_happiness = [np.exp(self._happiness_function(x_i)) for x_i in x]        
        
        Z = sum([np.exp(self._happiness_function(config)) for config in self.configs])
        
        if compute_log_likelihood == True:
            return sum([np.log(exp_happ)  - np.log(Z) for exp_happ in exponential_happiness])
        else: 
            return reduce(np.multiply, [exp_happ / Z for exp_happ in exponential_happiness])
        
    def gibbs_sampling(self, n_samples, n_burn_in=1000):
        samples = []
        samples.append([0 for _ in self.biases])
        
        def _gibbs_step(sample, i):
            z = sum([self.weights[i, j] * sample[j] for j in range(len(sample)) if j != i]) + self.biases[i]
            p = self._logistic_function(z)
            return np.random.binomial(n=1, p=p)
        
        for _ in range(n_samples + n_burn_in):
            sample = list(samples[-1])
            for i, _ in enumerate(sample):
                sample[i] = _gibbs_step(sample=sample, i=i)

            samples.append(sample)

        return np.array([sample for i, sample in enumerate(samples[n_burn_in:])])
```

Next, we define two functions to update the weights of the Boltzmann machine. The first function computes $\mathbb{E}_{\text{model}} [x_r x_s]$ and $\mathbb{E}_{\text{model}} [x_r]$ exactly, whereas the second function approximates those quantities using Gibbs sampling.


```python
def gradient_descent(data, model, learning_rate):
    model_distribution = [(np.array(config), model.likelihood(np.array([config]), compute_log_likelihood=False)) \
                           for config in model.configs]
    ### weights
    for i in range(model.weights.shape[0]):
        for j in range(i):
            w_positive = (data[:, i] * data[:, j]).mean()

            w_negative = sum([config[i] * config[j] * likelihood for config, likelihood in model_distribution])

            model.weights[i, j] += learning_rate * (w_positive - w_negative)
    
    ### biases
    for i in range(model.biases.shape[0]):
        b_positive = data[:, i].mean()
        
        b_negative = sum([config[i] * likelihood for config, likelihood in model_distribution])
        
        model.biases[i] += learning_rate * (b_positive - b_negative)
```


```python
def approximate_gradient_descent_using_gibbs_sampling(data, model, learning_rate, n_samples):
    gibbs_samples = model.gibbs_sampling(n_samples)

    ### weights
    for i in range(model.weights.shape[0]):
        for j in range(i):
            w_positive = (data[:, i] * data[:, j]).mean()

            w_negative = (gibbs_samples[:, i] * gibbs_samples[:, j]).mean()

            model.weights[i, j] += learning_rate * (w_positive - w_negative)
    
    ### biases
    for i in range(model.biases.shape[0]):
        b_positive = data[:, i].mean()
        
        b_negative = (gibbs_samples[:, i]).mean()
        
        model.biases[i] += learning_rate * (b_positive - b_negative)
```


```python
from itertools import product
from functools import reduce
from collections import defaultdict

import matplotlib.pyplot as plt
%matplotlib inline

n_units = 3
n_points = 100
n_steps = 1000
learning_rate = 0.1 
p = [.8, .1, .5]

### [[0, 0, 0], [0, 0, 1], [0, 1, 0], ..., [1, 1, 0], [1, 1, 1]]
all_configs = list(product([0, 1], repeat=n_units))

biases = np.random.randn(n_units)
### we initialize weights as a n_units by n_units, but actually just use the elements below the main diagonal
weights = np.random.randn(n_units, n_units)

boltzmann_machine = BoltzmannMachine(weights, biases, all_configs)
```

We are producing some toy data in 3 dimensions.


```python
data =  np.random.binomial(n=1, p=p, size=(n_points, n_units))

data_distribution = np.zeros((len(all_configs),))

for x in data:
    for j, config in enumerate(all_configs):
        if np.sum(x == config) == n_units:
            data_distribution[j] += 1 / n_points
```

We are now training the Boltzmann machine using $n_{\text{steps}}$ steps of gradient descent.


```python
for i in range(n_steps):
    gradient_descent(data, boltzmann_machine, learning_rate)

    if i % 100 == 0:
        log_likelihood = boltzmann_machine.likelihood(data)
        print(f'Epoch: {i:2} | Log-likelihood: {log_likelihood}')
```

    Epoch:  0 | Log-likelihood: -288.0336365242908
    Epoch: 100 | Log-likelihood: -154.3247149791821
    Epoch: 200 | Log-likelihood: -152.98770974403578
    Epoch: 300 | Log-likelihood: -152.66784399515788
    Epoch: 400 | Log-likelihood: -152.4571187752501
    Epoch: 500 | Log-likelihood: -152.30205110294634
    Epoch: 600 | Log-likelihood: -152.1863146368242
    Epoch: 700 | Log-likelihood: -152.09938160789244
    Epoch: 800 | Log-likelihood: -152.03375796159676
    Epoch: 900 | Log-likelihood: -151.98401084176203


We can not print the probability of each configuration learned by the Boltzmann machine.


```python
model_distribution = [(np.array(config), 
                       boltzmann_machine.likelihood(np.array([config]), 
                       compute_log_likelihood=False)) 
                      for config in boltzmann_machine.configs]

for i, _ in enumerate(boltzmann_machine.configs):
    print('Configuration:', all_configs[i], 
          'learned distribution: {0:.3f}'.format(model_distribution[i][1]), 
          'data distribution: {0:.3f}'.format(data_distribution[i]))
```

    Configuration: (0, 0, 0) learned distribution: 0.089 data distribution: 0.100
    Configuration: (0, 0, 1) learned distribution: 0.044 data distribution: 0.040
    Configuration: (0, 1, 0) learned distribution: 0.025 data distribution: 0.020
    Configuration: (0, 1, 1) learned distribution: 0.011 data distribution: 0.010
    Configuration: (1, 0, 0) learned distribution: 0.353 data distribution: 0.350
    Configuration: (1, 0, 1) learned distribution: 0.381 data distribution: 0.380
    Configuration: (1, 1, 0) learned distribution: 0.050 data distribution: 0.050
    Configuration: (1, 1, 1) learned distribution: 0.047 data distribution: 0.050


Lastly, we want to produce samples using Gibbs sampling and compare the sample distribution to the learned distribution/data distribution.


```python
n_samples = 10000
samples = boltzmann_machine.gibbs_sampling(n_samples=10000)

sampling_distribution = np.zeros((len(all_configs),))

for sample in samples:
    for j, config in enumerate(all_configs):
        if np.sum(sample == config) == n_units:
            sampling_distribution[j] += 1 / n_samples
```


```python
for i, _ in enumerate(boltzmann_machine.configs):
    print('Configuration:', all_configs[i], 
          'sampling distribution: {0:.3f}'.format(sampling_distribution[i]))
```

    Configuration: (0, 0, 0) sampling distribution: 0.083
    Configuration: (0, 0, 1) sampling distribution: 0.038
    Configuration: (0, 1, 0) sampling distribution: 0.027
    Configuration: (0, 1, 1) sampling distribution: 0.011
    Configuration: (1, 0, 0) sampling distribution: 0.344
    Configuration: (1, 0, 1) sampling distribution: 0.372
    Configuration: (1, 1, 0) sampling distribution: 0.065
    Configuration: (1, 1, 1) sampling distribution: 0.059


Let us now train a model using the approximate gradient descent method (and the same data).


```python
n_units = 3
n_points = 100
n_steps = 1000
learning_rate = 0.1 
p = np.random.uniform(0, 1, size=(n_units,))

all_configs = list(product([0, 1], repeat=n_units))

biases = np.random.randn(n_units)
weights = np.random.randn(n_units, n_units)

boltzmann_machine = BoltzmannMachine(weights, biases, all_configs)
```


```python
n_samples = 100

for i in range(n_steps):
    approximate_gradient_descent_using_gibbs_sampling(data, boltzmann_machine, learning_rate, n_samples)

    if i % 100 == 0:
        log_likelihood = boltzmann_machine.likelihood(data)
        print(f'Epoch: {i:2} | Log-likelihood: {log_likelihood}')
```

    Epoch:  0 | Log-likelihood: -258.1925192118695
    Epoch: 100 | Log-likelihood: -153.3524918808026
    Epoch: 200 | Log-likelihood: -153.7383024413453
    Epoch: 300 | Log-likelihood: -153.6838978949054
    Epoch: 400 | Log-likelihood: -154.05059012518075
    Epoch: 500 | Log-likelihood: -153.9162049494429
    Epoch: 600 | Log-likelihood: -154.15105041690703
    Epoch: 700 | Log-likelihood: -154.0103952779528
    Epoch: 800 | Log-likelihood: -154.17221830910893
    Epoch: 900 | Log-likelihood: -154.22467600409564



```python
model_distribution = [(np.array(config), 
                       boltzmann_machine.likelihood(np.array([config]), 
                       compute_log_likelihood=False)) 
                      for config in boltzmann_machine.configs]

for i, _ in enumerate(boltzmann_machine.configs):
    print('Configuration:', all_configs[i], 
          'learned distribution: {0:.3f}'.format(model_distribution[i][1]), 
          'data distribution: {0:.3f}'.format(data_distribution[i]))
```

    Configuration: (0, 0, 0) learned distribution: 0.055 data distribution: 0.100
    Configuration: (0, 0, 1) learned distribution: 0.025 data distribution: 0.040
    Configuration: (0, 1, 0) learned distribution: 0.013 data distribution: 0.020
    Configuration: (0, 1, 1) learned distribution: 0.005 data distribution: 0.010
    Configuration: (1, 0, 0) learned distribution: 0.380 data distribution: 0.350
    Configuration: (1, 0, 1) learned distribution: 0.429 data distribution: 0.380
    Configuration: (1, 1, 0) learned distribution: 0.046 data distribution: 0.050
    Configuration: (1, 1, 1) learned distribution: 0.047 data distribution: 0.050


2a1d7570a87f98e5abdaeb8b028dd667c4ea3e36
ipynb
Jupyter Notebook
boltzmann_machine.ipynb
timudk/generative_models_in_numpy
505871be358db047110d2b5b3a854abd81005733
MIT
boltzmann_machine.ipynb
timudk/generative_models_in_numpy
505871be358db047110d2b5b3a854abd81005733
MIT
boltzmann_machine.ipynb
timudk/generative_models_in_numpy
505871be358db047110d2b5b3a854abd81005733
MIT
Qwen/Qwen-72B
1. YES
2. YES