Repository: dirty-cat/dirty_cat
Pull request #153 (2020-12-03 12:02:36)

Adding Gamma poisson factorization
Changes:

- Added `gamma_poisson_factorization.py` which implements online Gamma-Poisson factorization for encoding string variables.
- Added the corresponding tests in`test_gamma_poisson_factorization`.
- Modified examples 02 and 03 to include this method.
- Updated `CHANGES.rst` and `index.rst` to describe this new method.


Patch:
diff --git a/.circleci/config.yml b/.circleci/config.yml
index 1e3e93c..168ec79 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -7,7 +7,7 @@ jobs:
     environment:
       - MINICONDA_PATH: ~/miniconda
       - CONDA_ENV_NAME: testenv
-      - PYTHON_VERSION: 3
+      - PYTHON_VERSION: 3.8
     steps:
       - checkout
       - run: ./build_tools/circle/checkout_merge_commit.sh
diff --git a/.gitignore b/.gitignore
index 9685ae1..e3e843f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -46,3 +46,6 @@ coverage.xml
 
 # Notebooks
 .ipynb_checkpoints/
+
+# VS Code
+.vscode
diff --git a/CHANGES.rst b/CHANGES.rst
index 39197a7..cf5528f 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,6 +1,11 @@
 Release 0.0.7
 =============
 
+* **GapEncoder**: Added online Gamma-Poisson factorization through the
+  :class:`GapEncoder` class. This method discovers latent categories formed
+  via combinations of substrings, and encodes string data as combinations of
+  these categories. To be used if interpretability is important.
+
 * **MinHashEncoder**: Added ``minhash_encoder.py`` and ``fast_hast.py`` files
   that implement minhash encoding through the ``MinHashEncoder`` class.
   This method allows for fast and scalable encoding of string categorical
diff --git a/dirty_cat/VERSION.txt b/dirty_cat/VERSION.txt
index bbdeab6..5a5831a 100644
--- a/dirty_cat/VERSION.txt
+++ b/dirty_cat/VERSION.txt
@@ -1,1 +1,1 @@
-0.0.5
+0.0.7
diff --git a/dirty_cat/__init__.py b/dirty_cat/__init__.py
index fef06a1..bb9c9e3 100644
--- a/dirty_cat/__init__.py
+++ b/dirty_cat/__init__.py
@@ -10,5 +10,7 @@ with open(version_file) as fh:
 from .similarity_encoder import SimilarityEncoder
 from .target_encoder import TargetEncoder
 from .minhash_encoder import MinHashEncoder
+from .gap_encoder import GapEncoder
 
