diff --git a/README.md b/README.md
index dc7c5dfaa..57b845edb 100644
--- a/README.md
+++ b/README.md
@@ -20,7 +20,11 @@ POT provides the following generic OT solvers (links to examples):
 
 * [OT Network Simplex solver](https://pythonot.github.io/auto_examples/plot_OT_1D.html) for the linear program/ Earth Movers Distance [1] .
 * [Conditional gradient](https://pythonot.github.io/auto_examples/plot_optim_OTreg.html) [6] and [Generalized conditional gradient](https://pythonot.github.io/auto_examples/plot_optim_OTreg.html) for regularized OT [7].
-* Entropic regularization OT solver with [Sinkhorn Knopp Algorithm](https://pythonot.github.io/auto_examples/plot_OT_1D.html) [2] , stabilized version [9] [10] [34], greedy Sinkhorn [22] and [Screening Sinkhorn [26] ](https://pythonot.github.io/auto_examples/plot_screenkhorn_1D.html).
+* Entropic regularization OT solver with [Sinkhorn Knopp
+  Algorithm](https://pythonot.github.io/auto_examples/plot_OT_1D.html) [2] ,
+  stabilized version [9] [10] [34], lazy CPU/GPU solver from geomloss [60] [61], greedy Sinkhorn [22] and [Screening
+  Sinkhorn [26]
+  ](https://pythonot.github.io/auto_examples/plot_screenkhorn_1D.html).
 * Bregman projections for [Wasserstein barycenter](https://pythonot.github.io/auto_examples/barycenters/plot_barycenter_lp_vs_entropic.html) [3], [convolutional barycenter](https://pythonot.github.io/auto_examples/barycenters/plot_convolutional_barycenter.html) [21]  and unmixing [4].
 * Sinkhorn divergence [23] and entropic regularization OT from empirical data.
 * Debiased Sinkhorn barycenters [Sinkhorn divergence barycenter](https://pythonot.github.io/auto_examples/barycenters/plot_debiased_barycenter.html) [37]
diff --git a/ot/solvers.py b/ot/solvers.py
index b65f1149d..a41762a5c 100644
--- a/ot/solvers.py
+++ b/ot/solvers.py
@@ -1001,6 +1001,26 @@ def solve_sample(X_a, X_b, a=None, b=None, metric='sqeuclidean', reg=None, reg_t
         # lazy OT plan
         lazy_plan = res.lazy_plan
 
+    We also have a very efficient solver with compiled CPU/CUDA code using
+    geomloss/PyKeOps that can be used with the following code:
+
+    .. code-block:: python
+
+        # automatic solver
+        res = ot.solve_sample(xa, xb, a, b, reg=1.0, method='geomloss')
+
+        # force O(n) memory efficient solver
+        res = ot.solve_sample(xa, xb, a, b, reg=1.0, method='geomloss_online')
+
+        # force pre-computed cost matrix
+        res = ot.solve_sample(xa, xb, a, b, reg=1.0, method='geomloss_tensorized')
+
+        # use multiscale solver
+        res = ot.solve_sample(xa, xb, a, b, reg=1.0, method='geomloss_multiscale')
+
+        # One can play with speed (small scaling factor) and precision (scaling close to 1)
+        res = ot.solve_sample(xa, xb, a, b, reg=1.0, method='geomloss', scaling=0.5)
+
     - **Quadratic regularized OT [17]** (when ``reg!=None`` and ``reg_type="L2"``):
 
     .. math::