From e885d78cc9608d791a9d1561d2f4e0b783ba0761 Mon Sep 17 00:00:00 2001 From: Kilian Fatras Date: Tue, 28 Aug 2018 17:24:07 -0700 Subject: debug sgd dual --- README.md | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) (limited to 'README.md') diff --git a/README.md b/README.md index 8e8dcd4..677a23b 100644 --- a/README.md +++ b/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]. @@ -26,12 +27,24 @@ 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 The library has been tested on Linux, MacOSX and Windows. It requires a 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) @@ -158,16 +171,7 @@ This toolbox benefit a lot from open source research and we would like to thank * [Nicolas Bonneel](http://liris.cnrs.fr/~nbonneel/) ( C++ code for EMD) * [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 Every contribution is welcome and should respect the [contribution guidelines](CONTRIBUTING.md). Each member of the project is expected to follow the [code of conduct](CODE_OF_CONDUCT.md). @@ -216,7 +220,7 @@ You can also post bug reports and feature requests in Github issues. Make sure t [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/pdf/1710.06276.pdf). Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics (AISTATS). +[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). [18] Genevay, A., Cuturi, M., Peyré, G. & Bach, F. (2016) [Stochastic Optimization for Large-scale Optimal Transport](arXiv preprint arxiv:1605.08527). Advances in Neural Information Processing Systems (2016). -- cgit v1.2.3 From b13feb07eaff4d971b749663652e5f8811c1992c Mon Sep 17 00:00:00 2001 From: Kilian Fatras Date: Tue, 28 Aug 2018 18:18:37 -0700 Subject: added gaussian test --- README.md | 7 +++++-- test/test_stochastic.py | 47 ++++++++++++++++++++++++++++++++++++++--------- 2 files changed, 43 insertions(+), 11 deletions(-) (limited to 'README.md') diff --git a/README.md b/README.md index 677a23b..1d3b097 100644 --- a/README.md +++ b/README.md @@ -24,6 +24,7 @@ It provides the following solvers: * Wasserstein Discriminant Analysis [11] (requires autograd + pymanopt). * Gromov-Wasserstein distances and barycenters ([13] and regularized [12]) * Stochastic Optimization for Large-scale Optimal Transport (semi-dual problem [18] and dual problem [19]) +* Non regularized free support Wasserstein barycenters [20]. Some demonstrations (both in Python and Jupyter Notebook format) are available in the examples folder. @@ -163,7 +164,7 @@ The contributors to this library are: * [Stanislas Chambon](https://slasnista.github.io/) * [Antoine Rolet](https://arolet.github.io/) * Erwan Vautier (Gromov-Wasserstein) -* [Kilian Fatras](https://kilianfatras.github.io/) (Stochastic optimization) +* [Kilian Fatras](https://kilianfatras.github.io/) This toolbox benefit a lot from open source research and we would like to thank the following persons for providing some code (in various languages): @@ -222,6 +223,8 @@ You can also post bug reports and feature requests in Github issues. Make sure t [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). -[18] Genevay, A., Cuturi, M., Peyré, G. & Bach, F. (2016) [Stochastic Optimization for Large-scale Optimal Transport](arXiv preprint arxiv:1605.08527). Advances in Neural Information Processing Systems (2016). +[18] Genevay, A., Cuturi, M., Peyré, G. & Bach, F. (2016) [Stochastic Optimization for Large-scale Optimal Transport](https://arxiv.org/abs/1605.08527). Advances in Neural Information Processing Systems (2016). [19] Seguy, V., Bhushan Damodaran, B., Flamary, R., Courty, N., Rolet, A.& Blondel, M. [Large-scale Optimal Transport and Mapping Estimation](https://arxiv.org/pdf/1711.02283.pdf). International Conference on Learning Representation (2018) + +[20] Cuturi, M. and Doucet, A. (2014) [Fast Computation of Wasserstein Barycenters](http://proceedings.mlr.press/v32/cuturi14.html). International Conference in Machine Learning diff --git a/test/test_stochastic.py b/test/test_stochastic.py index f315c88..88ad666 100644 --- a/test/test_stochastic.py +++ b/test/test_stochastic.py @@ -137,8 +137,8 @@ def test_stochastic_dual_sgd(): # test sgd n = 10 reg = 1 - numItermax = 300000 - batch_size = 8 + numItermax = 15000 + batch_size = 10 rng = np.random.RandomState(0) x = rng.randn(n, 2) @@ -151,9 +151,9 @@ def test_stochastic_dual_sgd(): # check constratints np.testing.assert_allclose( - u, G.sum(1), atol=1e-02) # cf convergence sgd + u, G.sum(1), atol=1e-04) # cf convergence sgd np.testing.assert_allclose( - u, G.sum(0), atol=1e-02) # cf convergence sgd + u, G.sum(0), atol=1e-04) # cf convergence sgd ############################################################################# @@ -168,10 +168,11 @@ def test_dual_sgd_sinkhorn(): # test all dual algorithms n = 10 reg = 1 - nb_iter = 300000 - batch_size = 8 + nb_iter = 150000 + batch_size = 10 rng = np.random.RandomState(0) +# Test uniform x = rng.randn(n, 2) u = ot.utils.unif(n) zero = np.zeros(n) @@ -184,8 +185,36 @@ def test_dual_sgd_sinkhorn(): # check constratints np.testing.assert_allclose( - zero, (G_sgd - G_sinkhorn).sum(1), atol=1e-02) # cf convergence sgd + zero, (G_sgd - G_sinkhorn).sum(1), atol=1e-04) # cf convergence sgd + np.testing.assert_allclose( + zero, (G_sgd - G_sinkhorn).sum(0), atol=1e-04) # cf convergence sgd + np.testing.assert_allclose( + G_sgd, G_sinkhorn, atol=1e-04) # cf convergence sgd + +# Test gaussian + n = 30 + n_source = n + n_target = n + reg = 1 + numItermax = 150000 + batch_size = 30 + + a = ot.datasets.get_1D_gauss(n_source, m=15, s=5) # m= mean, s= std + b = ot.datasets.get_1D_gauss(n_target, m=15, s=5) + X_source = np.arange(n_source,dtype=np.float64) + Y_target = np.arange(n_target,dtype=np.float64) + M = ot.dist(X_source.reshape((n_source, 1)), Y_target.reshape((n_target, 1))) + M /= M.max() + + G_sgd = ot.stochastic.solve_dual_entropic(a, b, M, reg, batch_size, + numItermax=nb_iter) + + G_sinkhorn = ot.sinkhorn(a, b, M, reg) + + # check constratints + np.testing.assert_allclose( + zero, (G_sgd - G_sinkhorn).sum(1), atol=1e-04) # cf convergence sgd np.testing.assert_allclose( - zero, (G_sgd - G_sinkhorn).sum(0), atol=1e-02) # cf convergence sgd + zero, (G_sgd - G_sinkhorn).sum(0), atol=1e-04) # cf convergence sgd np.testing.assert_allclose( - G_sgd, G_sinkhorn, atol=1e-02) # cf convergence sgd + G_sgd, G_sinkhorn, atol=1e-04) # cf convergence sgd -- cgit v1.2.3 From d99abf078537acf6cf49480b9790a9c450889031 Mon Sep 17 00:00:00 2001 From: Nicolas Courty Date: Fri, 7 Sep 2018 11:58:42 +0200 Subject: Wasserstein convolutional barycenter --- README.md | 4 +- data/duck.png | Bin 0 -> 5112 bytes data/heart.png | Bin 0 -> 5225 bytes data/redcross.png | Bin 0 -> 1683 bytes data/tooth.png | Bin 0 -> 4931 bytes examples/plot_convolutional_barycenter.py | 92 ++++++++++++++++++++++++++ ot/bregman.py | 106 ++++++++++++++++++++++++++++++ 7 files changed, 201 insertions(+), 1 deletion(-) create mode 100644 data/duck.png create mode 100644 data/heart.png create mode 100644 data/redcross.png create mode 100644 data/tooth.png create mode 100644 examples/plot_convolutional_barycenter.py (limited to 'README.md') diff --git a/README.md b/README.md index dded582..1105362 100644 --- a/README.md +++ b/README.md @@ -227,4 +227,6 @@ You can also post bug reports and feature requests in Github issues. Make sure t [19] Seguy, V., Bhushan Damodaran, B., Flamary, R., Courty, N., Rolet, A.& Blondel, M. [Large-scale Optimal Transport and Mapping Estimation](https://arxiv.org/pdf/1711.02283.pdf). International Conference on Learning Representation (2018) -[20] Cuturi, M. and Doucet, A. (2014) [Fast Computation of Wasserstein Barycenters](http://proceedings.mlr.press/v32/cuturi14.html). International Conference in Machine Learning \ No newline at end of file +[20] Cuturi, M. and Doucet, A. (2014) [Fast Computation of Wasserstein Barycenters](http://proceedings.mlr.press/v32/cuturi14.html). International Conference in Machine Learning + +[21] Solomon, J., De Goes, F., Peyré, G., Cuturi, M., Butscher, A., Nguyen, A. & Guibas, L. (2015). [Convolutional wasserstein distances: Efficient optimal transportation on geometric domains](https://dl.acm.org/citation.cfm?id=2766963). ACM Transactions on Graphics (TOG), 34(4), 66. \ No newline at end of file diff --git a/data/duck.png b/data/duck.png new file mode 100644 index 0000000..9181697 Binary files /dev/null and b/data/duck.png differ diff --git a/data/heart.png b/data/heart.png new file mode 100644 index 0000000..44a6385 Binary files /dev/null and b/data/heart.png differ diff --git a/data/redcross.png b/data/redcross.png new file mode 100644 index 0000000..8d0a6fa Binary files /dev/null and b/data/redcross.png differ diff --git a/data/tooth.png b/data/tooth.png new file mode 100644 index 0000000..cd92c9d Binary files /dev/null and b/data/tooth.png differ diff --git a/examples/plot_convolutional_barycenter.py b/examples/plot_convolutional_barycenter.py new file mode 100644 index 0000000..d231da9 --- /dev/null +++ b/examples/plot_convolutional_barycenter.py @@ -0,0 +1,92 @@ + +#%% +# -*- coding: utf-8 -*- +""" +============================================ +Convolutional Wasserstein Barycenter example +============================================ + +This example is designed to illustrate how the Convolutional Wasserstein Barycenter +function of POT works. +""" + +# Author: Nicolas Courty +# +# License: MIT License + + +import numpy as np +import pylab as pl +import ot + +############################################################################## +# Data preparation +# ---------------- +# +# The four distributions are constructed from 4 simple images + + +f1 = 1 - pl.imread('../data/redcross.png')[:, :, 2] +f2 = 1 - pl.imread('../data/duck.png')[:, :, 2] +f3 = 1 - pl.imread('../data/heart.png')[:, :, 2] +f4 = 1 - pl.imread('../data/tooth.png')[:, :, 2] + +A = [] +f1=f1/np.sum(f1) +f2=f2/np.sum(f2) +f3=f3/np.sum(f3) +f4=f4/np.sum(f4) +A.append(f1) +A.append(f2) +A.append(f3) +A.append(f4) +A=np.array(A) + +nb_images = 5 + +# those are the four corners coordinates that will be interpolated by bilinear +# interpolation +v1=np.array((1,0,0,0)) +v2=np.array((0,1,0,0)) +v3=np.array((0,0,1,0)) +v4=np.array((0,0,0,1)) + + +############################################################################## +# Barycenter computation and visualization +# ---------------------------------------- +# + +pl.figure(figsize=(10,10)) +pl.title('Convolutional Wasserstein Barycenters in POT') +cm='Blues' +# regularization parameter +reg=0.004 +for i in range(nb_images): + for j in range(nb_images): + pl.subplot(nb_images,nb_images,i*nb_images+j+1) + tx=float(i)/(nb_images-1) + ty=float(j)/(nb_images-1) + + # weights are constructed by bilinear interpolation + tmp1=(1-tx)*v1+tx*v2 + tmp2=(1-tx)*v3+tx*v4 + weights=(1-ty)*tmp1+ty*tmp2 + + if i==0 and j==0: + pl.imshow(f1,cmap=cm) + pl.axis('off') + elif i==0 and j==(nb_images-1): + pl.imshow(f3,cmap=cm) + pl.axis('off') + elif i==(nb_images-1) and j==0: + pl.imshow(f2,cmap=cm) + pl.axis('off') + elif i==(nb_images-1) and j==(nb_images-1): + pl.imshow(f4,cmap=cm) + pl.axis('off') + else: + # call to barycenter computation + pl.imshow(ot.convolutional_barycenter2d(A,reg,weights),cmap=cm) + pl.axis('off') +pl.show() \ No newline at end of file diff --git a/ot/bregman.py b/ot/bregman.py index c755f51..05f4d9d 100644 --- a/ot/bregman.py +++ b/ot/bregman.py @@ -918,6 +918,112 @@ def barycenter(A, M, reg, weights=None, numItermax=1000, else: return geometricBar(weights, UKv) +def convolutional_barycenter2d(A,reg,weights=None,numItermax = 10000, stopThr=1e-9, verbose=False, log=False): + """Compute the entropic regularized wasserstein barycenter of distributions A + where A is a collection of 2D images. + + The function solves the following optimization problem: + + .. math:: + \mathbf{a} = arg\min_\mathbf{a} \sum_i W_{reg}(\mathbf{a},\mathbf{a}_i) + + where : + + - :math:`W_{reg}(\cdot,\cdot)` is the entropic regularized Wasserstein distance (see ot.bregman.sinkhorn) + - :math:`\mathbf{a}_i` are training distributions (2D images) in the mast two dimensions of matrix :math:`\mathbf{A}` + - reg is the regularization strength scalar value + + The algorithm used for solving the problem is the Sinkhorn-Knopp matrix scaling algorithm as proposed in [21]_ + + Parameters + ---------- + A : np.ndarray (n,w,h) + n distributions (2D images) of size w x h + reg : float + Regularization term >0 + weights : np.ndarray (n,) + Weights of each image on the simplex (barycentric coodinates) + numItermax : int, optional + Max number of iterations + stopThr : float, optional + Stop threshol on error (>0) + verbose : bool, optional + Print information along iterations + log : bool, optional + record log if True + + + Returns + ------- + a : (w,h) ndarray + 2D Wasserstein barycenter + log : dict + log dictionary return only if log==True in parameters + + + References + ---------- + + .. [21] Solomon, J., De Goes, F., Peyré, G., Cuturi, M., Butscher, A., Nguyen, A. & Guibas, L. (2015). + Convolutional wasserstein distances: Efficient optimal transportation on geometric domains + ACM Transactions on Graphics (TOG), 34(4), 66 + + + """ + + if weights is None: + weights = np.ones(A.shape[0]) / A.shape[0] + else: + assert(len(weights) == A.shape[0]) + + if log: + log = {'err': []} + + b=np.zeros_like(A[0,:,:]) + U=np.ones_like(A) + KV=np.ones_like(A) + threshold = 1e-30 # in order to avoids numerical precision issues + + cpt = 0 + err=1 + + # build the convolution operator + t = np.linspace(0,1,A.shape[1]) + [Y,X] = np.meshgrid(t,t) + xi1 = np.exp(-(X-Y)**2/reg) + K = lambda x: np.dot(np.dot(xi1,x),xi1) + + while (err>stopThr and cpt