-__all__ = ['SimilarityEncoder', 'TargetEncoder', 'MinHashEncoder']
+__all__ = ['SimilarityEncoder', 'TargetEncoder', 'MinHashEncoder',
+           'GapEncoder']
diff --git a/dirty_cat/gap_encoder.py b/dirty_cat/gap_encoder.py
new file mode 100644
index 0000000..56cadea
--- /dev/null
+++ b/dirty_cat/gap_encoder.py
@@ -0,0 +1,639 @@
+"""
+Online Gamma-Poisson factorization of string arrays.
+The principle is as follows:
+    1. Given an input string array X, we build its bag-of-n-grams
+       representation V (n_samples, vocab_size).
+    2. Instead of using the n-grams counts as encodings, we look for low-
+       dimensional representations by modeling n-grams counts as linear
+       combinations of topics V = HW, with W (n_topics, vocab_size) the topics
+       and H (n_samples, n_topics) the associated activations.
+    3. Assuming that n-grams counts follow a Poisson law, we fit H and W to
+       maximize the likelihood of the data, with a Gamma prior for the
+       activations H to induce sparsity.
+    4. In practice, this is equivalent to a non-negative matrix factorization
+       with the Kullback-Leibler divergence as loss, and a Gamma prior on H.
+       We thus optimize H and W with the multiplicative update method.
+"""
+import numpy as np
+from distutils.version import LooseVersion
+from scipy import sparse
+from sklearn import __version__ as sklearn_version
+from sklearn.utils import check_random_state, gen_batches
+from sklearn.utils.extmath import row_norms, safe_sparse_dot
+from sklearn.base import BaseEstimator, TransformerMixin
+from sklearn.feature_extraction.text import CountVectorizer, HashingVectorizer
+from sklearn.cluster import KMeans
+from sklearn.neighbors import NearestNeighbors
+
+if LooseVersion(sklearn_version) < LooseVersion('0.22'):
+    from sklearn.cluster.k_means_ import _k_init
+elif LooseVersion(sklearn_version) < LooseVersion('0.24'):
+    from sklearn.cluster._kmeans import _k_init
+else:
+    from sklearn.cluster import kmeans_plusplus
+
+if LooseVersion(sklearn_version) < LooseVersion('0.22'):
+    from sklearn.decomposition.nmf import _beta_divergence
+else:
+    from sklearn.decomposition._nmf import _beta_divergence
+
+
+class GapEncoder(BaseEstimator, TransformerMixin):
+    """
+    This encoder can be understood as a continuous encoding on a set of latent
+    categories estimated from the data. The latent categories are built by
+    capturing combinations of substrings that frequently co-occur.
+    
+    The GapEncoder supports online learning on batches of data for
+    scalability through the partial_fit method.
+
+    Parameters
+    ----------
+
+    n_components : int, default=10
+        Number of latent categories used to model string data.
+
+    batch_size : int, default=128
+        Number of samples per batch.
+
+    gamma_shape_prior : float, default=1.1
+        Shape parameter for the Gamma prior distribution.
+
+    gamma_scale_prior : float, default=1.0
+        Scale parameter for the Gamma prior distribution.
+
+    rho : float, default=0.95
+        Weight parameter for the update of the W matrix.
+    
+    rescale_rho : bool, default=False
+        If true, use rho ** (batch_size / len(X)) instead of rho to obtain an
+        update rate per iteration that is independent of the batch size.
+
+    hashing : bool, default=False
+        If true, HashingVectorizer is used instead of CountVectorizer.
+        It has the advantage of being very low memory scalable to large
+        datasets as there is no need to store a vocabulary dictionary in
+        memory.
+
+    hashing_n_features : int, default=2**12
+        Number of features for the HashingVectorizer. Only relevant if
+        hashing=True.
+
+    init : str, default='k-means++'
+        Initialization method of the W matrix.
+        Options: {'k-means++', 'random', 'k-means'}.
+        If init='k-means++', we use the init method of sklearn.cluster.KMeans.
+        If init='random', topics are initialized with a Gamma distribution.
+        If init='k-means', topics are initialized with a KMeans on the n-grams
+        counts. This usually makes convergence faster but is a bit slower.
+
+    tol : float, default=1e-4
+        Tolerance for the convergence of the matrix W.
+
+    min_iter : int, default=2
+        Minimum number of iterations on the input data.
+
+    max_iter : int, default=5
+        Maximum number of iterations on the input data.
+
+    ngram_range : tuple, default=(2, 4)
+        The range of ngram length that will be used to build the
+        bag-of-n-grams representation of the input data.
+
+    analyzer : str, default='char'.
+        Analyzer parameter for the CountVectorizer/HashingVectorizer.
+        Options: {‘word’, ‘char’, ‘char_wb’}, describing whether the matrix V
+        to factorize should be made of word counts or character n-gram counts.
+        Option ‘char_wb’ creates character n-grams only from text inside word
+        boundaries; n-grams at the edges of words are padded with space.
+
+    add_words : bool, default=False
+        If true, add the words counts to the bag-of-n-grams representation
+        of the input data.
+
+    random_state : int or None, default=None
+        Pass an int for reproducible output across multiple function calls.
+
+    rescale_W : bool, default=True
+        If true, the weight matrix W is rescaled at each iteration
+        to have an l1 norm equal to 1 for each row.
+    
+    max_iter_e_step : int, default=20
+        Maximum number of iterations to adjust the activations h at each step.
+
+
+    Attributes
+    ----------
+
+    References
+    ----------
+    For a detailed description of the method, see
+    `Encoding high-cardinality string categorical variables
+    <https://hal.inria.fr/hal-02171256v4>`_ by Cerda, Varoquaux (2019).
+    """
+
+    def __init__(self, n_components=10, batch_size=128, gamma_shape_prior=1.1,
+                 gamma_scale_prior=1.0, rho=.95, rescale_rho=False,
+                 hashing=False, hashing_n_features=2**12, init='k-means++',
+                 tol=1e-4, min_iter=2, max_iter=5, ngram_range=(2, 4),
+                 analyzer='char', add_words=False, random_state=None,
+                 rescale_W=True, max_iter_e_step=20):
+
+        self.ngram_range = ngram_range
+        self.n_components = n_components
+        self.gamma_shape_prior = gamma_shape_prior  # 'a' parameter
+        self.gamma_scale_prior = gamma_scale_prior  # 'b' parameter
+        self.rho = rho
+        self.rho_ = self.rho
+        self.rescale_rho = rescale_rho
+        self.batch_size = batch_size
+        self.tol = tol
+        self.hashing = hashing
+        self.hashing_n_features = hashing_n_features
+        self.max_iter = max_iter
+        self.min_iter = min_iter
+        self.init = init
+        self.analyzer = analyzer
+        self.add_words = add_words
+        self.random_state = check_random_state(random_state)
+        self.rescale_W = rescale_W
+        self.max_iter_e_step = max_iter_e_step
+
+    def _init_vars(self, X):
+        """
+        Build the bag-of-n-grams representation V of X and initialize
+        the topics W.
+        """
+        # Init n-grams counts vectorizer
+        if self.hashing:
+            self.ngrams_count_ = HashingVectorizer(
+                 analyzer=self.analyzer, ngram_range=self.ngram_range,
+                 n_features=self.hashing_n_features,
+                 norm=None, alternate_sign=False)
+            if self.add_words: # Init a word counts vectorizer if needed
+                self.word_count_ = HashingVectorizer(
+                     analyzer='word',
+                     n_features=self.hashing_n_features,
+                     norm=None, alternate_sign=False)
+        else:
+            self.ngrams_count_ = CountVectorizer(
+                 analyzer=self.analyzer, ngram_range=self.ngram_range,
+                 dtype=np.float64)
+            if self.add_words:
+                self.word_count_ = CountVectorizer(dtype=np.float64)
+
+        # Init H_dict_ with empty dict to train from scratch
+        self.H_dict_ = dict()
+        # Build the n-grams counts matrix unq_V on unique elements of X
+        unq_X, lookup = np.unique(X, return_inverse=True)
+        unq_V = self.ngrams_count_.fit_transform(unq_X)
+        if self.add_words: # Add word counts to unq_V
+            unq_V2 = self.word_count_.fit_transform(unq_X)
+            unq_V = sparse.hstack((unq_V, unq_V2), format='csr')
+
+        if not self.hashing: # Build n-grams/word vocabulary
+            self.vocabulary = self.ngrams_count_.get_feature_names()
+            if self.add_words:
+                self.vocabulary = np.concatenate(
+                    (self.vocabulary, self.word_count_.get_feature_names()))
+
+        _, self.n_vocab = unq_V.shape
+        # Init the topics W given the n-grams counts V
+        self.W_, self.A_, self.B_ = self._init_w(unq_V[lookup], X)
+        # Init the activations unq_H of each unique input string
+        unq_H = _rescale_h(unq_V, np.ones((len(unq_X), self.n_components)))
+        # Update self.H_dict_ with unique input strings and their activations
+        self.H_dict_.update(zip(unq_X, unq_H))
+        if self.rescale_rho:
+            # Make update rate per iteration independant of the batch_size
+            self.rho_ = self.rho ** (self.batch_size / len(X))
+        return unq_X, unq_V, lookup
+
+    def _get_H(self, X):
+        """
+        Return the bag-of-n-grams representation of X.
+        """
+        H_out = np.empty((len(X), self.n_components))
+        for x, h_out in zip(X, H_out):
+            h_out[:] = self.H_dict_[x]
+        return H_out
+
+    def _init_w(self, V, X):
+        """
+        Initialize the topics W.
+        If self.init='k-means++', we use the init method of
+        sklearn.cluster.KMeans.
+        If self.init='random', topics are initialized with a Gamma
+        distribution.
+        If self.init='k-means', topics are initialized with a KMeans on the
+        n-grams counts.
+        """
+        if self.init == 'k-means++':
+            if LooseVersion(sklearn_version) < LooseVersion('0.24'):
+                W = _k_init(
+                    V, self.n_components,
+                    x_squared_norms=row_norms(V, squared=True),
+                    random_state=self.random_state,
+                    n_local_trials=None) + .1
+            else:
+                W, _ = kmeans_plusplus(
+                    V, self.n_components,
+                    x_squared_norms=row_norms(V, squared=True),
+                    random_state=self.random_state,
+                    n_local_trials=None)
+                W = W + .1 # To avoid restricting topics to few n-grams only
+        elif self.init == 'random':
+            W = self.random_state.gamma(
+                shape=self.gamma_shape_prior, scale=self.gamma_scale_prior,
+                size=(self.n_components, self.n_vocab))
+        elif self.init == 'k-means':
+            prototypes = get_kmeans_prototypes(
+                X, self.n_components, random_state=self.random_state)
+            W = self.ngrams_count_.transform(prototypes).A + .1
+            if self.add_words:
+                W2 = self.word_count_.transform(prototypes).A + .1
+                W = np.hstack((W, W2))
+            # if k-means doesn't find the exact number of prototypes
+            if W.shape[0] < self.n_components:
+                if LooseVersion(sklearn_version) < LooseVersion('0.24'):
+                    W2 = _k_init(
+                        V, self.n_components - W.shape[0],
+                        x_squared_norms=row_norms(V, squared=True),
+                        random_state=self.random_state,
+                        n_local_trials=None) + .1
+                else:
+                    W2, _ = kmeans_plusplus(
+                        V, self.n_components - W.shape[0],
+                        x_squared_norms=row_norms(V, squared=True),
+                        random_state=self.random_state,
+                        n_local_trials=None)
+                    W2 = W2 + .1
+                W = np.concatenate((W, W2), axis=0)
+        else:
+            raise AttributeError(
+                'Initialization method %s does not exist.' % self.init)
+        W /= W.sum(axis=1, keepdims=True)
+        A = np.ones((self.n_components, self.n_vocab)) * 1e-10
+        B = A.copy()
+        return W, A, B
+
+    def fit(self, X, y=None):
+        """
+        Fit the GapEncoder on batches of X.
+
+        Parameters
+        ----------
+        X : array-like, shape (n_samples, ) or (n_samples, 1)
+            The string data to fit the model on.
+        
+        Returns
+        -------
+        self
+        """
+        # Check input data shape
+        X = np.asarray(X)
+        assert X.ndim == 1 or (X.ndim == 2 and X.shape[1] == 1), f"ERROR:\
+        shape {X.shape} of input array is not supported."
+        if X.ndim == 2:
+            X = X[:, 0]
+        # Check if first item has str or np.str_ type
+        assert isinstance(X[0], str), "ERROR: Input data is not string."
+        # Make n-grams counts matrix unq_V
+        unq_X, unq_V, lookup = self._init_vars(X)
+        n_batch = (len(X) - 1) // self.batch_size + 1
+        del X
+        # Get activations unq_H
+        unq_H = self._get_H(unq_X)
+
+        for n_iter_ in range(self.max_iter):
+            # Loop over batches
+            for i, (unq_idx, idx) in enumerate(batch_lookup(
+              lookup, n=self.batch_size)):
+                if i == n_batch-1:
+                    W_last = self.W_.copy()
+                # Update the activations unq_H
+                unq_H[unq_idx] = _multiplicative_update_h(
+                    unq_V[unq_idx], self.W_, unq_H[unq_idx],
+                    epsilon=1e-3, max_iter=self.max_iter_e_step,
+                    rescale_W=self.rescale_W,
+                    gamma_shape_prior=self.gamma_shape_prior,
+                    gamma_scale_prior=self.gamma_scale_prior)
+                # Update the topics self.W_
+                _multiplicative_update_w(
+                    unq_V[idx], self.W_, self.A_, self.B_, unq_H[idx],
+                    self.rescale_W, self.rho_)
+
+                if i == n_batch-1:
+                    # Compute the norm of the update of W in the last batch
+                    W_change = np.linalg.norm(
+                        self.W_ - W_last) / np.linalg.norm(W_last)
+
+            if (W_change < self.tol) and (n_iter_ >= self.min_iter - 1):
+                break # Stop if the change in W is smaller than the tolerance
+
+        # Update self.H_dict_ with the learned encoded vectors (activations)
+        self.H_dict_.update(zip(unq_X, unq_H))
+        return self
+
+    def get_feature_names(self, n_labels=3):
+        
+        """
+        Returns the labels that best summarize the learned components/topics.
+        For each topic, labels with highest activations are selected.
+        
+        Parameters
+        ----------
+        
+        n_labels : int, default=3
+            The number of labels used to describe each topic.
+        
+        Returns
+        -------
+        
+        topic_labels : list of strings
+            The labels that best describe each topic.
+        
+        """
+        
+        vectorizer = CountVectorizer()
+        vectorizer.fit(list(self.H_dict_.keys()))
+        vocabulary = np.array(vectorizer.get_feature_names())
+        encoding = self.transform(np.array(vocabulary).reshape(-1))
+        encoding = abs(encoding)
+        encoding = encoding / np.sum(encoding, axis=1, keepdims=True)
+        n_components = encoding.shape[1]
+        topic_labels = []
+        for i in range(n_components):
+            x = encoding[:, i]
+            labels = vocabulary[np.argsort(-x)[:n_labels]]
+            topic_labels.append(labels)
+        topic_labels = [', '.join(label) for label in topic_labels]
+        return topic_labels
+
+    def score(self, X):
+        """
+        Returns the Kullback-Leibler divergence between the n-grams counts
+        matrix V of X, and its non-negative factorization HW.
+
+        Parameters
+        ----------
+        X : array-like (str), shape (n_samples, )
+            The data to encode.
+
+        Returns
+        -------
+        kl_divergence : float.
+            The Kullback-Leibler divergence.
+        """
+        # Build n-grams/word counts matrix
+        unq_X, lookup = np.unique(X, return_inverse=True)
+        unq_V = self.ngrams_count_.transform(unq_X)
+        if self.add_words:
+            unq_V2 = self.word_count_.transform(unq_X)
+            unq_V = sparse.hstack((unq_V, unq_V2), format='csr')
+
+        self._add_unseen_keys_to_H_dict(unq_X)
+        unq_H = self._get_H(unq_X)
+        # Given the learnt topics W, optimize the activations H to fit V = HW
+        for slice in gen_batches(n=unq_H.shape[0],
+                                 batch_size=self.batch_size):
+            unq_H[slice] = _multiplicative_update_h(
+                unq_V[slice], self.W_, unq_H[slice],
+                epsilon=1e-3, max_iter=self.max_iter_e_step,
+                rescale_W=self.rescale_W,
+                gamma_shape_prior=self.gamma_shape_prior,
+                gamma_scale_prior=self.gamma_scale_prior)
+        # Compute the KL divergence between V and HW
+        kl_divergence = _beta_divergence(
+            unq_V[lookup], unq_H[lookup], self.W_,
+            'kullback-leibler', square_root=False)
+        return kl_divergence
+
+    def partial_fit(self, X, y=None):
+        
+        """
+        Partial fit of the GapEncoder on X.
+        To be used in a online learning procedure where batches of data are
+        coming one by one.
+
+        Parameters
+        ----------
+        X : array-like, shape (n_samples, ) or (n_samples, 1)
+            The string data to fit the model on.
+        
+        Returns
+        -------
+        self
+        
+        """
+        
+        # Init H_dict_ with empty dict if it's the first call of partial_fit
+        if not hasattr(self, 'H_dict_'):
+            self.H_dict_ = dict()
+        # Check input data shape
+        X = np.asarray(X)
+        assert X.ndim == 1 or (X.ndim == 2 and X.shape[1] == 1), f"ERROR:\
+        shape {X.shape} of input array is not supported."
+        if X.ndim == 2:
+            X = X[:, 0]
+        # Check if first item has str or np.str_ type
+        assert isinstance(X[0], str), "ERROR: Input data is not string."
+        # Check if it is not the first batch
+        if hasattr(self, 'vocabulary'): # Update unq_X, unq_V with new batch
+            unq_X, lookup = np.unique(X, return_inverse=True)
+            unq_V = self.ngrams_count_.transform(unq_X)
+            if self.add_words:
+                unq_V2 = self.word_count_.transform(unq_X)
+                unq_V = sparse.hstack((unq_V, unq_V2), format='csr')
+
+            unseen_X = np.setdiff1d(unq_X, np.array([*self.H_dict_]))
+            unseen_V = self.ngrams_count_.transform(unseen_X)
+            if self.add_words:
+                unseen_V2 = self.word_count_.transform(unseen_X)
+                unseen_V = sparse.hstack((unseen_V, unseen_V2), format='csr')
+
+            if unseen_V.shape[0] != 0:
+                unseen_H = _rescale_h(
+                    unseen_V, np.ones((len(unseen_X), self.n_components)))
+                for x, h in zip(unseen_X, unseen_H):
+                    self.H_dict_[x] = h
+                del unseen_H
+            del unseen_X, unseen_V
+        else: # If it is the first batch, call _init_vars to init unq_X, unq_V
+            unq_X, unq_V, lookup = self._init_vars(X)
+
+        unq_H = self._get_H(unq_X)
+        # Update the activations unq_H
+        unq_H = _multiplicative_update_h(
+            unq_V, self.W_, unq_H,
+            epsilon=1e-3, max_iter=self.max_iter_e_step,
+            rescale_W=self.rescale_W,
+            gamma_shape_prior=self.gamma_shape_prior,
+            gamma_scale_prior=self.gamma_scale_prior)
+        # Update the topics self.W_
+        _multiplicative_update_w(
+            unq_V[lookup], self.W_, self.A_, self.B_,
+            unq_H[lookup], self.rescale_W, self.rho_)
+        # Update self.H_dict_ with the learned encoded vectors (activations)
+        self.H_dict_.update(zip(unq_X, unq_H))
+        return self
+
+    def _add_unseen_keys_to_H_dict(self, X):
+        """
+        Add activations of unseen string categories from X to H_dict.
+        """
+        unseen_X = np.setdiff1d(X, np.array([*self.H_dict_]))
+        if unseen_X.size > 0:
+            unseen_V = self.ngrams_count_.transform(unseen_X)
+            if self.add_words:
+                unseen_V2 = self.word_count_.transform(unseen_X)
+                unseen_V = sparse.hstack((unseen_V, unseen_V2), format='csr')
+
+            unseen_H = _rescale_h(
+                unseen_V, np.ones((unseen_V.shape[0], self.n_components)))
+            self.H_dict_.update(zip(unseen_X, unseen_H))
+
+    def transform(self, X):
+        """
+        Return the encoded vectors (activations) H of input strings in X.
+        Given the learnt topics W, the activations H are tuned to fit V = HW.
+
+        Parameters
+        ----------
+        X : array-like, shape (n_samples, ) or (n_samples, 1)
+            The string data to encode.
+
+        Returns
+        -------
+        H : 2-d array, shape (n_samples, n_topics)
+            Transformed input.
+        """
+        # Check input data shape
+        X = np.asarray(X)
+        assert X.ndim == 1 or (X.ndim == 2 and X.shape[1] == 1), f"ERROR:\
+        shape {X.shape} of input array is not supported."
+        if X.ndim == 2:
+            X = X[:, 0]
+        # Check if first item has str or np.str_ type
+        assert isinstance(X[0], str), "ERROR: Input data is not string."
+        unq_X = np.unique(X)
+        # Build the n-grams counts matrix V for the string data to encode
+        unq_V = self.ngrams_count_.transform(unq_X)
+        if self.add_words: # Add words counts
+            unq_V2 = self.word_count_.transform(unq_X)
+            unq_V = sparse.hstack((unq_V, unq_V2), format='csr')
+        # Add unseen strings in X to H_dict
+        self._add_unseen_keys_to_H_dict(unq_X)
+        unq_H = self._get_H(unq_X)
+        # Loop over batches
+        for slice in gen_batches(n=unq_H.shape[0],
+                                 batch_size=self.batch_size):
+            # Given the learnt topics W, optimize H to fit V = HW
+            unq_H[slice] = _multiplicative_update_h(
+                unq_V[slice], self.W_, unq_H[slice],
+                epsilon=1e-3, max_iter=100,
+                rescale_W=self.rescale_W,
+                gamma_shape_prior=self.gamma_shape_prior,
+                gamma_scale_prior=self.gamma_scale_prior)
+        # Store and return the encoded vectors of X
+        self.H_dict_.update(zip(unq_X, unq_H))
+        return self._get_H(X)
+
+
+def _rescale_W(W, A):
+    """
+    Rescale the topics W to have a L1-norm equal to 1.
+    """
+    s = W.sum(axis=1, keepdims=True)
+    W /= s
+    A /= s
+    return
+
+
+def _multiplicative_update_w(Vt, W, A, B, Ht, rescale_W, rho):
+    """
+    Multiplicative update step for the topics W.
+    """
+    A *= rho
+    A += W * safe_sparse_dot(Ht.T, Vt.multiply(np.dot(Ht, W) ** -1))
+    B *= rho
+    B += Ht.sum(axis=0).reshape(-1, 1)
+    np.divide(A, B, out=W)
+    if rescale_W:
+        _rescale_W(W, A)
+    return W, A, B
+
+
+def _rescale_h(V, H):
+    """
+    Rescale the activations H.
+    """
+    epsilon = 1e-10  # in case of a document having length=0
+    H *= np.maximum(epsilon, V.sum(axis=1).A)
+    H /= H.sum(axis=1, keepdims=True)
+    return H
+
+
+def _multiplicative_update_h(Vt, W, Ht, epsilon=1e-3, max_iter=10,
+                             rescale_W=False,
+                             gamma_shape_prior=1.1, gamma_scale_prior=1.):
+    """
+    Multiplicative update step for the activations H.
+    """
+    if rescale_W:
+        WT1 = 1 + 1 / gamma_scale_prior
+        W_WT1 = W / WT1
+    else:
+        WT1 = np.sum(W, axis=1) + 1 / gamma_scale_prior
+        W_WT1 = W / WT1.reshape(-1, 1)
+    const = (gamma_shape_prior - 1) / WT1
+    squared_epsilon = epsilon**2
+    for vt, ht in zip(Vt, Ht):
+        vt_ = vt.data
+        idx = vt.indices
+        W_WT1_ = W_WT1[:, idx]
+        W_ = W[:, idx]
+        squared_norm = 1
+        for n_iter_ in range(max_iter):
+            if squared_norm <= squared_epsilon:
+                break
+            aux = np.dot(W_WT1_, vt_ / np.dot(ht, W_))
+            ht_out = ht * aux + const
+            squared_norm = np.dot(
+                ht_out - ht, ht_out - ht) / np.dot(ht, ht)
+            ht[:] = ht_out
+    return Ht
+
+
+def batch_lookup(lookup, n=1):
+    """ Make batches of the lookup array. """
+    len_iter = len(lookup)
+    for idx in range(0, len_iter, n):
+        indices = lookup[slice(idx, min(idx + n, len_iter))]
+        unq_indices = np.unique(indices)
+        yield (unq_indices, indices)
+
+
+def get_kmeans_prototypes(X, n_prototypes, hashing_dim=128,
+                          ngram_range=(2, 4), sparse=False,
+                          sample_weight=None, random_state=None):
+    """
+    Computes prototypes based on:
+      - dimensionality reduction (via hashing n-grams)
+      - k-means clustering
+      - nearest neighbor
+    """
+    vectorizer = HashingVectorizer(analyzer='char', norm=None,
+                                   alternate_sign=False,
+                                   ngram_range=ngram_range,
+                                   n_features=hashing_dim)
+    projected = vectorizer.transform(X)
+    if not sparse:
+        projected = projected.toarray()
+    kmeans = KMeans(n_clusters=n_prototypes, random_state=random_state)
+    kmeans.fit(projected, sample_weight=sample_weight)
+    centers = kmeans.cluster_centers_
+    neighbors = NearestNeighbors()
+    neighbors.fit(projected)
+    indexes_prototypes = np.unique(neighbors.kneighbors(centers, 1)[-1])
+    return np.sort(X[indexes_prototypes])
diff --git a/dirty_cat/similarity_encoder.py b/dirty_cat/similarity_encoder.py
index aaad7ff..53562c6 100644
--- a/dirty_cat/similarity_encoder.py
+++ b/dirty_cat/similarity_encoder.py
@@ -1,3 +1,17 @@
+"""
+Similarity encoding of string arrays.
+The principle is as follows:
+    1. Given an input string array X = [x1, ..., xn] with k unique categories
+       [c1, ..., ck] and a similarity measure sim(s1, s2) between strings,
+       we define the encoded vector of xi as [sim(xi, c1), ... , sim(xi, ck)].
+       Similarity encoding of X results in a matrix with shape (n, k) that
+       captures morphological similarities between string entries.
+    2. To avoid dealing with high-dimensional encodings when k is high, we can
+       use d << k prototypes [p1, ..., pd] with which similarities will be
+       computed:  xi -> [sim(xi, p1), ..., sim(xi, pd)]. These prototypes can
+       be automatically sampled from the input data (most frequent categories,
+       KMeans) or provided by the user.
+"""
 import warnings
 
 import numpy as np
