ICA natural images und Binary RBM added

MelJan · Apr 22, 2017 · 1952ac7 · 1952ac7
1 parent d870209
commit 1952ac7
Show file tree

Hide file tree

Showing 8 changed files with 230 additions and 44 deletions.
diff --git a/docs/tutorial.rst b/docs/tutorial.rst
@@ -3,9 +3,13 @@ Tutorials
 .. toctree::
    :maxdepth: 2
 
-   Eigenfaces<tutorials/eigenfaces.rst>
-
    PCA 2D<tutorials/PCA_2D_example.rst>
 
+   Eigenfaces<tutorials/eigenfaces.rst>
+
    ICA 2D<tutorials/ICA_2D_example.rst>
 
+   ICA on natural images<tutorials/ICA_natural_images.rst>
+
+   ICA on natural images<tutorials/Binary_RBM_MNIST_small.rst>
+
diff --git a/docs/tutorials/Binary_RBM_MNIST_small.rst b/docs/tutorials/Binary_RBM_MNIST_small.rst
@@ -0,0 +1,130 @@
+Binary RBM
+==========================================================
+
+Example for training a centered and normal Binary Restricted Boltzmann machine on the MNIST handwritten digit dataset and the flipped version.
+The model is small enough to calculate the exact Log-Likelihood and annealed importance sampling and reverse annealed importance sampling are evaluated.
+
+For an analysis of advantage of centering in RBMs see `How to Center Binary Restricted Boltzmann Machines (Vol. 1311.1354). Melchior, J., Fischer, A., Wang, N., & Wiskott, L.. (2013). arXiv.org e-Print archive. <http://arxiv.org/pdf/1311.1354.pdf>`_
+
+Learned filters of a centered binary RBM on the MNIST dataset.
+
+.. figure:: images/BRBM_small_centered_weights.png
+   :scale: 75 %
+   :alt: 100 gray scale natural image patch examples
+
+Sampling results for some examples. The first row shows training data and the following rows are the results after one Gibbs-sampling step starting from the previous row.
+
+.. figure:: images/BRBM_small_centered_samples.png
+   :scale: 75 %
+   :alt: 100 gray scale natural image patch examples
+
+The Log-Likelihood is calculated using the exact Partition function, an annealed importance sampling estimation (optimistic) and reverse annealed importance sampling estimation (pessimistic).
+
+.. code-block:: Python
+
+   Training
+   Epoch    Recon. Error   Log likelihood  Expected End-Time
+   10       0.0519         -152.8375       2017-04-22 19:19:20.492937
+   20       0.0512         -144.6574       2017-04-22 19:19:18.518975
+   30       0.0510         -143.9597       2017-04-22 19:19:18.118906
+   End-time:                               2017-04-22 19:19:06.424132
+   Training time:          0:03:37.913801
+   True Partition:         332.061422631 (LL: -143.295691625)
+   AIS Partition:          331.856312898 (LL: -143.090581892)
+   reverse AIS Partition:  333.975682243 (LL: -145.209951237)
+
+
+The code can also be executed without centering by setting
+
+.. code-block:: python
+
+   update_offsets = 0.0
+
+Resulting in the following weights and sampling steps.
+
+.. figure:: images/BRBM_small_normal_weights.png
+   :scale: 75 %
+   :alt: 100 gray scale natural image patch examples
+
+Sampling results for some examples. The first row shows training data and the following rows are the result after one Gibbs-sampling step.
+
+.. figure:: images/BRBM_small_normal_samples.png
+   :scale: 75 %
+   :alt: 100 gray scale natural image patch examples
+
+The Log-Likelihood for this model is significantly worse (8 nats lower).
+
+.. code-block:: Python
+
+   Epoch    Recon. Error   Log likelihood  Expected End-Time
+   10       0.0550         -153.1939       2017-04-22 19:48:31.196758
+   20       0.0538         -158.0558       2017-04-22 19:48:30.991562
+   30       0.0534         -154.9978       2017-04-22 19:48:31.701043
+   End-time:                               2017-04-22 19:48:17.356501
+   Training time:          0:03:51.361354
+   True Partition:         209.93080041  (LL: -151.678071606)
+   AIS Partition:          209.922151068 (LL: -151.669422263)
+   reverse AIS Partition:  209.55108856  (LL: -151.298359755)
+
+Further, the models can be trained on the flipped version of MNIST (1-MNIST).
+
+.. code-block:: python
+
+   flipped = True
+
+.. _code:
+
+While the centerer model has a similar performance on the flipped version,
+
+.. code-block:: Python
+
+   Epoch    Recon. Error  Log likelihood  Expected End-Time
+   10       0.0520        -145.0683       2017-04-22 20:10:32.506455
+   20       0.0513        -143.1599       2017-04-22 20:10:41.450600
+   30       0.0510        -143.4857       2017-04-22 20:10:40.404681
+   End-time:                              2017-04-22 20:10:29.587139
+   Training time:         0:08:21.800137
+   True Partition:        330.041980724 (LL: -142.715842237)
+   AIS Partition:         329.78134331  (LL: -142.455204822)
+   reverse AIS Partition: 331.047724145 (LL: -143.721585658)
+
+.. figure:: images/BRBM_small_centered_weights_flipped.png
+   :scale: 75 %
+   :alt: 100 gray scale natural image patch examples
+
+Sampling results for some examples. The first row shows training data and the following rows are the result after one Gibbs-sampling step.
+
+.. figure:: images/BRBM_small_centered_samples_flipped.png
+   :scale: 75 %
+   :alt: 100 gray scale natural image patch examples
+
+the normal RBM fails to learn the distribution.
+
+.. code-block:: Python
+
+   Epoch    Recon. Error  Log likelihood  Expected End-Time
+   10       0.0588        -183.3271       2017-04-22 20:10:26.751454
+   20       0.0580        -206.4093       2017-04-22 20:10:35.824502
+   30       0.0576        -180.9389       2017-04-22 20:10:33.400136
+   End-time:                              2017-04-22 20:10:12.549085
+   Training time:         0:07:56.808453
+   True Partition:        3755.27224976 (LL: -193.127765428)
+   AIS Partition:         3755.25179256 (LL: -193.10730823)
+   reverse AIS Partition: 3755.18060949 (LL: -193.036125159)
+
+.. figure:: images/BRBM_small_normal_weights_flipped.png
+   :scale: 75 %
+   :alt: 100 gray scale natural image patch examples
+
+.. figure:: images/BRBM_small_normal_samples_flipped.png
+   :scale: 75 %
+   :alt: 100 gray scale natural image patch examples
+
+Source code
+***********
+
+.. figure:: images/download_icon.png
+   :scale: 20 %
+   :target: https://github.com/MelJan/PyDeep/blob/master/examples/ICA_natural_images.py
+
+.. literalinclude:: ../../examples/ICA_natural_images.py
diff --git a/docs/tutorials/ICA_2D_example.rst b/docs/tutorials/ICA_2D_example.rst
@@ -45,7 +45,7 @@ We can also back-project the ICA projection matrix back to the original space an
    :scale: 75 %
    :alt: Examples of ICA 2D
 
