diff --git a/bob/pad/face/extractor/VideoHistOfSparseCodes.py b/bob/pad/face/extractor/VideoHistOfSparseCodes.py
index 61747358cf625ba3043c521c46bcd78d800f0d3d..5ef9c47abfa02e42165562e8b5cafb10005b7c01 100644
--- a/bob/pad/face/extractor/VideoHistOfSparseCodes.py
+++ b/bob/pad/face/extractor/VideoHistOfSparseCodes.py
@@ -102,6 +102,21 @@ class VideoHistOfSparseCodes(Extractor, object):
         return frame_container
 
 
+    #==========================================================================
+    def reduce_features_number(self, list_of_arrays):
+        """
+        Reduce the number of features.
+        """
+
+        return_list = []
+
+        for item in list_of_arrays:
+
+            return_list.append( item[1][32:] )
+
+        return return_list
+
+
     #==========================================================================
     def __call__(self, frames):
         """
@@ -119,7 +134,9 @@ class VideoHistOfSparseCodes(Extractor, object):
             Histograms of sparse codes stored in the FrameContainer.
         """
 
-        histograms = self.comp_hist_of_sparse_codes(frames, self.method)
+#        histograms = self.comp_hist_of_sparse_codes(frames, self.method)
+
+        histograms = self.reduce_features_number(frames)
 
         frame_container = self.convert_sparse_codes_to_frame_container(histograms)
 
diff --git a/bob/pad/face/preprocessor/VideoSparseCoding.py b/bob/pad/face/preprocessor/VideoSparseCoding.py
index 944ca82e30ae792dd269d8ae927b61bbbaa01ce5..016ada30f992c02dd035079a45d394209a4eface 100644
--- a/bob/pad/face/preprocessor/VideoSparseCoding.py
+++ b/bob/pad/face/preprocessor/VideoSparseCoding.py
@@ -77,6 +77,11 @@ class VideoSparseCoding(Preprocessor, object):
         "mean" and "hist". This argument is valid only if ``extract_histograms_flag``
         is set to ``True``.
         Default: "hist".
