Merge pull request #50 from rflamary/smooth_ot

Smooth and Sparse OT

Author: rflamary

README.md

@@ -15,7 +15,8 @@ It provides the following solvers:
 
 * OT Network Flow solver for the linear program/ Earth Movers Distance [1].
 * Entropic regularization OT solver with Sinkhorn Knopp Algorithm [2] and stabilized version [9][10] with optional GPU implementation (requires cudamat).
-* Non regularized Wasserstein barycenters [16] with LP solver.
+* Smooth optimal transport solvers (dual and semi-dual) for KL and squared L2 regularizations [17].
+* Non regularized Wasserstein barycenters [16] with LP solver (only small scale).
 * Bregman projections for Wasserstein barycenter [3] and unmixing [4].
 * Optimal transport for domain adaptation with group lasso regularization [5]
 * Conditional gradient [6] and Generalized conditional gradient for regularized OT [7].
@@ -29,10 +30,11 @@ Some demonstrations (both in Python and Jupyter Notebook format) are available i
 
 If you use this toolbox in your research and find it useful, please cite POT using the following bibtex reference:
 ```
-@article{flamary2017pot,
-  title={POT Python Optimal Transport library},
-  author={Flamary, R{\'e}mi and Courty, Nicolas},
-  year={2017}
+@misc{flamary2017pot,
+title={POT Python Optimal Transport library},
+author={Flamary, R{'e}mi and Courty, Nicolas},
+url={https://github.com/rflamary/POT},
+year={2017}
 }
 ```
 
