Repository: dirty-cat/dirty_cat
Pull request #267 (2022-06-29 17:41:30)

ENH MinHash parallel
Compute the min hash transform method in parallel, as suggested by @alexis-cvetkov.

We no longer use the self.hash_dict attribute, so the fit method does nothing now.

LeoGrin: I don't understand why the coverage has changed. It seems that the function called by joblib.Parallel (`compute_hash`) is not counted in the coverage, but I may be reading codecov wrong.
GaelVaroquaux: I worry that the self.hash_dict was really useful to speed things up by avoiding recomputation of repeated entries
LeoGrin: We use np.unique, so for one transform we don't recompute repeated entries.
Using self.hash_dict would indeed speed things up if we transform several
inputs with common entries, using the same encoder. Do you think that's
likely to happen / worth the additional memory usage ?

On Wed, Jun 29, 2022, 20:59 Gael Varoquaux ***@***.***> wrote:

> I worry that the self.hash_dict was really useful to speed things up by
> avoiding recomputation of repeated entries
>
> —
> Reply to this email directly, view it on GitHub
> <https://github.com/dirty-cat/dirty_cat/pull/267#issuecomment-1170376511>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AK46V2FGQSUNMHVILR5MBODVRSMIDANCNFSM52GSCHYA>
> .
> You are receiving this because you authored the thread.Message ID:
> ***@***.***>
>

GaelVaroquaux: > Do you think that's likely to happen / worth the additional memory usage ?

Yes: it's very frequent. Typically, the entries are repeated many times.

LeoGrin: @GaelVaroquaux just want to make sure I understood what you were saying before putting the hash_dict back in the code.
LeoGrin: Following discussion with @GaelVaroquaux : using the same encoder to transform several Xs may happen in online learning settings, for instance with a big X.
I've put the hash_dict back in
GaelVaroquaux: > Benchmarking really quickly on my Mac M1 with 8 cores, the rows version is about twice faster than the batched version, but we should do a more serious benchmark.

On a dataset of which size (number of rows)?

jjerphan: As discussed with @jovan-stojanovic privately regarding understanding the current implementation performance: one can profile it using profilers like [`py-spy`](https://github.com/benfred/py-spy) and visually search for bottlenecks using [Speed-scope](https://github.com/jlfwong/speedscope)

In that case, one can get a `py-spy` profile for Speedcope using: 
```bash
py-spy record --native -o py-spy-profile.svg -f speedscope -- python ./script.py
```

where `script.py` is the script using `MinHash` on a suitable dataset for `dirty_cat`.

If `MinHash` relies on `joblib` for parallelization, it might be better to use a sole thread (by setting `n_jobs=1`) so as to better understand the execution when inspecting the profile.

After having uploaded the profile on [Speedscope's UI](https://speedscope.app/), you will be able to have access to various representations:
 - `"Left Heavy"` to get a flame graph recapitulating the total runtime of each functions, we the heavier first
 - `"Time Order"` to understand the order of function call in their context.
 - `"Sandwich"` for fine-grained information for each interfaces
jovan-stojanovic: The changelog test apparently works fine :) You should add the PR number to CHANGES.rst
LeoGrin: - I've run the benchmark (thanks Lilian) on Margaret, and it confirms that the batched version seems to always be better. I've also added to the benchmark a `batch_per_job` parameter, but it had no influence on the results.

- I've put the minhashencoder code with all options in the bench_minhash function.

- In the original minhashencoder file, I've only kept the batched version, and removed the parameter to choose.

Ready for review!
LilianBoulard: Just waiting for Jovan's review, and we'll be ready to merge :)

Patch:
diff --git a/.gitignore b/.gitignore
index 75917d3..dc2445f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -62,3 +62,6 @@ doc/RELEASE_PROCESS.rst
 doc/CONTRIBUTING.rst
 .DS_Store
 