+
+    ``comp_reconstruct_err_flag`` : :py:class:`bool`
+        If this flag is set to ``True`` resulting feature vector will be a
+        reconstruction error, not a histogram.
+        Default: ``False``.
     """
 
     #==========================================================================
@@ -89,6 +94,7 @@ class VideoSparseCoding(Preprocessor, object):
                  frame_step = 1,
                  extract_histograms_flag = False,
                  method = "hist",
+                 comp_reconstruct_err_flag = False,
                  **kwargs):
 
         super(VideoSparseCoding, self).__init__(block_size = block_size,
@@ -98,6 +104,7 @@ class VideoSparseCoding(Preprocessor, object):
                                                 dictionary_file_names = dictionary_file_names,
                                                 frame_step = frame_step,
                                                 extract_histograms_flag = extract_histograms_flag,
+                                                comp_reconstruct_err_flag = comp_reconstruct_err_flag,
                                                 method = method)
 
         self.block_size = block_size
@@ -108,6 +115,7 @@ class VideoSparseCoding(Preprocessor, object):
         self.frame_step = frame_step
         self.extract_histograms_flag = extract_histograms_flag
         self.method = method
+        self.comp_reconstruct_err_flag = comp_reconstruct_err_flag
 
         self.video_preprocessor = bob.bio.video.preprocessor.Wrapper()
 
@@ -799,6 +807,164 @@ class VideoSparseCoding(Preprocessor, object):
         return frame_container
 
 
+    #==========================================================================
+    def compute_patches_mean_squared_errors(self, sparse_codes, original_data, dictionary):
+        """
+        This function computes mean squared errors (MSE) for each feature (column)
+        in the reconstructed array of vectorized patches. The patches are reconstructed
+        given array of sparse codes and a dictionary.
+
+        **Parameters:**
+
+        ``sparse_codes`` : 2D :py:class:`numpy.ndarray`
+            An array of sparse codes. Each row contains a sparse code encoding a
+            vectorized patch. The dimensionality of the array:
+            (``n_samples`` x ``n_words_in_dictionary``).
+
+        ``original_data`` : 2D :py:class:`numpy.ndarray`
+            An array with original vectorized patches.
+            The dimensionality of the array:
+            (``n_samples`` x ``n_features_in_patch``).
+
+        ``dictionary`` : 2D :py:class:`numpy.ndarray`
+            A dictionary with vectorized visual words.
+            The dimensionality of the array:
+            (``n_words_in_dictionary`` x ``n_features_in_patch``).
+
+        **Returns:**
+
+        ``squared_errors`` : 1D :py:class:`numpy.ndarray`
+            MSE for each feature across all patches/samples.
+            The dimensionality of the array:
+            (``n_features_in_patch``, ).
+        """
+
+        recovered_data = np.dot(sparse_codes, dictionary)
+
+        n_samples = recovered_data.shape[0]
+
+        squared_error = 1.0 / n_samples * np.sum((recovered_data - original_data) ** 2, axis=0)
+
+        return squared_error
+
+
+    #==========================================================================
+    def compute_mse_for_all_patches_types(self, sparse_codes_list, original_data_list, dictionary_list):
+        """
+        This function computes mean squared errors (MSE) for all types of patches:
+        frontal, horizontal, and vertical. In this case the function
+        ``compute_patches_mean_squared_errors`` is called in a loop for all
+        values in the input lists.
+
+        **Parameters:**
+
+        ``sparse_codes_list`` : [2D :py:class:`numpy.ndarray`]
+            A list with arrays of sparse codes. Each row in the arrays contains a
+            sparse code encoding a vectorized patch of particular type.
+            The dimensionality of the each array:
+            (``n_samples`` x ``n_words_in_dictionary``).
+
+        ``original_data_list`` : [2D :py:class:`numpy.ndarray`]
+            A list of arrays with original vectorized patches of various types.
+            The dimensionality of the arrays might be different for various types
+            of the patches:
+            (``n_samples`` x ``n_features_in_patch_of_particular_type``).
+
+        ``dictionary_list`` : [2D :py:class:`numpy.ndarray`]
+            A list of dictionaries with vectorized visual words of various types.
+            The dimensionality of the arrays might be different for various types
+            of the patches:
+            (``n_words_in_dictionary`` x ``n_features_in_patch_of_particular_type``).
+
+        **Returns:**
+
+        ``squared_errors`` : 2D :py:class:`numpy.ndarray`
+            First row:
+            MSE of features for various types of patches concatenated into a single
+            vector.
+            Second row:
+            The same as above but MSE are sorted for each type of patches.
+            The dimensionality of the array:
+            (2 x ``n_features_in_patch_of_all_types``).
+        """
+
+        squared_errors = []
+
+        squared_errors_sorted = []
+
+        for sparse_codes, original_data, dictionary in zip(sparse_codes_list, original_data_list, dictionary_list):
+
+            squared_error = self.compute_patches_mean_squared_errors(sparse_codes, original_data, dictionary)
+
+            squared_error_sorted = np.sort(squared_error)
+
+            squared_errors.append(squared_error)
+
+            squared_errors_sorted.append(squared_error_sorted)
+
+        squared_errors = np.hstack(squared_errors)
+
+        squared_errors_sorted = np.hstack(squared_errors_sorted)
+
+        squared_errors = np.vstack([squared_errors, squared_errors_sorted])
+
+        return squared_errors
+
+
+    #==========================================================================
+    def compute_mse_for_all_stacks(self, video_codes_list, patches_list, dictionary_list):
+        """
+        Call ``compute_mse_for_all_patches_types`` for data coming from all stacks
+        of facial images.
+
+        **Parameters:**
+
+        ``video_codes_list`` : [ [2D :py:class:`numpy.ndarray`] ]
+            A list with ``frontal_video_codes``, ``horizontal_video_codes``, and
+            ``vertical_video_codes`` as returned by ``get_sparse_codes_for_list_of_patches``
+            method of this class.
+
+        ``patches_list`` : [ [2D :py:class:`numpy.ndarray`] ]
+            A list with ``frontal_patches``, ``horizontal_patches``, and
+            ``vertical_patches`` as returned by ``extract_patches_from_blocks``
+            method of this class.
+
+        ``dictionary_list`` : [2D :py:class:`numpy.ndarray`]
+            A list of dictionaries with vectorized visual words of various types.
+            The dimensionality of the arrays might be different for various types
+            of the patches:
+            (``n_words_in_dictionary`` x ``n_features_in_patch_of_particular_type``).
+
+        **Returns:**
+
+        ``squared_errors_list`` : [2D :py:class:`numpy.ndarray`]
+            A list of ``squared_errors`` as returned by ``compute_mse_for_all_patches_types``
+            method of this class.
+        """
+
+        fcs = video_codes_list[0]
+        hcs = video_codes_list[1]
+        vcs = video_codes_list[2]
+
+        fps = patches_list[0]
+        hps = patches_list[1]
+        vps = patches_list[2]
+
+        squared_errors_list = []
+
+        for fc, hc, vc, fp, hp, vp in zip(fcs, hcs, vcs, fps, hps, vps):
+
+            sparse_codes_list = [fc, hc, vc]
+
+            original_data_list = [fp, hp, vp]
+
+            squared_errors = self.compute_mse_for_all_patches_types(sparse_codes_list, original_data_list, dictionary_list)
+
+            squared_errors_list.append(squared_errors)
+
+        return squared_errors_list
+
+
     #==========================================================================
     def __call__(self, frames, annotations):
         """