@@ -215,3 +217,5 @@ You can also post bug reports and feature requests in Github issues. Make sure t
 [15] Peyré, G., & Cuturi, M. (2018). [Computational Optimal Transport](https://arxiv.org/pdf/1803.00567.pdf) .
 
 [16] Agueh, M., & Carlier, G. (2011). [Barycenters in the Wasserstein space](https://hal.archives-ouvertes.fr/hal-00637399/document). SIAM Journal on Mathematical Analysis, 43(2), 904-924.
+
+[17] Blondel, M., Seguy, V., & Rolet, A. (2018). [Smooth and Sparse Optimal Transport](https://arxiv.org/abs/1710.06276). Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics (AISTATS).

docs/source/all.rst

@@ -19,6 +19,11 @@ ot.bregman
 
 .. automodule:: ot.bregman
    :members:
+   
+ot.smooth
+-----
+.. automodule:: ot.smooth
+   :members:
 
 ot.gromov
 ----------

docs/source/readme.rst

@@ -14,7 +14,11 @@ It provides the following solvers:
    [1].
 -  Entropic regularization OT solver with Sinkhorn Knopp Algorithm [2]
    and stabilized version [9][10] with optional GPU implementation
-   (required cudamat).
+   (requires cudamat).
+-  Smooth optimal transport solvers (dual and semi-dual) for KL and
+   squared L2 regularizations [17].
+-  Non regularized Wasserstein barycenters [16] with LP solver (only
+   small scale).
 -  Bregman projections for Wasserstein barycenter [3] and unmixing [4].
 -  Optimal transport for domain adaptation with group lasso
    regularization [5]
@@ -29,6 +33,21 @@ It provides the following solvers:
 Some demonstrations (both in Python and Jupyter Notebook format) are
 available in the examples folder.
 
+Using and citing the toolbox
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+If you use this toolbox in your research and find it useful, please cite
+POT using the following bibtex reference:
+
+::
+
+    @misc{flamary2017pot,
+    title={POT Python Optimal Transport library},
+    author={Flamary, R{'e}mi and Courty, Nicolas},
+    url={https://github.com/rflamary/POT},
+    year={2017}
+    }
+
 Installation
 ------------
 
@@ -37,7 +56,7 @@ C++ compiler for using the EMD solver and relies on the following Python
 modules:
 
 -  Numpy (>=1.11)
--  Scipy (>=0.17)
+-  Scipy (>=1.0)
 -  Cython (>=0.23)
 -  Matplotlib (>=1.5)
 
@@ -212,20 +231,6 @@ languages):
 -  `Marco Cuturi <http://marcocuturi.net/>`__ (Sinkhorn Knopp in
    Matlab/Cuda)
 
-Using and citing the toolbox
-----------------------------
-
-If you use this toolbox in your research and find it useful, please cite
-POT using the following bibtex reference:
-
-::
-
-    @article{flamary2017pot,
-      title={POT Python Optimal Transport library},
-      author={Flamary, R{\'e}mi and Courty, Nicolas},
-      year={2017}
-    }
-
 Contributions and code of conduct
 ---------------------------------
 
@@ -320,6 +325,15 @@ Journal of Optimization Theory and Applications Vol 43.
 [15] Peyré, G., & Cuturi, M. (2018). `Computational Optimal
 Transport <https://arxiv.org/pdf/1803.00567.pdf>`__ .
 
+[16] Agueh, M., & Carlier, G. (2011). `Barycenters in the Wasserstein
+space <https://hal.archives-ouvertes.fr/hal-00637399/document>`__. SIAM
+Journal on Mathematical Analysis, 43(2), 904-924.
+
+[17] Blondel, M., Seguy, V., & Rolet, A. (2018). `Smooth and Sparse
+Optimal Transport <https://arxiv.org/abs/1710.06276>`__. Proceedings of
+the Twenty-First International Conference on Artificial Intelligence and
+Statistics (AISTATS).
+
 .. |PyPI version| image:: https://badge.fury.io/py/POT.svg
    :target: https://badge.fury.io/py/POT
 .. |Anaconda Cloud| image:: https://anaconda.org/conda-forge/pot/badges/version.svg

examples/plot_OT_1D_smooth.py

@@ -0,0 +1,110 @@
+# -*- coding: utf-8 -*-
+"""
+===========================
+1D smooth optimal transport
+===========================
+
+This example illustrates the computation of EMD, Sinkhorn and smooth OT plans
+and their visualization.
+
+"""
+
+# Author: Remi Flamary <remi.flamary@unice.fr>
+#
+# License: MIT License
+
+import numpy as np
+import matplotlib.pylab as pl
+import ot
+import ot.plot
+from ot.datasets import make_1D_gauss as gauss
+
+##############################################################################
+# Generate data
+# -------------
+
+
+#%% parameters
+
+n = 100  # nb bins
+
+# bin positions
+x = np.arange(n, dtype=np.float64)
+
+# Gaussian distributions
+a = gauss(n, m=20, s=5)  # m= mean, s= std
+b = gauss(n, m=60, s=10)
+
+# loss matrix
+M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1)))
+M /= M.max()
+
+
+##############################################################################
+# Plot distributions and loss matrix
+# ----------------------------------
+
+#%% plot the distributions
+
+pl.figure(1, figsize=(6.4, 3))
+pl.plot(x, a, 'b', label='Source distribution')
+pl.plot(x, b, 'r', label='Target distribution')
+pl.legend()
+
+#%% plot distributions and loss matrix
+
+pl.figure(2, figsize=(5, 5))
+ot.plot.plot1D_mat(a, b, M, 'Cost matrix M')
+
+##############################################################################
+# Solve EMD
+# ---------
+
+
+#%% EMD
+
+G0 = ot.emd(a, b, M)
+
+pl.figure(3, figsize=(5, 5))
+ot.plot.plot1D_mat(a, b, G0, 'OT matrix G0')
+
+##############################################################################
+# Solve Sinkhorn
+# --------------
+
+
+#%% Sinkhorn
+
+lambd = 2e-3
+Gs = ot.sinkhorn(a, b, M, lambd, verbose=True)
+
+pl.figure(4, figsize=(5, 5))
+ot.plot.plot1D_mat(a, b, Gs, 'OT matrix Sinkhorn')
+
+pl.show()
+
+##############################################################################
+# Solve Smooth OT
+# --------------
+
+
+#%% Smooth OT with KL regularization
+
+lambd = 2e-3
+Gsm = ot.smooth.smooth_ot_dual(a, b, M, lambd, reg_type='kl')
+
+pl.figure(5, figsize=(5, 5))
+ot.plot.plot1D_mat(a, b, Gsm, 'OT matrix Smooth OT KL reg.')
+
+pl.show()
+
+
+#%% Smooth OT with KL regularization
+
+lambd = 1e-1
+Gsm = ot.smooth.smooth_ot_dual(a, b, M, lambd, reg_type='l2')
+
+pl.figure(6, figsize=(5, 5))
+ot.plot.plot1D_mat(a, b, Gsm, 'OT matrix Smooth OT l2 reg.')
+
+pl.show()

ot/__init__.py

@@ -18,6 +18,8 @@
 from . import datasets
 from . import da
 from . import gromov
+from . import smooth
+
 
 # OT functions
 from .lp import emd, emd2

ot/lp/__init__.py

@@ -18,6 +18,8 @@
 from ..utils import parmap
 from .cvx import barycenter
 
+__all__=['emd', 'emd2', 'barycenter', 'cvx']
+
 
 def emd(a, b, M, numItermax=100000, log=False):
     """Solves the Earth Movers distance problem and returns the OT matrix