@@ -273,7 +287,9 @@ class SimilarityEncoder(OneHotEncoder):
         return values[:self.n_prototypes]
 
     def fit(self, X, y=None):
-        """Fit the CategoricalEncoder to X.
+        """
+        Fit the SimilarityEncoder to X.
+
         Parameters
         ----------
         X : array-like, shape [n_samples, n_features]
@@ -389,7 +405,8 @@ class SimilarityEncoder(OneHotEncoder):
         return self
 
     def transform(self, X, fast=True):
-        """Transform X using specified encoding scheme.
+        """
+        Transform X using specified encoding scheme.
 
         Parameters
         ----------
diff --git a/doc/conf.py b/doc/conf.py
index 757b608..1d33a35 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -229,7 +229,7 @@ sphinx_gallery_conf = {
         'dirty_cat': None,
         'numpy': 'http://docs.scipy.org/doc/numpy',
         'scipy': 'http://docs.scipy.org/doc/scipy/reference',
-        'pandas': 'http://pandas.pydata.org/pandas-docs/stable/',
+        'pandas': 'http://pandas.pydata.org/pandas-docs/stable',
         #'seaborn': 'http://seaborn.pydata.org/',
         #'matplotlib': 'http://matplotlib.org/',
         'sklearn': 'http://scikit-learn.org/0.20',
diff --git a/doc/index.rst b/doc/index.rst
index fb9fe6a..1dce332 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -17,6 +17,10 @@ The :class:`SimilarityEncoder` is a drop-in replacement for
 If speed and scalability are an issue, the :class:`MinHashEncoder` provides
 a fast encoding method.
 
+If interpretability is important,the :class:`GapEncoder` is a good
+alternative, as it can be interpreted as one-hot encoding, where each encoding
+dimension corresponds to a topic that summarizes the substrings captured.
+
 For a detailed description of the problem of encoding dirty categorical data,
 see `Similarity encoding for learning with dirty categorical variables
 <https://hal.inria.fr/hal-01806175>`_ [1]_ and `Encoding high-cardinality
@@ -47,6 +51,7 @@ Encoders
    :template: class.rst
    :nosignatures:
 
+   GapEncoder
    MinHashEncoder
    SimilarityEncoder
    TargetEncoder
diff --git a/examples/02_fit_predict_plot_employee_salaries.py b/examples/02_fit_predict_plot_employee_salaries.py
index cb05105..e831533 100644
--- a/examples/02_fit_predict_plot_employee_salaries.py
+++ b/examples/02_fit_predict_plot_employee_salaries.py
@@ -55,7 +55,8 @@ clean_columns = {
 # We then choose the categorical encoding methods we want to benchmark
 # and the dirty categorical variable:
 
-encoding_methods = ['one-hot', 'target', 'similarity', 'minhash']
+encoding_methods = ['one-hot', 'target', 'similarity', 'minhash',
+                    'gap']
 dirty_column = 'employee_position_title'
 #########################################################################
 
@@ -67,13 +68,15 @@ dirty_column = 'employee_position_title'
 
 from sklearn.preprocessing import FunctionTransformer
 from sklearn.preprocessing import OneHotEncoder
-from dirty_cat import SimilarityEncoder, TargetEncoder, MinHashEncoder
+from dirty_cat import SimilarityEncoder, TargetEncoder, MinHashEncoder,\
+    GapEncoder
 
 encoders_dict = {
     'one-hot': OneHotEncoder(handle_unknown='ignore', sparse=False),
     'similarity': SimilarityEncoder(similarity='ngram'),
     'target': TargetEncoder(handle_unknown='ignore'),
     'minhash': MinHashEncoder(n_components=100),
+    'gap': GapEncoder(n_components=100),
     'numerical': FunctionTransformer(None)}
 
 # We then create a function that takes one key of our ``encoders_dict``,
diff --git a/examples/03_fit_predict_plot_midwest_survey.py b/examples/03_fit_predict_plot_midwest_survey.py
index e1dd715..eb6e8fe 100644
--- a/examples/03_fit_predict_plot_midwest_survey.py
+++ b/examples/03_fit_predict_plot_midwest_survey.py
@@ -64,12 +64,14 @@ y = df[target_column].values.ravel()
 # -------------------------------------------
 # We first import the right encoders to transform our clean/dirty data:
 from sklearn.preprocessing import FunctionTransformer, OneHotEncoder
-from dirty_cat import SimilarityEncoder, MinHashEncoder
+from dirty_cat import SimilarityEncoder, MinHashEncoder,\
+    GapEncoder
 
 encoder_dict = {
     'one-hot': OneHotEncoder(handle_unknown='ignore', sparse=False),
     'similarity': SimilarityEncoder(similarity='ngram'),
     'minhash': MinHashEncoder(),
+    'gap': GapEncoder(),
     'num': FunctionTransformer(None)
 }
 ##############################################################################
@@ -111,7 +113,7 @@ from sklearn.model_selection import StratifiedKFold
 
 cv = StratifiedKFold(n_splits=3, random_state=12, shuffle=True)
 all_scores = {}
-for method in ['one-hot', 'similarity', 'minhash']:
+for method in ['one-hot', 'similarity', 'minhash', 'gap']:
     pipeline = make_pipeline(method)
     # Now predict the census region of each participant
     scores = cross_val_score(pipeline, df, y, cv=cv)
diff --git a/examples/04_feature_interpretation_gap_encoder.py b/examples/04_feature_interpretation_gap_encoder.py
new file mode 100644
index 0000000..2fddfc8
--- /dev/null
+++ b/examples/04_feature_interpretation_gap_encoder.py
@@ -0,0 +1,89 @@
+"""
+Feature interpretation with the GapEncoder
+==========================================
+
+We illustrate here how categorical encodings obtained with the GapEncoder
+can be interpreted in terms of latent topics. We use as example the
+`employee salaries <https://catalog.data.gov/dataset/employee-salaries-2016>`_
+dataset, and encode the column *Employee Position Title*, that contains dirty
+categorical data.
+
+"""
+
+################################################################################
+# Data Importing
+# --------------
+#
+# We first download the dataset:
+from dirty_cat.datasets import fetch_employee_salaries
+employee_salaries = fetch_employee_salaries()
+print(employee_salaries['DESCR'])
+
+
+################################################################################
+# Then we load it:
+import pandas as pd
+df = employee_salaries['data']
+
+################################################################################
+# Now, we retrieve the dirty column to encode:
+dirty_column = 'employee_position_title'
+X_dirty = df[dirty_column]
+print(X_dirty.head(), end='\n\n')
+print(f'Number of dirty entries = {len(X_dirty)}')
+
+################################################################################
+# Encoding dirty job titles
+# -------------------------
+#
+# We first create an instance of the GapEncoder with n_components=10:
+from dirty_cat import GapEncoder
+enc = GapEncoder(n_components=10, random_state=42)
+
+################################################################################
+# Then we fit the model on the dirty categorical data and transform it to
+# obtain encoded vectors of size 10:
+X_enc = enc.fit_transform(X_dirty)
+print(f'Shape of encoded vectors = {X_enc.shape}')
+
+################################################################################
+# Interpreting encoded vectors
+# ----------------------------
+#
+# The GapEncoder can be understood as a continuous encoding on a set of latent
+# topics estimated from the data. The latent topics are built by
+# capturing combinations of substrings that frequently co-occur, and encoded
+# vectors correspond to their activations.
+# To interpret these latent topics, we select for each of them a few labels
+# from the input data with the highest activations.
+# In the example below we select 3 labels to summarize each topic.
+
+topic_labels = enc.get_feature_names(n_labels=3)
+for k in range(len(topic_labels)):
+    labels = topic_labels[k]
+    print(f'Topic n°{k}: {labels}')
+
+################################################################################
+# As expected, topics capture labels that frequently co-occur. For instance,
+# the labels *firefighter*, *rescuer*, *rescue* appear together in
+# *Firefigther/Rescuer III*, or *Fire/Rescue Lieutenant*.
+# 
+# This enables us to understand the encoding of different samples
+
+import matplotlib.pyplot as plt
+encoded_labels = enc.transform(X_dirty[:20])
+plt.figure(figsize=(8,10))
+plt.imshow(encoded_labels)
+plt.xlabel('Latent topics', size=12)
+plt.xticks(range(0, 10), labels=topic_labels, rotation=50, ha='right')
+plt.ylabel('Data entries', size=12)
+plt.yticks(range(0, 20), labels=X_dirty[:20])
+plt.colorbar().set_label(label='Topic activations', size=12)
+plt.tight_layout()
+plt.show()
+
+################################################################################
+# As we can see, each dirty category encodes on a small number of topics,
+# These can thus be reliably used to summarize each topic, which are in
+# effect latent categories captured from the data.
+
diff --git a/examples/04_dimension_reduction_and_performance.py b/examples/06_dimension_reduction_and_performance.py
similarity index 97%
rename from examples/04_dimension_reduction_and_performance.py
rename to examples/06_dimension_reduction_and_performance.py
index c5bb36f..d8b4a66 100644
--- a/examples/04_dimension_reduction_and_performance.py
+++ b/examples/06_dimension_reduction_and_performance.py
@@ -6,6 +6,10 @@ Here we discuss how to apply efficiently SimilarityEncoder to larger
 datasets: reducing the number of reference categories to "prototypes",
 either chosen as the most frequent categories, or with kmeans clustering.
 
+Note that the :class:`GapEncoder` naturally does data reduction and comes
+with online estimation. As a result is it more scalable than the
+SimilarityEncoder, and should be prefered in large-scale settings.
+
 """
 # Avoid the warning in scikit-learn's LogisticRegression for the change
 # in the solver


Test patch:
diff --git a/dirty_cat/test/test_gap_encoder.py b/dirty_cat/test/test_gap_encoder.py
new file mode 100644
index 0000000..db3f926
--- /dev/null
+++ b/dirty_cat/test/test_gap_encoder.py
@@ -0,0 +1,122 @@
+import time
+import numpy as np
+import pytest
+
+from sklearn.datasets import fetch_20newsgroups
+from dirty_cat import GapEncoder
+
+@pytest.mark.parametrize("hashing, init, analyzer, add_words", [
+    (False, 'k-means++', 'word', True),
+    (True, 'random', 'char', False),
+    (True, 'k-means', 'char_wb', True)
+])
+def test_gap_encoder(hashing, init, analyzer, add_words, n_samples=70):
+    X_txt = fetch_20newsgroups(subset='train')['data']
+    X = X_txt[:n_samples]
+    n_components = 10
+    # Test output shape
+    encoder = GapEncoder(
+        n_components=n_components, hashing=hashing, init=init,
+        analyzer=analyzer, add_words=add_words,
+        random_state=42, rescale_W=True)
+    encoder.fit(X)
+    y = encoder.transform(X)
+    assert y.shape == (n_samples, n_components), str(y.shape)
+    assert len(set(y[0])) == n_components
+
+    # Test L1-norm of topics W.
+    l1_norm_W = np.abs(encoder.W_).sum(axis=1)
+    np.testing.assert_array_almost_equal(
+        l1_norm_W, np.ones(n_components))
+
+    # Test same seed return the same output
+    encoder = GapEncoder(
+        n_components=n_components, hashing=hashing, init=init,
+        analyzer=analyzer, add_words=add_words,
+        random_state=42)
+    encoder.fit(X)
+    y2 = encoder.transform(X)
+    np.testing.assert_array_equal(y, y2)
+    return
+
+
+def test_input_type():
+    # Numpy array
+    X = np.array(['alice', 'bob'])
+    enc = GapEncoder(n_components=2, random_state=42)
+    X_enc_array = enc.fit_transform(X)
+    # List
+    X = ['alice', 'bob']
+    enc = GapEncoder(n_components=2, random_state=42)
+    X_enc_list = enc.fit_transform(X)
+    # Check if the encoded vectors are the same
+    np.testing.assert_array_equal(X_enc_array, X_enc_list)
+    return
+
+
+def test_partial_fit(n_samples=70):
+    X_txt = fetch_20newsgroups(subset='train')['data']
+    X = X_txt[:n_samples]
+    # Gap encoder with fit on one batch
+    enc = GapEncoder(random_state=42, batch_size=n_samples, max_iter=1)
+    X_enc = enc.fit_transform(X)
+    # Gap encoder with partial fit
+    enc = GapEncoder(random_state=42)
+    enc.partial_fit(X)
+    X_enc_partial = enc.transform(X)
+    # Check if the encoded vectors are the same
+    np.testing.assert_array_equal(X_enc, X_enc_partial)
+    return
+
+
+def test_get_feature_names(n_samples=70):
+    X_txt = fetch_20newsgroups(subset='train')['data']
+    X = X_txt[:n_samples]
+    enc = GapEncoder()
+    enc.fit(X)
+    topic_labels = enc.get_feature_names()
+    # Check number of labels
+    assert len(topic_labels) == enc.n_components
+    return
+    
+
+def profile_encoder(Encoder, init):
+    # not an unit test
+    
+    from dirty_cat import datasets
+    employee_salaries = datasets.fetch_employee_salaries()
+    data = employee_salaries['data']
+    X = data['employee_position_title'].tolist()
+    t0 = time.time()
+    encoder = Encoder(n_components=50, init=init)
+    encoder.fit(X)
+    y = encoder.transform(X)
+    assert y.shape == (len(X), 50)
+    eta = time.time() - t0
+    return eta
+
+
+if __name__ == '__main__':
+    print('start test_gap_encoder')
+    test_gap_encoder(True, 'k-means++', 'char', False)
+    print('test_gap_encoder passed')
+    print('start test_input_type')
+    test_input_type()
+    print('test_input_type passed')
+    print('start test_partial_fit')
+    test_partial_fit()
+    print('test_partial_fit passed')
+    print('start test_get_feature_names')
+    test_get_feature_names()
+    print('test_get_feature_names passed')
+    
+    for _ in range(3):
+        print('time profile_encoder(GapEncoder, init="k-means++")')
+        print("{:.4} seconds".format(profile_encoder(
+            GapEncoder, init='k-means++')))
+    for _ in range(3):
+        print('time profile_encoder(GapEncoder, init="k-means")')
+        print("{:.4} seconds".format(profile_encoder(
+            GapEncoder, init='k-means')))
+    
+    print('Done')
diff --git a/dirty_cat/test/test_minhash_encoder.py b/dirty_cat/test/test_minhash_encoder.py
index b1afa02..c968816 100644
--- a/dirty_cat/test/test_minhash_encoder.py
+++ b/dirty_cat/test/test_minhash_encoder.py
@@ -148,12 +148,14 @@ if __name__ == '__main__':
 
     for _ in range(3):
         print('time profile_encoder(MinHashEncoder, hashing=fast)')
-        print("{:.4} seconds".format(profile_encoder(MinHashEncoder, hashing='fast')))
+        print("{:.4} seconds".format(
+            profile_encoder(MinHashEncoder, hashing='fast')))
     for _ in range(3):
         print('time profile_encoder(MinHashEncoder, hashing=fast) with minmax')
         print("{:.4} seconds".format(profile_encoder(MinHashEncoder,
                                      hashing='fast', minmax_hash=True)))
     print('time profile_encoder(MinHashEncoder, hashing=murmur)')
-    print("{:.4} seconds".format(profile_encoder(MinHashEncoder, hashing='murmur')))
+    print("{:.4} seconds".format(
+        profile_encoder(MinHashEncoder, hashing='murmur')))
 
     print('Done')