+# Pkl files for benchmarks
+benchmarks/*.pkl
+
diff --git a/CHANGES.rst b/CHANGES.rst
index 41bb3f0..3c1af98 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -9,8 +9,8 @@ Major changes
 * New experimental feature: joining tables using :func:`fuzzy_join` by approximate key matching. Matches are based
   on string similarities and the nearest neighbors matches are found for each category. :pr:`291` by :user:`Jovan Stojanovic <jovan-stojanovic>` and :user:`Leo Grinsztajn <LeoGrin>`
 * **datasets.fetching**: contains a new function :func:`fetch_world_bank_indicator` that can be used to download any indicator from the World Bank Open Data platform. It only needs the indicator ID that can be found on the website. :pr:`291` by :user:`Jovan Stojanovic <jovan-stojanovic>`
-* Unnecessary API has been made private: everything (files, functions, classes)
-  starting with an underscore shouldn't be imported in your code. :pr:`331` by :user:`Lilian Boulard <LilianBoulard>`
+* Unnecessary API has been made private: everything (files, functions, classes) starting with an underscore shouldn't be imported in your code. :pr:`331` by :user:`Lilian Boulard <LilianBoulard>`
+* The MinHashEncoder now supports a `n_jobs` parameter to parallelize the hashes computation. :pr:`267` by `Leo Grinsztajn <LeoGrin>` and `Lilian Boulard <LilianBoulard>`.
 
 Minor changes
 -------------
diff --git a/benchmarks/bench_minhash.py b/benchmarks/bench_minhash.py
new file mode 100644
index 0000000..bbffa6f
--- /dev/null
+++ b/benchmarks/bench_minhash.py
@@ -0,0 +1,460 @@
+#################
+# The MinHashEncoder version used for the benchmark
+# On the main branch, we only kept the best version of the MinHashEncoder
+# which is the batched version
+# with batch_per_job=1
+# (the batch_per_job parameter has no effect on the results)
+#################
+
+from typing import Callable, Collection, Dict, List, Literal, Tuple
+
+import numpy as np
+from joblib import Parallel, delayed, effective_n_jobs
+from sklearn.base import BaseEstimator, TransformerMixin
+from sklearn.utils import gen_even_slices, murmurhash3_32
+
+from dirty_cat._fast_hash import ngram_min_hash
+from dirty_cat._string_distances import get_unique_ngrams
+from dirty_cat._utils import LRUDict, check_input
+
+NoneType = type(None)
+
+
+# Ignore lines too long, as links in the docstring cannot be cut.
+# flake8: noqa: E501
+
+
+class MinHashEncoder(BaseEstimator, TransformerMixin):
+    """
+    Encode string categorical features as a numeric array, minhash method
+    applied to ngram decomposition of strings based on ngram decomposition
+    of the string.
+
+    Parameters
+    ----------
+    n_components : int, default=30
+        The number of dimension of encoded strings. Numbers around 300 tend to
+        lead to good prediction performance, but with more computational cost.
+    ngram_range : typing.Tuple[int, int], default=(2, 4)
+        The lower and upper boundary of the range of n-values for different
+        n-grams to be extracted. All values of n such that min_n <= n <= max_n.
+        will be used.
+    hashing : typing.Literal["fast", "murmur"], default=fast
+        Hashing function. fast is faster but
+        might have some concern with its entropy.
+    minmax_hash : bool, default=False
+        if True, return min hash and max hash concatenated.
+    handle_missing : typing.Literal["error", "zero_impute"], default=zero_impute
+        Whether to raise an error or encode missing values (NaN) with
+        vectors filled with zeros.
+    batch: bool, default=False
+        If False, parallelize the computation of the hash of each unique
+        element. If True, parallelize the computation of the hash of each
+        batch of elements.
+    batch_per_job: int, default=1
+        Number of batches to be processed in each job.
+    n_jobs : int, default=None
+        The number of jobs to run in parallel.
+        The hash computations for all unique elements are parallelized.
+        None means 1 unless in a
+        `joblib.parallel_backend context <https://joblib.readthedocs.io/en/latest/parallel.html>`_.
+        -1 means using all processors.
+        See `Scikit-learn Glossary <https://scikit-learn.org/stable/glossary.html#term-n_jobs>`_
+        for more details.
+
+    Attributes
+    ----------
+    hash_dict_ : LRUDict
+        Computed hashes.
+    
+    Examples
+    --------
+    >>> enc = MinHashEncoder(n_components=5)
+
+    Let's encode the following non-normalized data:
+
+    >>> X = [['paris, FR'], ['Paris'], ['London, UK'], ['London']]
+
+    >>> enc.fit(X)
+    MinHashEncoder()
+
+    The encoded data with 5 components are:
+
+    >>> enc.transform(X)
+    array([[-1.78337518e+09, -1.58827021e+09, -1.66359234e+09,
+            -1.81988679e+09, -1.96259387e+09],
+           [-8.48046971e+08, -1.76657887e+09, -1.55891205e+09,
+            -1.48574446e+09, -1.68729890e+09],
+           [-1.97582893e+09, -2.09500033e+09, -1.59652117e+09,
+            -1.81759383e+09, -2.09569333e+09],
+           [-1.97582893e+09, -2.09500033e+09, -1.53072052e+09,
+            -1.45918266e+09, -1.58098831e+09]])
+
+    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).
+
+    """
+
+    hash_dict_: LRUDict
+
+    _capacity: int = 2**10
+
+    def __init__(
+        self,
+        n_components: int = 30,
+        ngram_range: Tuple[int, int] = (2, 4),
+        hashing: Literal["fast", "murmur"] = "fast",
+        minmax_hash: bool = False,
+        handle_missing: Literal["error", "zero_impute"] = "zero_impute",
+        batch: bool = False,
+        batch_per_job: int = 1,
+        n_jobs: int = None,
+    ):
+        self.ngram_range = ngram_range
+        self.n_components = n_components
+        self.hashing = hashing
+        self.minmax_hash = minmax_hash
+        self.handle_missing = handle_missing
+        self.batch = batch
+        self.batch_per_job = batch_per_job
+        self.n_jobs = n_jobs
+
+    def _more_tags(self) -> Dict[str, List[str]]:
+        """
+        Used internally by sklearn to ease the estimator checks.
+        """
+        return {"X_types": ["categorical"]}
+
+    def _get_murmur_hash(self, string: str) -> np.array:
+        """
+        Encode a string using murmur hashing function.
+
+        Parameters
+        ----------
+        string : str
+            The string to encode.
+
+        Returns
+        -------
+        ndarray of shape (n_components, )
+            The encoded string.
+        """
+        min_hashes = np.ones(self.n_components) * np.infty
+        grams = get_unique_ngrams(string, self.ngram_range)
+        if len(grams) == 0:
+            grams = get_unique_ngrams(" Na ", self.ngram_range)
+        for gram in grams:
+            hash_array = np.array(
+                [
+                    murmurhash3_32("".join(gram), seed=d, positive=True)
+                    for d in range(self.n_components)
+                ]
+            )
+            min_hashes = np.minimum(min_hashes, hash_array)
+        return min_hashes / (2**32 - 1)
+
+    def _get_fast_hash(self, string: str) -> np.array:
+        """
+        Encode a string with fast hashing function.
+        fast hashing supports both min_hash and minmax_hash encoding.
+
+        Parameters
+        ----------
+        string : str
+            The string to encode.
+
+        Returns
+        -------
+        ndarray of shape (n_components, )
+            The encoded string, using specified encoding scheme.
+        """
+        if self.minmax_hash:
+            return np.concatenate(
+                [
+                    ngram_min_hash(string, self.ngram_range, seed, return_minmax=True)
+                    for seed in range(self.n_components // 2)
+                ]
+            )
+        else:
+            return np.array(
+                [
+                    ngram_min_hash(string, self.ngram_range, seed)
+                    for seed in range(self.n_components)
+                ]
+            )
+
+    def _compute_hash(
+        self, string: str, hash_func: Callable[[str], np.ndarray]
+    ) -> np.ndarray:
+        """Function called to compute the hash of a string.
+
+        Check if the string is in the hash dictionary, if not, scompute the hash using
+        the specified hashing function and add it to the dictionary.
+
+        Parameters
+        ----------
+        string : str
+            The string to encode.
+        hash_func : callable
+            Hashing function to use on the string.
+
+        Returns
+        -------
+        np.array of shape (n_components, )
+            The encoded string, using specified encoding scheme.
+        """
+        if string not in self.hash_dict_:
+            if string == "NAN":  # true if x is a missing value
+                self.hash_dict_[string] = np.zeros(self.n_components)
+            else:
+                self.hash_dict_[string] = hash_func(string)
+        return self.hash_dict_[string]
+
+    def _compute_hash_batched(
+        self, batch: Collection[str], hash_func: Callable[[str], np.ndarray]
+    ):
+        """Function called to compute the hashes of a batch of strings.
+
+        Check if the string is in the hash dictionary, if not, compute the hash using
+        the specified hashing function and add it to the dictionary.
+
+        Parameters
+        ----------
+        batch : iterable of str
+            The batch of strings to encode.
+        hash_func : callable
+            Hashing function to use on the string.
+
+        Returns
+        -------
+        np.array of shape (n_samples, n_components)
+            The encoded strings, using specified encoding scheme.
+        """
+        res = np.zeros((len(batch), self.n_components))
+        for i, string in enumerate(batch):
+            if string not in self.hash_dict_:
+                if string == "NAN":  # true if x is a missing value
+                    self.hash_dict_[string] = np.zeros(self.n_components)
+                else:
+                    self.hash_dict_[string] = hash_func(string)
+            res[i] = self.hash_dict_[string]
+        return res
+
+    def fit(self, X, y=None) -> "MinHashEncoder":
+        """
+        Fit the MinHashEncoder to X. In practice, just initializes a dictionary
+        to store encodings to speed up computation.
+
+        Parameters
+        ----------
+        X : array-like, shape (n_samples, ) or (n_samples, 1)
+            The string data to encode. Only here for compatibility.
+        y : None
+            Unused, only here for compatibility.
+
+        Returns
+        -------
+        MinHashEncoder
+            The fitted MinHashEncoder instance.
+        """
+        if self.hashing not in ["fast", "murmur"]:
+            raise ValueError(
+                f"Got hashing={self.hashing!r}, "
+                'but expected any of {"fast", "murmur"}. '
+            )
+        if self.handle_missing not in ["error", "zero_impute"]:
+            raise ValueError(
+                f"Got handle_missing={self.handle_missing!r}, but expected "
+                'any of {"error", "zero_impute"}. '
+            )
+        self.hash_dict_ = LRUDict(capacity=self._capacity)
+        return self
+
+    def transform(self, X) -> np.array:
+        """
+        Transform X using specified encoding scheme.
+
+        Parameters
+        ----------
+        X : array-like, shape (n_samples, ) or (n_samples, 1)
+            The string data to encode.
+
+        Returns
+        -------
+        ndarray of shape (n_samples, n_components)
+            Transformed input.
+        """
+        X = check_input(X)
+        if self.minmax_hash:
+            if self.n_components % 2 != 0:
+                raise ValueError(
+                    "n_components should be even when using"
+                    f"minmax_hash encoding, got {self.n_components}"
+                )
+        if self.hashing == "murmur":
+            if self.minmax_hash:
+                raise ValueError(
+                    "minmax_hash encoding is not supported"
+                    "with the murmur hashing function"
+                )
+        if self.handle_missing not in ["error", "zero_impute"]:
+            raise ValueError(
+                "handle_missing should be either "
+                f"'error' or 'zero_impute', got {self.handle_missing!r}"
+            )
+
+        # Handle missing values
+        missing_mask = (
+            ~(X == X)  # Find np.nan
+            | (X == None)  # Find None. Note: `X is None` doesn't work.
+            | (X == "")  # Find empty strings
+        )
+
+        if missing_mask.any():  # contains at least one missing value
+            if self.handle_missing == "error":
+                raise ValueError(
+                    "Found missing values in input data; set "
+                    "handle_missing='zero_impute' to encode with missing values"
+                )
+            elif self.handle_missing == "zero_impute":
+                # NANs will be replaced by zeroes in _compute_hash
+                X[missing_mask] = "NAN"
+
+        if self.hashing == "fast":
+            hash_func = self._get_fast_hash
+        elif self.hashing == "murmur":
+            hash_func = self._get_murmur_hash
+        else:
+            raise ValueError(
+                "Hashing function should be either 'fast' or 'murmur', "
+                f"got {self.hashing!r}"
+            )
+
+        # Compute the hashes for unique values
+        unique_x, indices_x = np.unique(X, return_inverse=True)
+        n_jobs = effective_n_jobs(self.n_jobs)
+
+        if self.batch:
+            unique_x_trans = Parallel(n_jobs=n_jobs)(
+                delayed(self._compute_hash_batched)(
+                    unique_x[idx_slice],
+                    hash_func,
+                )
+                for idx_slice in gen_even_slices(len(unique_x), 
+                n_jobs * self.batch_per_job)
+            )
+            # Match the hashes of the unique value to the original values
+            X_out = np.concatenate(unique_x_trans)[indices_x].reshape(
+                len(X), X.shape[1] * self.n_components
+            )
+        else:
+            unique_x_trans = Parallel(n_jobs=n_jobs)(
+                delayed(self._compute_hash)(x, hash_func) for x in unique_x
+            )
+            # Match the hashes of the unique value to the original values
+            X_out = np.stack(unique_x_trans)[indices_x].reshape(
+                len(X), X.shape[1] * self.n_components
+            )
+
+        return X_out.astype(np.float64)  # The output is an int32 before conversion
+
+
+import pickle
+import pandas as pd
+import seaborn as sns
+import matplotlib.pyplot as plt
+from argparse import ArgumentParser
+
+from pathlib import Path
+
+from utils import monitor, parse_func_repr, find_result, default_parser
+from dirty_cat.tests.utils import generate_data
+
+benchmark_name = "minhash_batch_comparison"
+
+
+@monitor(
+    memory=True,
+    time=True,
+    parametrize={
+        "dataset_size": ["medium"],
+        "batched": [True, False],
+        "n_jobs": [1, 4, 8, 16, 32, 64],
+        "batch_per_job": [1, 2, 4],
+    },
+    save_as=benchmark_name,
+    repeat=10,
+)
+def benchmark(
+    dataset_size: str,
+    batched: bool,
+    n_jobs: int,
+    batch_per_job: int,
+) -> None:
+    X = data[dataset_size]
+    MinHashEncoder(batch=batched, n_jobs=n_jobs, 
+    batch_per_job=batch_per_job).fit(X).transform(X)
+
+
+def plot(res: pd.DataFrame):
+    sns.set_theme(style="ticks", palette="pastel")
+
+    rows = []
+    for i, ser in res.iterrows():
+        times = eval(str(ser["time"]))
+        memories = eval(str(ser["memory"]))
+        _, _, kwargs = parse_func_repr(ser["call"])
+        for time, memory in zip(times, memories):
+            rows.append((kwargs["batched"], kwargs["batch_per_job"], 
+            kwargs["n_jobs"], time, memory))
+
+    df = pd.DataFrame(rows, columns=["batched", "batch_per_job", 
+    "n_jobs", "time", "memory"])
+
+    # Create a new columns merging batched and batch_per_job
+    # If batch is False, ignore batch_per_job
+    df["config"] = df.apply(
+        lambda row: f"batched={row['batched']}, batch_per_job={row['batch_per_job']}"
+        if row["batched"] == 'True'
+        else "batched=False",
+        axis=1,
+    )
+
+    sns.boxplot(x="n_jobs", y="time", hue="config", data=df)
+    plt.show()
+
+
+if __name__ == "__main__":
+    _args = ArgumentParser(
+        description="Benchmark for the batch feature of the MinHashEncoder.",
+        parents=[default_parser],
+    ).parse_args()
+
+    # Generate the data if not already on disk, and keep them in memory.
+    data = {}  # Will hold the datasets in memory.
+    _data_info = {
+        "small": 10_000,
+        "medium": 100_000,
+    }
+    for name, size in _data_info.items():
+        data_file = Path(f"data_{name}.pkl")
+        if data_file.is_file():
+            with data_file.open("rb") as fl:
+                data.update({name: pickle.load(fl)})
+        else:
+            with data_file.open("wb") as fl:
+                _gen = generate_data(size).reshape(-1, 1)
+                pickle.dump(_gen, fl)
+                data.update({name: _gen})
+
+    if _args.run:
+        df = benchmark()
+    else:
+        result_file = find_result(benchmark_name)
+        df = pd.read_csv(result_file)
+
+    if _args.plot:
+        plot(df)
diff --git a/benchmarks/results/README.md b/benchmarks/results/README.md
new file mode 100644
index 0000000..4858d46
--- /dev/null
+++ b/benchmarks/results/README.md
@@ -0,0 +1,5 @@
+# Results
+
+## Format
+
+Results are saved with the format ``<name>-<YYYYMMDD>-<commit>.csv``.
diff --git a/benchmarks/results/minhash_batch_comparison-20221213-15855d6a5ba287859b8ef3f9c4428d8fd1ec6e0d.csv b/benchmarks/results/minhash_batch_comparison-20221213-15855d6a5ba287859b8ef3f9c4428d8fd1ec6e0d.csv
new file mode 100644
index 0000000..91b4b38
--- /dev/null
+++ b/benchmarks/results/minhash_batch_comparison-20221213-15855d6a5ba287859b8ef3f9c4428d8fd1ec6e0d.csv
@@ -0,0 +1,37 @@
+call,time,memory
+"benchmark(dataset_size=medium, batched=True, n_jobs=1, batch_per_job=1)","[71.7798215909861, 70.72907463507727, 71.03302158997394, 72.48582117399201, 72.24784656590782, 72.2887754640542, 72.58428386412561, 72.27035080105998, 72.6238446731586, 72.47761388099752]","[214.21227550506592, 214.27148246765137, 214.27683544158936, 214.28184509277344, 214.2869691848755, 214.2921543121338, 214.29727840423584, 214.3024024963379, 214.30752658843994, 214.3127727508545]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=1, batch_per_job=2)","[66.2392854259815, 67.91648005112074, 68.2309153759852, 68.10801578010432, 67.818716380978, 67.35926558589563, 68.37697630305775, 68.16613689390942, 68.34200280206278, 68.1619816957973]","[214.1971616744995, 214.20277404785156, 214.20819568634033, 214.21353340148926, 214.21855068206787, 214.22362899780273, 214.22864627838135, 214.23366355895996, 214.23868083953857, 214.2438201904297]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=1, batch_per_job=4)","[66.83258968195878, 67.59531567618251, 68.39756254991516, 66.72520995908417, 65.89418164314702, 66.32795236585662, 66.05194070399739, 66.08375209686346, 65.92125198896974, 66.0887539810501]","[214.1971616744995, 214.2028350830078, 214.20836353302002, 214.21370124816895, 214.20516204833984, 214.208233833313, 214.21124458312988, 214.21425533294678, 214.21726608276367, 214.22241306304932]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=4, batch_per_job=1)","[9.992393533000723, 8.598808948881924, 8.539826513966545, 8.563422998879105, 8.5266357162036, 8.551135550951585, 8.528889795066789, 8.603934412822127, 8.557675641961396, 8.567077609943226]","[214.1971616744995, 237.34835720062256, 260.25675678253174, 243.08962059020996, 265.9920015335083, 288.8993282318115, 311.8058500289917, 334.7131128311157, 357.61919593811035, 254.56191635131836]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=4, batch_per_job=2)","[8.61273738881573, 8.59429851709865, 8.596364103024825, 8.602652452886105, 8.562634913949296, 8.583562623010948, 8.565711832838133, 8.592467660084367, 8.642348730936646, 8.560480661923066]","[214.1972074508667, 237.11753845214844, 260.03616428375244, 237.12703704833984, 260.0396203994751, 282.9640064239502, 305.88242626190186, 260.06193923950195, 282.97411823272705, 305.8922481536865]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=4, batch_per_job=4)","[8.712641555117443, 8.761632398935035, 8.749574074056, 8.722248969133943, 8.705631776014343, 8.72091040504165, 8.714637140976265, 8.740190399810672, 8.724591991864145, 8.68041093996726]","[214.1972074508667, 237.13661670684814, 260.07965183258057, 244.32647895812988, 267.2655735015869, 255.81029510498047, 278.7473306655884, 264.42577743530273, 287.36153984069824, 275.9062042236328]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=8, batch_per_job=1)","[6.419639005092904, 4.446685581002384, 6.0397503359708935, 6.714352491078898, 6.616947069996968, 7.068617801880464, 7.249980088789016, 7.248966332990676, 7.220058102160692, 7.254382906947285]","[214.1972074508667, 237.1735725402832, 260.0879096984863, 283.0056734085083, 305.9242286682129, 260.10417079925537, 283.01622772216797, 305.93353366851807, 328.85325145721436, 277.3041639328003]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=8, batch_per_job=2)","[6.7550203781574965, 6.866674510994926, 7.020281277131289, 7.12872219691053, 6.6853322891984135, 6.924615056952462, 7.0306664479430765, 7.142228520009667, 6.79599877493456, 6.739760363008827]","[214.1972074508667, 237.1434564590454, 260.0799732208252, 260.0921382904053, 283.027792930603, 281.6062488555908, 304.5414276123047, 303.1205167770386, 326.05605697631836, 237.18389701843262]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=8, batch_per_job=4)","[6.626536522060633, 6.592130372067913, 6.653645798098296, 6.58796801394783, 6.388609720859677, 6.582637360086665, 6.700353967025876, 6.0681178260128945, 6.293093926040456, 6.689890488982201]","[214.1972074508667, 237.18165683746338, 247.96743869781494, 259.50123023986816, 271.0041456222534, 285.3778839111328, 298.3162031173706, 312.6886901855469, 328.4975814819336, 343.59184074401855]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=16, batch_per_job=1)","[6.376049616141245, 3.732655407860875, 4.044250387931243, 4.049081849865615, 4.108462610049173, 4.107385375071317, 4.281491765053943, 4.04321722406894, 4.077744923066348, 3.96216689911671]","[214.1972074508667, 237.2104253768921, 237.22570991516113, 260.1607847213745, 260.17352294921875, 283.1133909225464, 281.692232131958, 304.6257629394531, 304.63868618011475, 327.5725860595703]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=16, batch_per_job=2)","[3.864930802024901, 4.318091581808403, 4.157828219002113, 4.0143457821104676, 3.997693278826773, 4.093358827987686, 3.9836495220661163, 4.169352222932503, 4.3794266388285905, 4.120034449035302]","[214.1971616744995, 237.20193195343018, 257.32004737854004, 277.43759059906006, 297.5552864074707, 317.7040796279907, 337.82298469543457, 357.2237071990967, 237.29335498809814, 255.97091674804688]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=16, batch_per_job=4)","[4.159369282890111, 3.9657786379102618, 3.901250201044604, 3.472345104906708, 3.6302966971416026, 3.7497342589776963, 3.493280997965485, 3.6301434109918773, 3.5428597962018102, 3.6531644260976464]","[214.19717693328857, 237.29165744781494, 254.6153507232666, 273.7376022338867, 292.1402778625488, 231.58457279205322, 249.26528453826904, 268.3882761001587, 287.1500988006592, 305.55174446105957]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=32, batch_per_job=1)","[5.809042064938694, 2.531181629979983, 2.5675521921366453, 2.392015988007188, 2.541553267976269, 2.4252031629439443, 2.3804814720060676, 2.501053577987477, 2.4517790980171412, 2.4895756288897246]","[214.1972074508667, 237.29968929290771, 252.39299201965332, 271.79650688171387, 291.1957950592041, 310.59706687927246, 329.99656867980957, 349.39600944519043, 368.7954502105713, 388.19506645202637]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=32, batch_per_job=2)","[2.634293040027842, 2.4802560589741915, 2.5602449481375515, 2.61546884290874, 2.494540259009227, 2.5182681342121214, 2.4511514641344547, 2.5862338228616863, 2.486740163061768, 2.551073926035315]","[214.1971616744995, 237.30239963531494, 260.3798999786377, 278.4223222732544, 298.62550163269043, 318.46828269958496, 338.31251335144043, 260.45961570739746, 237.4181718826294, 260.5137004852295]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=32, batch_per_job=4)","[2.3421540178824216, 2.3203594910446554, 2.341076180106029, 2.3535482508596033, 2.4552976719569415, 2.317839816911146, 2.325829871930182, 2.3533592959865928, 2.31673446088098, 2.370306008029729]","[214.1972074508667, 231.27839279174805, 249.81013011932373, 271.9742259979248, 242.0189971923828, 262.54346656799316, 284.16064643859863, 254.0118818283081, 274.3552665710449, 296.1530656814575]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=64, batch_per_job=1)","[4.569100285880268, 3.4434140380471945, 2.5561872948892415, 2.536384748062119, 2.4932553980033845, 2.5488017410971224, 2.484845104860142, 2.5588361639529467, 2.468256894964725, 2.483591098105535]","[214.1971616744995, 237.39584159851074, 260.4836187362671, 283.5699005126953, 237.4486789703369, 254.8939085006714, 274.7360563278198, 293.8602533340454, 313.34483337402344, 332.10920429229736]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=64, batch_per_job=2)","[2.4020017068833113, 2.4049791598226875, 2.4657973321154714, 2.4932025589514524, 2.4242782338988036, 2.4349299429450184, 2.4682194851338863, 2.424371248111129, 2.4207907870877534, 2.3780359169468284]","[214.1971616744995, 237.4597454071045, 260.712233543396, 237.46996879577637, 260.72054958343506, 239.66408443450928, 262.91400146484375, 251.30936813354492, 274.55061054229736, 297.80428886413574]"
+"benchmark(dataset_size=medium, batched=True, n_jobs=64, batch_per_job=4)","[2.5223827590234578, 2.463605107041076, 2.4444141660351306, 2.503975408151746, 2.453904598020017, 2.4427910877857357, 2.4992811230476946, 2.5613689851015806, 2.7549221641384065, 2.4777374118566513]","[214.1972074508667, 237.7064723968506, 237.831880569458, 242.69738864898682, 264.90131759643555, 254.70816802978516, 268.8981580734253, 263.7617025375366, 285.05506896972656, 237.87235069274902]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=1, batch_per_job=1)","[81.9005116620101, 84.12667609006166, 84.48668759106658, 82.52440749807283, 84.33229907299392, 83.91563078900799, 83.53606450394727, 83.3229243028909, 83.4883614720311, 83.65322434785776]","[214.1972074508667, 214.22509956359863, 214.22706031799316, 214.22735786437988, 214.2282886505127, 214.22845649719238, 214.22850227355957, 214.228609085083, 214.2286548614502, 214.22882270812988]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=1, batch_per_job=2)","[81.63410263694823, 83.94753866200335, 83.83449222310446, 83.85091615584679, 83.60108938603662, 83.74655381008051, 83.70933479792438, 83.65041901101358, 83.94523491105065, 83.7649101971183]","[214.1972532272339, 214.21961975097656, 214.2212142944336, 214.22132110595703, 214.22148895263672, 214.2216567993164, 214.2217025756836, 214.22174835205078, 214.22179412841797, 214.22196197509766]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=1, batch_per_job=4)","[80.1717243748717, 83.28627449390478, 84.89595309901051, 84.0289533289615, 83.5839789728634, 83.637083777925, 83.6863488140516, 83.56131676700898, 83.67058414593339, 83.96652988507412]","[214.1972532272339, 214.21961975097656, 214.22115325927734, 214.22126007080078, 214.22142791748047, 214.2216567993164, 214.2217025756836, 214.22174835205078, 214.22179412841797, 214.22196197509766]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=4, batch_per_job=1)","[13.828872990095988, 12.610539071029052, 12.713231677887961, 12.66187134385109, 12.670638750074431, 12.614489786094055, 12.158315289067104, 12.280425124103203, 12.518137340899557, 12.298528810963035]","[214.1972532272339, 219.16725826263428, 220.7945785522461, 221.59795761108398, 220.79603576660156, 219.19250106811523, 218.75469589233398, 236.84133529663086, 236.8415756225586, 236.84194564819336]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=4, batch_per_job=2)","[12.361622188938782, 12.386297964956611, 12.436157736927271, 12.6223130710423, 12.273970633978024, 12.40097718895413, 12.59156500804238, 12.48220199206844, 12.502566652139649, 12.620846503879875]","[214.1971616744995, 235.57226181030273, 235.19981002807617, 233.5980987548828, 233.5986671447754, 232.39707565307617, 230.3914680480957, 228.7782440185547, 226.7861099243164, 225.58290100097656]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=4, batch_per_job=4)","[12.524993978906423, 12.615317760966718, 12.659746569115669, 12.659019827842712, 12.778477458981797, 12.583421828923747, 12.817983750952408, 12.479396475944668, 12.84060882893391, 12.536849381169304]","[214.1971616744995, 222.3494873046875, 222.3649787902832, 222.36820220947266, 223.5717315673828, 222.37197875976562, 222.37269973754883, 222.3477554321289, 222.37321853637695, 221.17063522338867]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=8, batch_per_job=1)","[12.212124633835629, 10.334988625021651, 9.897353612817824, 10.057689802953973, 9.868589780060574, 9.582644291920587, 10.306122920941561, 10.168412658153102, 9.12369826901704, 10.27396033401601]","[214.1971616744995, 243.41093730926514, 235.39967346191406, 229.15591049194336, 221.27819442749023, 235.73096466064453, 238.0448989868164, 228.9281120300293, 225.71957397460938, 246.6971778869629]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=8, batch_per_job=2)","[10.128737315069884, 9.820236494066194, 10.145853808848187, 9.686696131015196, 10.86166800512001, 9.629130020970479, 10.03125728899613, 9.1216912670061, 10.092704117996618, 9.36665769200772]","[214.1971616744995, 231.49182891845703, 226.0737648010254, 216.7766876220703, 235.2882194519043, 230.89139938354492, 229.50551223754883, 221.2546272277832, 243.60263442993164, 232.06815338134766]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=8, batch_per_job=4)","[9.629435261944309, 9.06287576816976, 9.76156359212473, 9.901172891026363, 9.298036752035841, 10.096929386025295, 9.72219236381352, 9.93314010091126, 10.15223568212241, 9.94538566796109]","[214.1971616744995, 221.0069351196289, 235.28239822387695, 230.47023010253906, 225.6230354309082, 248.24976348876953, 234.10868072509766, 233.74130249023438, 225.67836380004883, 218.7703514099121]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=16, batch_per_job=1)","[13.102114401059225, 9.55464343400672, 9.99851996009238, 10.199892106000334, 9.763394037028775, 9.965421123895794, 9.915003299014643, 9.895154868951067, 8.901183387963101, 10.199274356011301]","[214.1971616744995, 232.56503868103027, 245.2825803756714, 235.83797359466553, 225.93012523651123, 235.81061458587646, 229.1607847213745, 238.83708095550537, 232.58728313446045, 242.66749477386475]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=16, batch_per_job=2)","[9.91651511611417, 9.744320638012141, 10.302636603126302, 9.793761774199083, 10.447240883018821, 9.809983982006088, 8.742376008071005, 9.312882976140827, 9.90905388793908, 9.764009084086865]","[214.1971616744995, 222.58468627929688, 235.73210525512695, 226.26560974121094, 238.76808547973633, 229.12213134765625, 238.76989364624023, 227.57976531982422, 241.89985275268555, 230.3344497680664]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=16, batch_per_job=4)","[9.871994839049876, 8.97750040399842, 8.659272033022717, 9.923719692975283, 9.819190895883366, 10.091103576123714, 9.176302388776094, 9.12500288314186, 9.706911295186728, 10.242298427037895]","[214.1971616744995, 231.24578857421875, 245.1859474182129, 235.35590362548828, 222.64037704467773, 235.5442657470703, 225.86821365356445, 235.54329681396484, 226.0712547302246, 239.98149871826172]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=32, batch_per_job=1)","[15.313967250986025, 11.086807538988069, 11.638119252864271, 11.017481844173744, 12.048000205075368, 12.495711809955537, 12.307685813168064, 11.896957163931802, 12.394495534012094, 12.626006761100143]","[214.1971616744995, 245.9509916305542, 236.2114963531494, 226.8898696899414, 241.5643310546875, 233.5929822921753, 227.00633907318115, 220.28198528289795, 246.6159381866455, 240.48241806030273]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=32, batch_per_job=2)","[12.521056222030893, 11.173682172084227, 12.375536919105798, 11.180217234883457, 12.070607520174235, 11.774650245904922, 12.067700330168009, 12.548655533930287, 12.411118066171184, 11.716269023017958]","[214.1971616744995, 227.33865356445312, 242.0353889465332, 233.54219436645508, 227.97179889678955, 221.5887851715088, 247.3842430114746, 237.48651123046875, 232.1358346939087, 225.63086986541748]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=32, batch_per_job=4)","[11.238187476061285, 11.408645990071818, 11.363594421185553, 12.129189256113023, 11.753697399049997, 11.738398899091408, 12.538123136851937, 11.482952048070729, 11.560027429834008, 11.490204673958942]","[214.1971616744995, 234.82097625732422, 245.82017135620117, 236.2020492553711, 224.9736442565918, 218.31307125091553, 241.60422039031982, 234.63870334625244, 223.1679573059082, 236.14303588867188]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=64, batch_per_job=1)","[15.420067131053656, 13.773178801173344, 13.53922093193978, 17.184249887010083, 14.836532048881054, 17.421873118029907, 14.634188409196213, 15.463111788965762, 14.053519502980635, 13.865251197945327]","[214.1971616744995, 228.69264602661133, 236.57373237609863, 241.68749237060547, 227.63704109191895, 238.33738708496094, 224.12960815429688, 234.6660614013672, 230.1973533630371, 233.27981185913086]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=64, batch_per_job=2)","[14.3959746430628, 17.011672231834382, 16.294065695954487, 16.575546917971224, 13.808970778947696, 14.45977657311596, 13.634430977050215, 15.287360797869042, 15.643827926833183, 14.460907170781866]","[214.1971616744995, 231.35739135742188, 218.9172477722168, 236.74262237548828, 232.03815460205078, 238.65717887878418, 243.77961921691895, 224.71175575256348, 214.99603080749512, 236.29651641845703]"
+"benchmark(dataset_size=medium, batched=False, n_jobs=64, batch_per_job=4)","[14.079790257150307, 15.261586973210797, 16.608553173020482, 14.337346370099112, 14.052023435942829, 15.21568244509399, 17.222819039830938, 14.844108325894922, 16.11783075495623, 14.915371278999373]","[214.1971616744995, 223.70212364196777, 215.6741123199463, 233.20819282531738, 219.37741470336914, 224.4832878112793, 231.53609561920166, 219.67430114746094, 244.18278694152832, 222.2429962158203]"
diff --git a/benchmarks/utils/__init__.py b/benchmarks/utils/__init__.py
new file mode 100644
index 0000000..23299f9
--- /dev/null
+++ b/benchmarks/utils/__init__.py
@@ -0,0 +1,3 @@
+from .monitor import monitor, parse_func_repr, repr_func
+from ._various import find_result, find_results, choose_file
+from ._argparser import default_parser
diff --git a/benchmarks/utils/_argparser.py b/benchmarks/utils/_argparser.py
new file mode 100644
index 0000000..7ab0e26
--- /dev/null
+++ b/benchmarks/utils/_argparser.py
@@ -0,0 +1,22 @@
+from argparse import ArgumentParser
+
+# Inherit this parser in your file, like so:
+# >>> _parser = ArgumentParser(parents=[default_parser])
+
+default_parser = ArgumentParser(add_help=False)
+
+default_parser.add_argument(
+    "--run",
+    help="Runs the benchmark.",
+    action="store_true",
+)
+default_parser.add_argument(
+    "--plot",
+    help=(
+        "Plots the results. If '--run' is specified, plots those. "
+        "Otherwise, searches for result files in the results' directory. "
+        "If it finds none, exits. If it finds only one, displays it. "
+        "If it finds multiple, prompts the user to choose one. "
+    ),
+    action="store_true",
+)
diff --git a/benchmarks/utils/_various.py b/benchmarks/utils/_various.py
new file mode 100644
index 0000000..a55979f
--- /dev/null
+++ b/benchmarks/utils/_various.py
@@ -0,0 +1,58 @@
+import pandas as pd
+
+from typing import List
+from pathlib import Path
+
+
+def find_result(bench_name: str) -> Path:
+    return choose_file(find_results(bench_name))
+
+
+def find_results(bench_name: str) -> List[Path]:
+    """
+    Returns the list of results in the results' directory.
+    """
+    results_dir = Path(__file__).parent.parent / "results"
+    return [
+        file
+        for file in results_dir.iterdir()
+        if file.stem.startswith(bench_name) and file.suffix == ".csv"
+    ]
+
+
+def choose_file(results: List[Path]) -> Path:
+    """
+    Given a list of files, chooses one based on these rules:
+    - If there are no files to choose from, exit the program
+    - If there's only one, return this one
+    - If there are multiple, prompt the user to choose one
+    """
+    if len(results) == 0:
+        print("No results file to choose from, exiting...")
+        exit()
+    elif len(results) == 1:
+        return results[0]
+    else:
+        for i, file in enumerate(results):
+            # Read the result file to get its dimensions
+            df = pd.read_csv(file)
+            # Try and get the number of repeats.
+            if "time" in df.columns:
+                repeat = len(eval(df["time"].iloc[0]))
+            elif "memory" in df.columns:
+                repeat = len(eval(df["memory"].iloc[0]))
+            else:
+                repeat = 1
+
+            bench_name, date, commit = file.stem.split("-")
+            print(
+                f"{i + 1}) "
+                f"{date[:4]}-{date[4:6]}-{date[6:]} - "
+                f"{df.shape[0]}x{repeat} experiments "
+                f"(git {commit[:8]})"
+            )
+        choice = input("Choose the result to display: ")
+        if not choice.isnumeric() or (int(choice) - 1) not in range(len(results)):
+            print(f"Invalid choice {choice!r}, exiting.")
+            exit()
+        return results[int(choice) - 1]
diff --git a/benchmarks/utils/monitor.py b/benchmarks/utils/monitor.py
new file mode 100644
index 0000000..a73b8a2
--- /dev/null
+++ b/benchmarks/utils/monitor.py
@@ -0,0 +1,213 @@
+import tracemalloc
+import pandas as pd
+
+from typing import Callable, Collection, Any, Optional, Tuple, Dict, List, Union
+from time import perf_counter
+from itertools import product as _product
+from collections import defaultdict
+from datetime import datetime
+from warnings import warn
+from pathlib import Path
+from git import Repo
+
+
+def repr_func(f: Callable, args: tuple, kwargs: dict) -> str:
+    """
+    Takes a function (f) and its arguments (args, kwargs), and
+    returns it represented as "f(*args, **kwargs)".
+    For example, with ``f=do_smth``, ``args=(10, 5)`` and
+    ``kwargs={"keyboard": "qwerty"}``,
+    returns "do_smth(10, 5, keyboard=qwerty)".
+    Can be parsed with the function `parse_func_repr` below.
+    """
+    str_args = ", ".join(args)
+    str_kwargs = ", ".join(f"{k}={v}" for k, v in kwargs.items())
+    return f"{f.__name__}({', '.join(st for st in [str_args, str_kwargs] if st)})"
+
+
+def parse_func_repr(representation: str) -> Tuple[str, tuple, dict]:
+    """
+    Takes the representation of a function and its arguments created by
+    `repr_func`, and returns the function name,
+    the positional arguments as a tuple,
+    and the keyword arguments as a dictionary.
+    """
+    func_name, args_repr = representation[:-1].split("(", 1)
+    args = []
+    kwargs = {}
+    for arg in args_repr.split(", "):
+        if "=" in arg:
+            keyword, argument = arg.split("=", 1)
+            if keyword.isidentifier():
+                kwargs.update({keyword: argument})
+                continue
+        args.append(arg)
+    return func_name, tuple(args), kwargs
+
+
+def monitor(
+    memory: bool,
+    time: bool,
+    parametrize: Union[Collection[Collection[Any]], Dict[str, Collection[Any]]],
+    repeat: int = 1,
+    save_as: Optional[str] = None,
+) -> Callable[..., Callable[..., pd.DataFrame]]:
+    """Decorator used to monitor the execution of a function.
+
+    The decorated function should return nothing (even if it does, nothing will
+    be passed through by this decorator).
+    Executions are sequential, so it's usually pretty long to run!
+
+    Parameters
+    ----------
+
+    memory : bool
+        Whether the RAM usage should be monitored throughout the function execution.
+        Note: monitored in the main thread.
+    time : bool
+        Whether the time the function took to execute should be measured.
+        Note: if `memory` is also set, consider that as the memory profiler runs
+        in the main thread, the timings will be different from an execution
+        without the memory monitoring.
+    parametrize : a collection of a collection of parameters
+        Specifies the parameter matrix to be used on the function.
+        These can either be passed as positional arguments (e.g., list of list
+        of parameters), or as keyword arguments (e.g., dictionary of list of
+        parameters).
+        Note: when specified, ignores the parameters passed to the function.
+    repeat : int
+        How many times we want to repeat the execution of the function for more
+        representative time and memory extracts.
+    save_as : str, optional
+        Can be specified as a benchmark name for the results to be automatically
+        saved on disk.
+        E.g. "supervectorizer_tuning"
+        If None (the default), results are not saved on disk, only returned.
+
+    Returns
+    -------
+
+    Callable[..., Callable[..., pd.DataFrame]]
+        A double-nested callable that returns a DataFrame of the results.
+
+    """
+    if not memory and not time:
+        warn(
+            "Parameters 'memory' and 'time' are both set to False ; "
+            "there is therefore nothing to monitor for, returning empty. ",
+            stacklevel=2,
+        )
+        return lambda *_, **__: lambda *_, **__: pd.DataFrame()
+
+    def decorator(func: Callable[[Any], None]):
+        """Only catches the decorated function."""
+
+        def wrapper(*_, **__) -> pd.DataFrame:
+            """
+            Parameters
+            ----------
+
+            *_ : Tuple[Any]
+                Arguments passed by the function call. Ignored.
+                Use `parametrize` instead.
+            **__ : Dict[str, Any]
+                Keyword arguments passed by the function call. Ignored.
+                Use `parametrize` instead.
+
+            """
+
+            def product(
+                iterables: Union[
+                    Collection[Collection[Any]], Dict[str, Collection[Any]]
+                ],
+            ):
+                """
+                Wrapper for ``itertools.product`` in order to better
+                handle dictionaries.
+                """
+                if isinstance(iterables, dict):
+                    for pair in _product(*iterables.values()):
+                        yield (), dict(zip(iterables.keys(), pair))
+                else:
+                    for params in _product(*iterables):
+                        return params, {}
+
+            def exec_func(*args, **kwargs) -> Dict[str, List[float]]:
+                """
+                Parameters
+                ----------
+
+                *args : Tuple[Any]
+                    Arguments to pass to the function.
+                **kwargs : Dict[str, Any]
+                    Keyword arguments to pass to the function.
+
+                Returns
+                -------
+                Dict[str, List[float]]
+                    a mapping of monitored resource name to a list of values:
+                    - "time": how long the execution took, in seconds.
+                      Only present if ``time=True``.
+                    - "memory": The number of MB of memory used.
+                      Only present if ``memory=True``.
+                    The size of these lists is equal to `repeat`.
+
+                """
+                _monitored = defaultdict(lambda: [])
+
+                # Global initialization of monitored values
+                if memory:
+                    tracemalloc.start()
+
+                for n in range(repeat):
+                    # Initialize loop monitored values
+                    if memory:
+                        tracemalloc.reset_peak()
+                    if time:
+                        t0 = perf_counter()
+
+                    func(*args, **kwargs)
+
+                    # Collect and store loop monitored values
+                    if time:
+                        _monitored["time"].append(perf_counter() - t0)
+                    if memory:
+                        _, peak = tracemalloc.get_traced_memory()
+                        _monitored["memory"].append(peak / (1024**2))
+
+                # Global cleanup of monitored values
+                if memory:
+                    tracemalloc.stop()
+
+                return dict(_monitored)
+
+            def to_df(
+                results: Dict[str, Dict[str, List[float]]],
+            ) -> pd.DataFrame:
+                """
+                Converts the result of the benchmark to a pandas DataFrame.
+                """
+                index = list(results.keys())
+                columns = results[index[0]].keys()  # Use the first one as sample
+                data = [[results[idx][col] for col in columns] for idx in index]
+                return pd.DataFrame(data, index=index, columns=columns)
+
+            results = {
+                repr_func(func, args, kwargs): exec_func(*args, **kwargs)
+                for args, kwargs in product(parametrize)
+            }
+            df = to_df(results)
+
+            if save_as is not None:
+                save_dir = Path(__file__).parent.parent / "results"
+                save_dir.mkdir(exist_ok=True)
+                now = datetime.now()
+                commit = Repo(Path(__file__).parent.parent.parent).active_branch.commit
+                file = f"{save_as}-{now.year}{now.month}{now.day}-{commit}.csv"
+                df.to_csv(save_dir / file, index_label="call")
+
+            return df
+
+        return wrapper
+
+    return decorator
diff --git a/dirty_cat/_datetime_encoder.py b/dirty_cat/_datetime_encoder.py
index 035debe..29e2ff5 100644
--- a/dirty_cat/_datetime_encoder.py
+++ b/dirty_cat/_datetime_encoder.py
@@ -186,7 +186,7 @@ class DatetimeEncoder(BaseEstimator, TransformerMixin):
 
         return self
 
-    def transform(self, X, y=None) -> np.array:
+    def transform(self, X, y=None) -> np.ndarray:
         """Transform X by replacing each datetime column with corresponding numerical features.
 
         Parameters
@@ -198,7 +198,7 @@ class DatetimeEncoder(BaseEstimator, TransformerMixin):
 
         Returns
         -------
-        np.array, shape (n_samples, n_features_out_)
+        np.ndarray, shape (n_samples, n_features_out_)
             Transformed input.
         """
         check_is_fitted(
diff --git a/dirty_cat/_gap_encoder.py b/dirty_cat/_gap_encoder.py
index fe2e4fc..61302f2 100644
--- a/dirty_cat/_gap_encoder.py
+++ b/dirty_cat/_gap_encoder.py
@@ -46,7 +46,7 @@ class GapEncoderColumn(BaseEstimator, TransformerMixin):
     """See GapEncoder's docstring."""
 
     rho_: float
-    H_dict_: Dict[np.array, np.array]
+    H_dict_: Dict[np.ndarray, np.ndarray]
 
     def __init__(
         self,
@@ -88,7 +88,7 @@ class GapEncoderColumn(BaseEstimator, TransformerMixin):
         self.rescale_W = rescale_W
         self.max_iter_e_step = max_iter_e_step
 
-    def _init_vars(self, X) -> Tuple[np.array, np.array, np.array]:
+    def _init_vars(self, X) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
         """
         Build the bag-of-n-grams representation V of X and initialize
         the topics W.
diff --git a/dirty_cat/_minhash_encoder.py b/dirty_cat/_minhash_encoder.py
index fee46bf..9518c02 100644
--- a/dirty_cat/_minhash_encoder.py
+++ b/dirty_cat/_minhash_encoder.py
@@ -16,11 +16,12 @@ probability of having same encoding values. These encodings thus capture
 morphological similarities between strings.
 """
 
-from typing import Dict, List, Literal, Tuple
+from typing import Callable, Collection, Dict, List, Literal, Tuple
 
 import numpy as np
+from joblib import Parallel, delayed, effective_n_jobs
 from sklearn.base import BaseEstimator, TransformerMixin
-from sklearn.utils import murmurhash3_32
+from sklearn.utils import gen_even_slices, murmurhash3_32
 
 from ._fast_hash import ngram_min_hash
 from ._string_distances import get_unique_ngrams
@@ -29,6 +30,10 @@ from ._utils import LRUDict, check_input
 NoneType = type(None)
 
 
+# Ignore lines too long, as links in the docstring cannot be cut.
+# flake8: noqa: E501
+
+
 class MinHashEncoder(BaseEstimator, TransformerMixin):
     """
     Encode string categorical features as a numeric array, minhash method
@@ -52,6 +57,14 @@ class MinHashEncoder(BaseEstimator, TransformerMixin):
     handle_missing : typing.Literal["error", "zero_impute"], default=zero_impute
         Whether to raise an error or encode missing values (NaN) with
         vectors filled with zeros.
+    n_jobs : int, default=None
+        The number of jobs to run in parallel.
+        The hash computations for all unique elements are parallelized.
+        None means 1 unless in a
+        `joblib.parallel_backend context <https://joblib.readthedocs.io/en/latest/parallel.html>`_.
+        -1 means using all processors.
+        See `Scikit-learn Glossary <https://scikit-learn.org/stable/glossary.html#term-n_jobs>`_
+        for more details.
 
     Attributes
     ----------
@@ -100,12 +113,14 @@ class MinHashEncoder(BaseEstimator, TransformerMixin):
         hashing: Literal["fast", "murmur"] = "fast",
         minmax_hash: bool = False,
         handle_missing: Literal["error", "zero_impute"] = "zero_impute",
+        n_jobs: int = None,
     ):
         self.ngram_range = ngram_range
         self.n_components = n_components
         self.hashing = hashing
         self.minmax_hash = minmax_hash
         self.handle_missing = handle_missing
+        self.n_jobs = n_jobs
 
     def _more_tags(self) -> Dict[str, List[str]]:
         """
@@ -121,12 +136,6 @@ class MinHashEncoder(BaseEstimator, TransformerMixin):
         ----------
         string : str
             The string to encode.
-        n_components : int
-            The number of dimension of encoded string.
-        ngram_range : typing.Tuple[int, int]
-            The lower and upper boundaries of the range of n-values for
-            different n-grams to be extracted.
-            All values of n such that min_n <= n <= max_n.
 
         Returns
         -------
@@ -177,6 +186,36 @@ class MinHashEncoder(BaseEstimator, TransformerMixin):
                 ]
             )
 
+    def _compute_hash_batched(
+        self, batch: Collection[str], hash_func: Callable[[str], np.ndarray]
+    ):
+        """Function called to compute the hashes of a batch of strings.
+
+        Check if the string is in the hash dictionary, if not, compute the hash using
+        the specified hashing function and add it to the dictionary.
+
+        Parameters
+        ----------
+        batch : iterable of str
+            The batch of strings to encode.
+        hash_func : callable
+            Hashing function to use on the string.
+
+        Returns
+        -------
+        np.array of shape (n_samples, n_components)
+            The encoded strings, using specified encoding scheme.
+        """
+        res = np.zeros((len(batch), self.n_components))
+        for i, string in enumerate(batch):
+            if string not in self.hash_dict_:
+                if string == "NAN":  # true if x is a missing value
+                    self.hash_dict_[string] = np.zeros(self.n_components)
+                else:
+                    self.hash_dict_[string] = hash_func(string)
+            res[i] = self.hash_dict_[string]
+        return res
+
     def fit(self, X, y=None) -> "MinHashEncoder":
         """
         Fit the MinHashEncoder to X. In practice, just initializes a dictionary
@@ -185,7 +224,7 @@ class MinHashEncoder(BaseEstimator, TransformerMixin):
         Parameters
         ----------
         X : array-like, shape (n_samples, ) or (n_samples, 1)
-            The string data to encode.
+            The string data to encode. Only here for compatibility.
         y : None
             Unused, only here for compatibility.
 
@@ -204,7 +243,6 @@ class MinHashEncoder(BaseEstimator, TransformerMixin):
                 f"Got handle_missing={self.handle_missing!r}, but expected "
                 'any of {"error", "zero_impute"}. '
             )
-        self.count = 0
         self.hash_dict_ = LRUDict(capacity=self._capacity)
         return self
 
@@ -224,59 +262,67 @@ class MinHashEncoder(BaseEstimator, TransformerMixin):
         """
         X = check_input(X)
         if self.minmax_hash:
-            assert (
-                self.n_components % 2 == 0
-            ), "n_components should be even when minmax_hash=True. "
+            if self.n_components % 2 != 0:
+                raise ValueError(
+                    "n_components should be even when using"
+                    f"minmax_hash encoding, got {self.n_components}"
+                )
         if self.hashing == "murmur":
-            assert (
-                not self.minmax_hash
-            ), 'minmax_hash is not implemented with hashing="murmur". '
+            if self.minmax_hash:
+                raise ValueError(
+                    "minmax_hash encoding is not supported"
+                    "with the murmur hashing function"
+                )
+        if self.handle_missing not in ["error", "zero_impute"]:
+            raise ValueError(
+                "handle_missing should be either "
+                f"'error' or 'zero_impute', got {self.handle_missing!r}"
+            )
 
-        # TODO: Parallelize
-        is_nan_idx = False
+        # Handle missing values
+        missing_mask = (
+            ~(X == X)  # Find np.nan
+            | (X == None)  # Find None. Note: `X is None` doesn't work.
+            | (X == "")  # Find empty strings
+        )
+
+        if missing_mask.any():  # contains at least one missing value
+            if self.handle_missing == "error":
+                raise ValueError(
+                    "Found missing values in input data; set "
+                    "handle_missing='zero_impute' to encode with missing values"
+                )
+            elif self.handle_missing == "zero_impute":
+                # NANs will be replaced by zeroes in _compute_hash
+                X[missing_mask] = "NAN"
 
         if self.hashing == "fast":
-            X_out = np.zeros((len(X[:]), self.n_components * X.shape[1]))
-            counter = self.n_components
-            for k in range(X.shape[1]):
-                X_in = X[:, k].reshape(-1)
-                for i, x in enumerate(X_in):
-                    if (
-                        isinstance(x, float)  # true if x is a missing value
-                        or x is None
-                    ):
-                        is_nan_idx = True
-                    elif len(x) == 0:
-                        is_nan_idx = True
-                    elif x not in self.hash_dict_:
-                        X_out[i, k * self.n_components : counter] = self.hash_dict_[
-                            x
-                        ] = self._get_fast_hash(x)
-                    else:
-                        X_out[i, k * self.n_components : counter] = self.hash_dict_[x]
-                counter += self.n_components
-        if self.hashing == "murmur":
-            X_out = np.zeros((len(X[:]), self.n_components * X.shape[1]))
-            counter = self.n_components
-            for k in range(X.shape[1]):
-                X_in = X[:, k].reshape(-1)
-                for i, x in enumerate(X_in):
-                    if isinstance(x, float) or x is None:
-                        is_nan_idx = True
-                    elif len(x) == 0:
-                        is_nan_idx = True
-                    elif x not in self.hash_dict_:
-                        X_out[i, k * self.n_components : counter] = self.hash_dict_[
-                            x
-                        ] = self._get_murmur_hash(x)
-                    else:
-                        X_out[i, k * self.n_components : counter] = self.hash_dict_[x]
-                counter += self.n_components
-
-        if is_nan_idx and self.handle_missing == "error":
+            hash_func = self._get_fast_hash
+        elif self.hashing == "murmur":
+            hash_func = self._get_murmur_hash
+        else:
             raise ValueError(
-                "Found missing values in input data; set "
-                "handle_missing='zero_impute' "
-                "to encode with missing values. "
+                "Hashing function should be either 'fast' or 'murmur', "
+                f"got {self.hashing!r}"
             )
-        return X_out
+
+        # Compute the hashes for unique values
+        unique_x, indices_x = np.unique(X, return_inverse=True)
+        n_jobs = effective_n_jobs(self.n_jobs)
+
+        # Compute the hashes in parallel on n_jobs batches
+        unique_x_trans = Parallel(n_jobs=n_jobs)(
+            delayed(self._compute_hash_batched)(
+                unique_x[idx_slice],
+                hash_func,
+            )
+            for idx_slice in gen_even_slices(len(unique_x), 
+            n_jobs)
+        )
+
+        # Match the hashes of the unique value to the original values
+        X_out = np.concatenate(unique_x_trans)[indices_x].reshape(
+            len(X), X.shape[1] * self.n_components
+        )
+
+        return X_out.astype(np.float64)  # The output is an int32 before conversion
diff --git a/dirty_cat/_similarity_encoder.py b/dirty_cat/_similarity_encoder.py
index 9b1a99b..bbc3927 100644
--- a/dirty_cat/_similarity_encoder.py
+++ b/dirty_cat/_similarity_encoder.py
@@ -33,12 +33,12 @@ from ._utils import parse_version
 
 
 def _ngram_similarity_one_sample_inplace(
-    x_count_vector: np.array,
-    vocabulary_count_matrix: np.array,
+    x_count_vector: np.ndarray,
+    vocabulary_count_matrix: np.ndarray,
     str_x: str,
-    vocabulary_ngram_counts: np.array,
+    vocabulary_ngram_counts: np.ndarray,
     se_dict: dict,
-    unq_X: np.array,
+    unq_X: np.ndarray,
     i: int,
     ngram_range: Tuple[int, int],
 ) -> None:
@@ -47,17 +47,17 @@ def _ngram_similarity_one_sample_inplace(
 
     Parameters
     ----------
-    x_count_vector : np.array
+    x_count_vector : np.ndarray
         Count vector of the sample based on the ngrams of the vocabulary
-    vocabulary_count_matrix : np.array
+    vocabulary_count_matrix : np.ndarray
         Count vector of the vocabulary based on its ngrams
     str_x: str
         The actual sample string
-    vocabulary_ngram_counts : np.array
+    vocabulary_ngram_counts : np.ndarray
         Number of ngrams for each unique element of the vocabulary
     se_dict : dict
         Dictionary containing the similarities for each x in unq_X
-    unq_X : np.array
+    unq_X : np.ndarray
         The arrays of all unique samples
     i : str
         The index of x_count_vector in the csr count matrix
@@ -86,7 +86,7 @@ def ngram_similarity(
     ngram_range: Tuple[int, int],
     hashing_dim: int,
     dtype: type = np.float64,
-) -> np.array:
+) -> np.ndarray:
     """
     Similarity encoding for dirty categorical variables:
     Given two arrays of strings, returns the similarity encoding matrix
@@ -139,7 +139,7 @@ def ngram_similarity(
     return np.nan_to_num(out, copy=False)
 
 
-def get_prototype_frequencies(prototypes: np.array) -> np.array:
+def get_prototype_frequencies(prototypes: np.ndarray) -> np.array:
     """
     Computes the frequencies of the values contained in prototypes
     Reverse sorts the array by the frequency
diff --git a/dirty_cat/_target_encoder.py b/dirty_cat/_target_encoder.py
index f829661..a82ddfa 100644
--- a/dirty_cat/_target_encoder.py
+++ b/dirty_cat/_target_encoder.py
@@ -55,7 +55,7 @@ class TargetEncoder(BaseEstimator, TransformerMixin):
     ----------
     n_features_in_: int
         Number of features in the data seen during fit.
-    categories_ : typing.List[np.array]
+    categories_ : typing.List[np.ndarray]
         The categories of each feature determined during fitting
         (in order corresponding with output of ``transform``).
 
@@ -69,7 +69,7 @@ class TargetEncoder(BaseEstimator, TransformerMixin):
     --------
     >>> enc = TargetEncoder(handle_unknown='ignore')
     >>> X = [['male'], ['Male'], ['Female'], ['male'], ['Female']]
-    >>> y = np.array([1, 2, 3, 4, 5])
+    >>> y = np.ndarray([1, 2, 3, 4, 5])
 
     >>> enc.fit(X, y)
     TargetEncoder(handle_unknown='ignore')
@@ -94,12 +94,12 @@ class TargetEncoder(BaseEstimator, TransformerMixin):
 
     n_features_in_: int
     _label_encoders_: List[LabelEncoder]
-    categories_: List[np.array]
+    categories_: List[np.ndarray]
     n_: int
 
     def __init__(
         self,
-        categories: Union[Literal["auto"], List[Union[List[str], np.array]]] = "auto",
+        categories: Union[Literal["auto"], List[Union[List[str], np.ndarray]]] = "auto",
         clf_type: Literal["regression", "binary-clf", "multiclass-clf"] = "binary-clf",
         dtype: type = np.float64,
         handle_unknown: Literal["error", "ignore"] = "error",
@@ -211,7 +211,7 @@ class TargetEncoder(BaseEstimator, TransformerMixin):
         self.k_ = {j: len(self.counter_[j]) for j in self.counter_}
         return self
 
-    def transform(self, X) -> np.array:
+    def transform(self, X) -> np.ndarray:
         """
         Transform X using the specified encoding scheme.
 
@@ -222,7 +222,7 @@ class TargetEncoder(BaseEstimator, TransformerMixin):
 
         Returns
         -------
-        2-d np.array
+        2-d np.ndarray
             Transformed input.
         """
         check_is_fitted(self, attributes=["n_features_in_"])
diff --git a/dirty_cat/_utils.py b/dirty_cat/_utils.py
index 0cd2aee..d578303 100644
--- a/dirty_cat/_utils.py
+++ b/dirty_cat/_utils.py
@@ -41,7 +41,7 @@ class LRUDict:
         return key in self.cache
 
 
-def check_input(X) -> np.array:
+def check_input(X) -> np.ndarray:
     """
     Check input with sklearn standards.
     Also converts X to a numpy array if not already.


Test patch:
diff --git a/dirty_cat/tests/test_minhash_encoder.py b/dirty_cat/tests/test_minhash_encoder.py
index a416046..e88d380 100644
--- a/dirty_cat/tests/test_minhash_encoder.py
+++ b/dirty_cat/tests/test_minhash_encoder.py
@@ -1,17 +1,21 @@
 import random
 from string import ascii_lowercase
 
+import joblib
 import numpy as np
 import pandas as pd
 import pytest
+from sklearn.utils._testing import assert_array_equal, skip_if_no_parallel
 
 from dirty_cat import MinHashEncoder
 
+from .utils import generate_data
+
 
 @pytest.mark.parametrize(
     "hashing, minmax_hash", [("fast", True), ("fast", False), ("murmur", False)]
 )
-def test_MinHashEncoder(hashing, minmax_hash) -> None:
+def test_minhash_encoder(hashing, minmax_hash) -> None:
     X = np.array(["al ice", "b ob", "bob and alice", "alice and bob"])[:, None]
     # Test output shape
     encoder = MinHashEncoder(n_components=2, hashing=hashing)
@@ -76,19 +80,7 @@ def test_input_type() -> None:
     "hashing, minmax_hash", [("fast", True), ("fast", False), ("murmur", False)]
 )
 def test_encoder_params(hashing, minmax_hash) -> None:
-    MAX_LIMIT = 255  # extended ASCII Character set
-    i = 0
-    str_list = []
-    for i in range(100):
-        random_string = "aa"
-        for _ in range(100):
-            random_integer = random.randint(0, MAX_LIMIT)
-            random_string += chr(random_integer)
-            if random_integer < 50:
-                random_string += "  "
-        i += 1
-        str_list += [random_string]
-    X = np.array(str_list).reshape(100, 1)
+    X = generate_data(n_samples=20)
     enc = MinHashEncoder(
         n_components=50, hashing=hashing, minmax_hash=minmax_hash, ngram_range=(3, 3)
     )
@@ -169,3 +161,84 @@ def test_cache_overflow() -> None:
     y = encoder.fit_transform(raw_data)
 
     assert len(y[y == -1.0]) == 0
+
+
+@skip_if_no_parallel
+def test_parallelism():
+    # Test that parallelism works
+    encoder = MinHashEncoder(n_components=3, n_jobs=1)
+    X = np.array(["a", "b", "c", "d", "e", "f", "g", "h"])[:, None]
+    y = encoder.fit_transform(X)
+    for n_jobs in [None, 2, -1]:
+        encoder = MinHashEncoder(n_components=3, n_jobs=n_jobs)
+        y_parallel = encoder.fit_transform(X)
+        assert_array_equal(y, y_parallel)
+
+    # Test with threading backend
+    encoder = MinHashEncoder(n_components=3, n_jobs=2)
+    with joblib.parallel_backend("threading"):
+        y_threading = encoder.fit_transform(X)
+    assert_array_equal(y, y_threading)
+    assert encoder.n_jobs == 2
+
+
+DEFAULT_JOBLIB_BACKEND = joblib.parallel.get_active_backend()[0].__class__
+
+
+class DummyBackend(DEFAULT_JOBLIB_BACKEND):  # type: ignore
+    """
+    A dummy backend used to check that specifying a backend works
+    in MinHashEncoder.
+    The `count` attribute is used to check that the backend is used.
+    Copied from https://github.com/scikit-learn/scikit-learn/blob/36958fb240fbe435673a9e3c52e769f01f36bec0/sklearn/ensemble/tests/test_forest.py  # noqa
+    """
+    def __init__(self, *args, **kwargs):
+        self.count = 0
+        super().__init__(*args, **kwargs)
+
+    def start_call(self):
+        self.count += 1
+        return super().start_call()
+
+
+joblib.register_parallel_backend("testing", DummyBackend)
+
+
+@skip_if_no_parallel
+def test_backend_respected():
+    """
+    Test that the joblib backend is used.
+    Copied from https://github.com/scikit-learn/scikit-learn/blob/36958fb240fbe435673a9e3c52e769f01f36bec0/sklearn/ensemble/tests/test_forest.py  # noqa
+    """
+    # Test that parallelism works
+    encoder = MinHashEncoder(n_components=3, n_jobs=2)
+    X = np.array(["a", "b", "c", "d", "e", "f", "g", "h"])[:, None]
+
+    with joblib.parallel_backend("testing") as (ba, n_jobs):
+        encoder.fit_transform(X)
+        
+    assert ba.count > 0
+
+
+def test_correct_arguments():
+    # Test that the correct arguments are passed to the hashing function
+    X = np.array(["a", "b", "c", "d", "e", "f", "g", "h"])[:, None]
+    # Write an incorrect value for the `hashing` argument
+    with pytest.raises(ValueError, match=r"expected any of"):
+        encoder = MinHashEncoder(n_components=3, hashing="incorrect")
+        encoder.fit_transform(X)
+
+    # Write an incorrect value for the `handle_missing` argument
+    with pytest.raises(ValueError, match=r"expected any of"):
+        encoder = MinHashEncoder(n_components=3, handle_missing="incorrect")
+        encoder.fit_transform(X)
+
+    # Use minmax_hash with murmur hashing
+    with pytest.raises(ValueError, match=r"minmax_hash encoding is not supported"):
+        encoder = MinHashEncoder(n_components=2, minmax_hash=True, hashing="murmur")
+        encoder.fit_transform(X)
+
+    # Use minmax_hash with an odd number of components
+    with pytest.raises(ValueError, match=r"n_components should be even"):
+        encoder = MinHashEncoder(n_components=3, minmax_hash=True)
+        encoder.fit_transform(X)
diff --git a/dirty_cat/tests/test_super_vectorizer.py b/dirty_cat/tests/test_super_vectorizer.py
index 75bbf05..ce600c5 100644
--- a/dirty_cat/tests/test_super_vectorizer.py
+++ b/dirty_cat/tests/test_super_vectorizer.py
@@ -77,7 +77,7 @@ def _get_dirty_dataframe() -> pd.DataFrame:
     )
 
 
-def _get_numpy_array() -> np.array:
+def _get_numpy_array() -> np.ndarray:
     return np.array(
         [
             ["15", "56", pd.NA, "12", ""],
diff --git a/dirty_cat/tests/utils.py b/dirty_cat/tests/utils.py
index 77ca446..baf55bd 100644
--- a/dirty_cat/tests/utils.py
+++ b/dirty_cat/tests/utils.py
@@ -11,13 +11,12 @@ def generate_data(
 ) -> np.ndarray:
     if random_state is not None:
         random.seed(random_state)
-
     MAX_LIMIT = 255  # extended ASCII Character set
     str_list = []
     for i in range(n_samples):
         random_string = "category "
         for _ in range(100):
-            random_integer = random.randint(0, MAX_LIMIT)
+            random_integer = random.randint(1, MAX_LIMIT)
             random_string += chr(random_integer)
             if random_integer < 50:
                 random_string += "  "