-For a real-world application see the eigenfaces example.
+For a real-world application see the ICA_Natural_Images example.
 
 .. _code:
 

diff --git a/docs/tutorials/ICA_natural_images.rst b/docs/tutorials/ICA_natural_images.rst
@@ -0,0 +1,39 @@
+Independent Component Analysis on a natural image patches
+==========================================================
+
+Example for Independent Component Analysis (`ICA <https://en.wikipedia.org/wiki/Principal_component_analysis>`_)
+on natural image patches. The independent components (columns of the ICA projection matrix) of natural image patches are edge detector filters.
+
+See ICA_2D_example for a theoretical introduction.
+
+Visualization of 100 examples of the gray scale natural image dataset.
+
+.. figure:: images/ICA_natural_images_data.png
+   :scale: 75 %
+   :alt: 100 gray scale natural image patch examples
+
+The corresponding whitened image patches.
+
+.. figure:: images/ICA_natural_images_data_whitened.png
+   :scale: 75 %
+   :alt: 100 gray scale natural image patch examples
+
+The learned filters/independent components learned from the whitened natural image patches.
+
+.. figure:: images/ICA_natural_images_filter.png
+   :scale: 75 %
+   :alt: 100 gray scale natural image patch examples
+
+See `Gaussian-binary restricted Boltzmann machines for modeling natural image statistics. Melchior, J., Wang, N., & Wiskott, L.. (2017). PLOS ONE, 12(2), 1–24. <http://doi.org/10.1371/journal.pone.0171015>`_
+for a analysis of ICA and GRBM.
+
+.. _code:
+
+Source code
+***********
+
+.. figure:: images/download_icon.png
+   :scale: 20 %
+   :target: https://github.com/MelJan/PyDeep/blob/master/examples/ICA_natural_images.py
+
+.. literalinclude:: ../../examples/ICA_natural_images.py
diff --git a/docs/tutorials/eigenfaces.rst b/docs/tutorials/eigenfaces.rst
@@ -46,4 +46,4 @@ Source code
    :scale: 20 %
    :target: https://github.com/MelJan/PyDeep/blob/master/examples/eigenfaces.py
 
