# Restricted Boltzmann Machine

The restricted Boltzman Machine model is the Joint Probability Distribution which is specified by the Energy Function :

\begin{equation}
P(v,h) = \frac{1}{Z} e^{-E(v,h)}
\end{equation}

The energy function for the RBM is stated as follows:

\begin{equation}
E(v,h) = -b^{T} v - c^{T} h - v^{T} W h 
\end{equation}

We also have the Partition Function Z which is the normalizing constant.

\begin{equation}
Z = \Sigma_{v} \Sigma_{h} e^{-E(v,h)}
\end{equation}


In Boltzmann Machines, the partition function Z is intractable and hence implies that the normalized Joint Probability Distribution _P(v)_ is also intractable to evaluate. Even though this is the case, the bipartitie graph structure of the RBM has a special property that the visible and hidden units are conditionally independent, given one another. 

\begin{equation}
P(h|v) = \frac{P(h,v)}{P(v)}
       = \frac{1}{P(v)} \frac{1}{Z} exp\{b^{T} v + c^{T} h + v^{T} W h\}
\end{equation}

\begin{equation}
    = \frac{1}{Z'} exp\{\Sigma_j c_j h_j + \Sigma_j v^T W_j h_j \}
\end{equation}

\begin{equation}
    = \frac{1}{Z'} \Pi_j exp \{ c_j h_j + v^T W_j h_j\}
\end{equation}

\begin{equation}
P(h_j = 1,v) = \frac{\hat{P}(h_j = 1,v)}{\hat{P}(h_j = 0,v) + \hat{P}(h_j = 1,v)}
             = \frac{exp\{c_j + v^T W_j \}}{exp\{0\} + exp\{c_j + v^T W_j \}}
\end{equation}

Therefore:
\begin{equation}
P(h_j = 1|v) = \sigma(c_j + v^T W_j )
\end{equation}

Similary from Eq 8 we can say that:
\begin{equation}
P(v|h) = \frac{1}{Z'} \Pi_k exp \{b_k + h^T W_k \}
\end{equation}
Therefore:
\begin{equation}
P(v_k = 1|h) = \sigma(b_k + h^T W_k)
\end{equation}


```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import os
from tensorflow.examples.tutorials.mnist import input_data


if not os.path.exists('outRBM/'):
    os.makedirs('outRBM/')

mnist = input_data.read_data_sets('../MNIST_data', one_hot=True)
X_dim = mnist.train.images.shape[1]
y_dim = mnist.train.labels.shape[1]

mb_size = 16
h_dim = 20

W = np.random.randn(X_dim, h_dim) * 0.001
a = np.random.randn(h_dim) * 0.001
b = np.random.randn(X_dim) * 0.001


def sigm(x):
    return 1/(1 + np.exp(-x))


def infer(X):
    # mb_size x x_dim -> mb_size x h_dim
    return sigm(X @ W)


def generate(H):
    # mb_size x h_dim -> mb_size x x_dim
    return sigm(H @ W.T)


# Here we find the Contrastive Divergence
# ----------------------
# Approximate the log partition gradient Gibbs sampling

alpha = 0.1
K = 15  # Num. of Gibbs sampling step

for t in range(1, 1001):
    X_mb = (mnist.train.next_batch(mb_size)[0] > 0.5).astype(np.float)
    g = 0
    g_a = 0
    g_b = 0

    for v in X_mb:
        # E[h|v,W]
        h = infer(v)

        # Gibbs sampling steps
        # --------------------
        v_prime = np.copy(v)

        for k in range(K):
            # h ~ p(h|v,W)
            h_prime = np.random.binomial(n=1, p=infer(v_prime))
            # v ~ p(v|h,W)
            v_prime = np.random.binomial(n=1, p=generate(h_prime))

        # E[h|v',W]
        h_prime = infer(v_prime)

        # Compute data gradient
        grad_w = np.outer(v, h) - np.outer(v_prime, h_prime)
        grad_a = h - h_prime
        grad_b = v - v_prime

        # Accumulate minibatch gradient
        g += grad_w
        g_a += grad_a
        g_b += grad_b

    # Monte carlo gradient
    g *= 1 / mb_size
    g_a *= 1 / mb_size
    g_b *= 1 / mb_size

    # Update to maximize
    W += alpha * g
    a += alpha * g_a
    b += alpha * g_b


# Visualization
# -------------

def plot(samples, size, name):
    size = int(size)
    fig = plt.figure(figsize=(4, 4))
    gs = gridspec.GridSpec(4, 4)
    gs.update(wspace=0.05, hspace=0.05)

    for i, sample in enumerate(samples):
        ax = plt.subplot(gs[i])
        plt.axis('off')
        ax.set_xticklabels([])
        ax.set_yticklabels([])
        ax.set_aspect('equal')
        plt.imshow(sample.reshape(size, size), cmap='Greys_r')

    plt.savefig('outRBM/{}.png'.format(name), bbox_inches='tight')
    plt.close(fig)


X = (mnist.test.next_batch(mb_size)[0] > 0.5).astype(np.float)

H = np.random.binomial(n=1, p=infer(X))
plot(H, np.sqrt(h_dim), 'H')

X_recon = (generate(H) > 0.5).astype(np.float)
plot(X_recon, np.sqrt(X_dim), 'V')
```

# Variational Autoencoder 

In Variational AutoEncoder we convert the input into a latent space and calculate the mean and standard deviation. By combining this, we get a new distribution of latent space which provides a better output than a regular autoencoder. 

We shall consider the following definitions when deriving the VAE.

    1. X      - The data to be modeled
    2. z      - Latent Variable
    3. P(X)   - Probability distribution of Data
    4. P(z)   - Probability Distribution of Latent Variable
    P(X|z)    - Probability Distribution of the generated Data given the Latent Variable.

In VAE, we try to find the latent space _z_ using the data _X_. Hence we try to infer _P(z)_ using _P(z|X)_. But we do not know _P(z|X)_, hence we use Variational Inference and approach it as an optimization problem. We do this by actual distribution _P(z|X)_ using a simpler distribution such as Gaussian and then find the difference between the two distributions using KL Divergence.

For inferring _P(z|X)_ using _Q(z|X)_, we have the KL Divergence as :
\begin{equation}
D_{KL}[Q(z|X)∥P(z|X)] = \Sigma_z Q(z|X) log \frac{Q(z|X)}{P(z|X)}
                      = E[log Q(z|X) - log P(z|X)]
\end{equation}



Using Bayes' Rule, we can expand _P(z|X)_ as the following:

\begin{equation}
D_{KL}[Q(z|X)∥P(z|X)] = E[log Q(z|X) - log \frac{P(X|z) P(z)}{P(X)}]
= E[log Q(z|X) − log P(X|z) − log P(z) + log P(X)]
\end{equation}

The expection we are finding is over _z_, hence independent of _x_. Therefore, we we move _x_ out of the expectation. 

\begin{equation}
D_{KL}[Q(z|X)∥P(z|X)] = E[log Q(z|X) − log P(X|z) − log P(z) ] + log P(X)
\end{equation}
\begin{equation}
D_{KL}[Q(z|X)∥P(z|X)] - log P(X)  = E[log Q(z|X) − log P(X|z) − log P(z)]
\end{equation}

This equation can then be written as another KL Divergence. That is shown below. 

\begin{equation}
log P(X) - D_{KL}[Q(z|X)∥P(z|X)] = E[log P(X|z) - (log Q(z|X) − log P(z))]
= E[log P(X|z)] - E[log Q(z|X) - log P(z)]
\end{equation}


This is the final objective functions that we arrive from the derivation:

\begin{equation}
log P(X) - D_{KL}[Q(z|X)∥P(z|X)] = E[log P(X|z)] - D_{KL}[Q(z|X)∥P(z)]
\end{equation}

## Vanilla VAE



```python
import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import os
from tensorflow.examples.tutorials.mnist import input_data


mnist = input_data.read_data_sets('../MNIST_data', one_hot=True)
mb_size = 64
z_dim = 
X_dim = mnist.train.images.shape[1]
y_dim = mnist.train.labels.shape[1]
h_dim = 256
c = 0
lr = 1e-3

def plot(samples):
    fig = plt.figure(figsize=(4, 4))
    gs = gridspec.GridSpec(4, 4)
    gs.update(wspace=0.05, hspace=0.05)

    for i, sample in enumerate(samples):
        ax = plt.subplot(gs[i])
        plt.axis('off')
        ax.set_xticklabels([])
        ax.set_yticklabels([])
        ax.set_aspect('equal')
        plt.imshow(sample.reshape(28, 28), cmap='Greys_r')

    return fig


def xavier_init(size):
    in_dim = size[0]
    xavier_stddev = 1. / tf.sqrt(in_dim / 2.)
    return tf.random_normal(shape=size, stddev=xavier_stddev)


# =============================== Q(z|X) ======================================

X = tf.placeholder(tf.float32, shape=[None, X_dim])
z = tf.placeholder(tf.float32, shape=[None, z_dim])

Q_W1 = tf.Variable(xavier_init([X_dim, h_dim]))
Q_b1 = tf.Variable(tf.zeros(shape=[h_dim]))

Q_W2_mu = tf.Variable(xavier_init([h_dim, z_dim]))
Q_b2_mu = tf.Variable(tf.zeros(shape=[z_dim]))

Q_W2_sigma = tf.Variable(xavier_init([h_dim, z_dim]))
Q_b2_sigma = tf.Variable(tf.zeros(shape=[z_dim]))


def Q(X):
    h = tf.nn.relu(tf.matmul(X, Q_W1) + Q_b1)
    z_mu = tf.matmul(h, Q_W2_mu) + Q_b2_mu
    z_logvar = tf.matmul(h, Q_W2_sigma) + Q_b2_sigma
    return z_mu, z_logvar


def sample_z(mu, log_var):
    eps = tf.random_normal(shape=tf.shape(mu))
    return mu + tf.exp(log_var / 2) * eps


# =============================== P(X|z) ======================================

P_W1 = tf.Variable(xavier_init([z_dim, h_dim]))
P_b1 = tf.Variable(tf.zeros(shape=[h_dim]))

P_W2 = tf.Variable(xavier_init([h_dim, X_dim]))
P_b2 = tf.Variable(tf.zeros(shape=[X_dim]))


def P(z):
    h = tf.nn.relu(tf.matmul(z, P_W1) + P_b1)
    logits = tf.matmul(h, P_W2) + P_b2
    prob = tf.nn.sigmoid(logits)
    return prob, logits


# =============================== TRAINING ====================================

z_mu, z_logvar = Q(X)
z_sample = sample_z(z_mu, z_logvar)
_, logits = P(z_sample)

# Sampling from random z
X_samples, _ = P(z)

# E[log P(X|z)]
recon_loss = tf.reduce_sum(tf.nn.sigmoid_cross_entropy_with_logits(logits=logits, labels=X), 1)
# D_KL(Q(z|X) || P(z)); calculate in closed form as both dist. are Gaussian
kl_loss = 0.5 * tf.reduce_sum(tf.exp(z_logvar) + z_mu**2 - 1. - z_logvar, 1)
# VAE loss
vae_loss = tf.reduce_mean(recon_loss + kl_loss)

solver = tf.train.AdamOptimizer().minimize(vae_loss)

sess = tf.Session()
sess.run(tf.global_variables_initializer())
```

    WARNING:tensorflow:From <ipython-input-2-7e44334aa31a>:9: read_data_sets (from tensorflow.contrib.learn.python.learn.datasets.mnist) is deprecated and will be removed in a future version.
    Instructions for updating:
    Please use alternatives such as official/mnist/dataset.py from tensorflow/models.
    WARNING:tensorflow:From /Users/akash/anaconda3/lib/python3.7/site-packages/tensorflow/contrib/learn/python/learn/datasets/mnist.py:260: maybe_download (from tensorflow.contrib.learn.python.learn.datasets.base) is deprecated and will be removed in a future version.
    Instructions for updating:
    Please write your own downloading logic.
    WARNING:tensorflow:From /Users/akash/anaconda3/lib/python3.7/site-packages/tensorflow/contrib/learn/python/learn/datasets/mnist.py:262: extract_images (from tensorflow.contrib.learn.python.learn.datasets.mnist) is deprecated and will be removed in a future version.
    Instructions for updating:
    Please use tf.data to implement this functionality.
    Extracting ../MNIST_data/train-images-idx3-ubyte.gz
    WARNING:tensorflow:From /Users/akash/anaconda3/lib/python3.7/site-packages/tensorflow/contrib/learn/python/learn/datasets/mnist.py:267: extract_labels (from tensorflow.contrib.learn.python.learn.datasets.mnist) is deprecated and will be removed in a future version.
    Instructions for updating:
    Please use tf.data to implement this functionality.
    Extracting ../MNIST_data/train-labels-idx1-ubyte.gz
    WARNING:tensorflow:From /Users/akash/anaconda3/lib/python3.7/site-packages/tensorflow/contrib/learn/python/learn/datasets/mnist.py:110: dense_to_one_hot (from tensorflow.contrib.learn.python.learn.datasets.mnist) is deprecated and will be removed in a future version.
    Instructions for updating:
    Please use tf.one_hot on tensors.
    Extracting ../MNIST_data/t10k-images-idx3-ubyte.gz
    Extracting ../MNIST_data/t10k-labels-idx1-ubyte.gz
    WARNING:tensorflow:From /Users/akash/anaconda3/lib/python3.7/site-packages/tensorflow/contrib/learn/python/learn/datasets/mnist.py:290: DataSet.__init__ (from tensorflow.contrib.learn.python.learn.datasets.mnist) is deprecated and will be removed in a future version.
    Instructions for updating:
    Please use alternatives such as official/mnist/dataset.py from tensorflow/models.
    WARNING:tensorflow:From /Users/akash/anaconda3/lib/python3.7/site-packages/tensorflow/python/framework/op_def_library.py:263: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.
    Instructions for updating:
    Colocations handled automatically by placer.
    WARNING:tensorflow:From /Users/akash/anaconda3/lib/python3.7/site-packages/tensorflow/python/ops/math_ops.py:3066: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.
    Instructions for updating:
    Use tf.cast instead.





```python
if not os.path.exists('out/'):
    os.makedirs('out/')

i = 0

for it in range(1000000):
    X_mb, _ = mnist.train.next_batch(mb_size)

    _, loss = sess.run([solver, vae_loss], feed_dict={X: X_mb})

    if it % 1000 == 0:
        print('Iter: {}'.format(it))
        print('Loss: {:.4}'. format(loss))
        print()

        samples = sess.run(X_samples, feed_dict={z: np.random.randn(16, z_dim)})
        print (type(samples))
        #fig = plot(samples)
        #samples('out/{}.png'.format(str(i).zfill(3)), bbox_inches='tight')
        i += 1
        #plt.close(fig)
```



# References : 
    1.https://wiseodd.github.io/techblog/2016/12/10/variational-autoencoder/
    2.https://github.com/wiseodd/generative-models/blob/master/RBM/rbm_binary_cd.py
    3.https://towardsdatascience.com/intuitively-understanding-variational-autoencoders-1bfe67eb5daf
    4.http://anotherdatum.com/vae2.html


```python

```

e649070de19651ba88ec1d9a022fab89db5816f7
ipynb
Jupyter Notebook
Assignment 5.ipynb
Mgosi/Pattern-Recognition
e4a51b41e3ac0e64456adb629da2e8d8825c6b12
MIT
Assignment 5.ipynb
Mgosi/Pattern-Recognition
e4a51b41e3ac0e64456adb629da2e8d8825c6b12
MIT
Assignment 5.ipynb
Mgosi/Pattern-Recognition
e4a51b41e3ac0e64456adb629da2e8d8825c6b12
MIT
Qwen/Qwen-72B
1. YES
2. YES