@@ -859,18 +1025,37 @@ class VideoSparseCoding(Preprocessor, object):
         # Download the dictionaries:
         dictionary_frontal, dictionary_horizontal, dictionary_vertical = self.load_the_dictionaries(self.dictionary_file_names)
 
+        # Select subset of patches if ``frame_step`` > 1:
+        frontal_patches_subset = frontal_patches[::self.frame_step]
+        horizontal_patches_subset = horizontal_patches[::self.frame_step]
+        vertical_patches_subset = vertical_patches[::self.frame_step]
+
         # Compute sparse codes for all patches of all types:
-        frontal_video_codes = self.get_sparse_codes_for_list_of_patches(frontal_patches[::self.frame_step], dictionary_frontal)
-        horizontal_video_codes = self.get_sparse_codes_for_list_of_patches(horizontal_patches[::self.frame_step], dictionary_horizontal)
-        vertical_video_codes = self.get_sparse_codes_for_list_of_patches(vertical_patches[::self.frame_step], dictionary_vertical)
+        frontal_video_codes = self.get_sparse_codes_for_list_of_patches(frontal_patches_subset, dictionary_frontal)
+        horizontal_video_codes = self.get_sparse_codes_for_list_of_patches(horizontal_patches_subset, dictionary_horizontal)
+        vertical_video_codes = self.get_sparse_codes_for_list_of_patches(vertical_patches_subset, dictionary_vertical)
+
+        if self.comp_reconstruct_err_flag:
+
+            video_codes_list = [frontal_video_codes, horizontal_video_codes, vertical_video_codes]
+
+            patches_list = [frontal_patches_subset, horizontal_patches_subset, vertical_patches_subset]
+
+            dictionary_list = [dictionary_frontal, dictionary_horizontal, dictionary_vertical]
+
+            squared_errors_list = self.compute_mse_for_all_stacks(video_codes_list, patches_list, dictionary_list)
+
+            frame_container = self.convert_arrays_to_frame_container(squared_errors_list)
+
+        else:
 
-        frame_container = self.convert_sparse_codes_to_frame_container([frontal_video_codes, horizontal_video_codes, vertical_video_codes])
+            frame_container = self.convert_sparse_codes_to_frame_container([frontal_video_codes, horizontal_video_codes, vertical_video_codes])
 
-        if self.extract_histograms_flag: # in this case histograms will be extracted in the preprocessor , no feature extraction is needed then
+            if self.extract_histograms_flag: # in this case histograms will be extracted in the preprocessor , no feature extraction is needed then
 
-            histograms = self.comp_hist_of_sparse_codes(frame_container, self.method)
+                histograms = self.comp_hist_of_sparse_codes(frame_container, self.method)
 
-            frame_container = self.convert_arrays_to_frame_container(histograms)
+                frame_container = self.convert_arrays_to_frame_container(histograms)
 
         return frame_container