-.. literalinclude:: ../../examples/eigenfaces.py
+.. literalinclude:: ../../examples/Eigenfaces.py
diff --git a/examples/small_binary_RBM_MNIST.py → examples/Binary_RBM_MNIST_small.py b/examples/small_binary_RBM_MNIST.py → examples/Binary_RBM_MNIST_small.py
@@ -44,21 +44,45 @@
 # Set random seed (optional)
 numx.random.seed(42)
 
+# normal RBM
+#update_offsets = 0.0
+# centered RBM
+update_offsets = 0.01
+
+# Flipped/Inverse MNIST
+#flipped = True
+# Normal MNIST
+flipped = False
+
 # Input and hidden dimensionality
 v1 = v2 = 28
 h1 = h2 = 4
 
 # Load data , get it from 'deeplearning.net/data/mnist/mnist.pkl.gz'
 train_data = io.load_mnist("../../data/mnist.pkl.gz", True)[0]
 
+# Flip the dataset if chosen
+if flipped:
+    train_data = 1-train_data
+
 # Training paramters
 batch_size = 100
 epochs = 39
 
-# Create trainer and model
-rbm = model.BinaryBinaryRBM(number_visibles=v1 * v2,
-                            number_hiddens=h1 * h2,
-                            data=train_data)
+# Create centered or normal model
+if update_offsets <= 0.0:
+    rbm = model.BinaryBinaryRBM(number_visibles=v1 * v2,
+                                number_hiddens=h1 * h2,
+                                data=train_data,
+                                initial_visible_offsets=0.0,
+                                initial_hidden_offsets=0.0)
+else:
+    rbm = model.BinaryBinaryRBM(number_visibles=v1 * v2,
+                                number_hiddens=h1 * h2,
+                                data=train_data,
+                                initial_visible_offsets='AUTO',
+                                initial_hidden_offsets='AUTO')
+# Create trainer
 trainer = trainer.PCD(rbm, batch_size)
 
 # Measuring time
@@ -75,14 +99,18 @@
     # Loop over all batches
     for b in range(0, train_data.shape[0], batch_size):
         batch = train_data[b:b + batch_size, :]
-        trainer.train(data=batch, epsilon=0.05)
+        trainer.train(data=batch,
+                      epsilon=0.05,
+                      update_visible_offsets=update_offsets,
+                      update_hidden_offsets=update_offsets)
 
     # Calculate Log-Likelihood, reconstruction error and expected end time every 10th epoch
     if epoch % 10 == 0:
         logZ = estimator.partition_function_factorize_h(rbm)
         ll = numx.mean(estimator.log_likelihood_v(rbm, logZ, train_data))
         re = numx.mean(estimator.reconstruction_error(rbm, train_data))
-        print('{}\t\t{:.4f}\t\t\t{:.4f}\t\t\t{}'.format(epoch, re, ll, measurer.get_expected_end_time(epoch, epochs)))
+        print('{}\t\t{:.4f}\t\t\t{:.4f}\t\t\t{}'.format(
+            epoch, re, ll, measurer.get_expected_end_time(epoch, epochs)))
     else:
         print(epoch)
 

diff --git a/examples/eigenfaces.py → examples/Eigenfaces.py b/examples/eigenfaces.py → examples/Eigenfaces.py
diff --git a/examples/ica_on_natural_images.py → examples/ICA_natural_images.py b/examples/ica_on_natural_images.py → examples/ICA_natural_images.py
@@ -1,10 +1,10 @@
-""" Example for the Independent component analysis on natural image patches.
+""" Example for the Independent Component Analysis (ICA) on natural image patches.
 
     :Version:
         1.1.0
 
     :Date:
-        08.04.2017
+        22.04.2017
 
     :Author:
         Jan Melchior
@@ -32,27 +32,30 @@
         along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 """
-# Import PCA, numpy, input output functions, and visualization functions
+
+# Import ZCA, ICA, numpy, input output functions, and visualization functions
 import numpy as numx
 from pydeep.preprocessing import ICA, ZCA
 import pydeep.misc.io as io
 import pydeep.misc.visualization as vis
 
-# Set the random seed (optional, if stochastic processes are involved we always get the same results)
+# Set the random seed
+# (optional, if stochastic processes are involved we always get the same results)
 numx.random.seed(42)
 
-# Load the data
+# Load the data (automatic download if it does not exist at given path)
 data = io.load_natural_image_patches('../../data/NaturalImage.mat')
 
 # Specify image width and height for displaying
 width = height = 14
 
-# Create a ZCA node to whiten the data and train it (you could also use PCA whitened=True)
+# Create a ZCA node to whiten the data and train it
+# (you could also use PCA whitened=True)
 zca = ZCA(input_dim=width * height)
 zca.train(data=data)
 
-# ZCA projects the whitened data back to the original space, thus does not perform a
-# dimensionality reduction but a whitening in the original space
+# ZCA projects the whitened data back to the original space, thus does not
+# perform a dimensionality reduction but a whitening in the original space
 whitened_data = zca.project(data)
 
 # Create a ZCA node and train it (you could also use PCA whitened=True)
@@ -62,45 +65,27 @@
           convergence=1.0,
           status=True)
 
-# Show eigenvectors of the covariance matrix
-eigenvectors = vis.tile_matrix_rows(matrix=zca.projection_matrix,
-                                    tile_width=width,
-                                    tile_height=height,
-                                    num_tiles_x=width,
-                                    num_tiles_y=height,
-                                    border_size=1,
-                                    normalized=True)
-vis.imshow_matrix(matrix=eigenvectors,
-                  windowtitle='Eigenvectors of the covariance matrix')
-
 # Show whitened images
-images = vis.tile_matrix_rows(matrix=data[0:width*height].T,
+images = vis.tile_matrix_rows(matrix=data[0:100].T,
                               tile_width=width,
                               tile_height=height,
-                              num_tiles_x=width,
-                              num_tiles_y=height,
+                              num_tiles_x=10,
+                              num_tiles_y=10,
                               border_size=1,
                               normalized=True)
 vis.imshow_matrix(matrix=images,
-                  windowtitle='Some image patches')
+                  windowtitle='First 100 image patches')
 
 # Show some whitened images
-images = vis.tile_matrix_rows(matrix=whitened_data[0:width*height].T,
+images = vis.tile_matrix_rows(matrix=whitened_data[0:100].T,
                               tile_width=width,
                               tile_height=height,
-                              num_tiles_x=width,
-                              num_tiles_y=height,
+                              num_tiles_x=10,
+                              num_tiles_y=10,
                               border_size=1,
                               normalized=True)
 vis.imshow_matrix(matrix=images,
-                  windowtitle='Some whitened image patches')
-
-# Plot the cumulative sum of teh Eigenvalues.
-eigenvalue_sum = numx.cumsum(zca.eigen_values)
-vis.imshow_plot(matrix=eigenvalue_sum,
-                windowtitle="Cumulative sum of Eigenvalues")
-vis.xlabel("Eigenvalue index")
-vis.ylabel("Sum of Eigenvalues 0 to index")
+                  windowtitle='First 100 image patches whitened')
 
 # Show the ICA filters/bases
 ica_filters = vis.tile_matrix_rows(matrix=ica.projection_matrix,