From 8d288976d398749c5261daca22ee4d772bd5b489 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Thu, 31 May 2018 11:08:01 +0200 Subject: add smooth.py + first tests --- ot/smooth.py | 405 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 405 insertions(+) create mode 100644 ot/smooth.py (limited to 'ot/smooth.py') diff --git a/ot/smooth.py b/ot/smooth.py new file mode 100644 index 0000000..b8c0a9a --- /dev/null +++ b/ot/smooth.py @@ -0,0 +1,405 @@ +#Copyright (c) 2018, Mathieu Blondel +#All rights reserved. +# +#Redistribution and use in source and binary forms, with or without +#modification, are permitted provided that the following conditions are met: +# +#1. Redistributions of source code must retain the above copyright notice, this +#list of conditions and the following disclaimer. +# +#2. Redistributions in binary form must reproduce the above copyright notice, +#this list of conditions and the following disclaimer in the documentation and/or +#other materials provided with the distribution. +# +#THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +#ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +#WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +#IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +#INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT +#NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, +#OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +#LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +#OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF +#THE POSSIBILITY OF SUCH DAMAGE. + +# Author: Mathieu Blondel +# Remi Flamary + +""" +Implementation of +Smooth and Sparse Optimal Transport. +Mathieu Blondel, Vivien Seguy, Antoine Rolet. +In Proc. of AISTATS 2018. +https://arxiv.org/abs/1710.06276 +""" + +import numpy as np +from scipy.optimize import minimize + + +def projection_simplex(V, z=1, axis=None): + """ Projection of x onto the simplex, scaled by z + + P(x; z) = argmin_{y >= 0, sum(y) = z} ||y - x||^2 + z: float or array + If array, len(z) must be compatible with V + axis: None or int + - axis=None: project V by P(V.ravel(); z) + - axis=1: project each V[i] by P(V[i]; z[i]) + - axis=0: project each V[:, j] by P(V[:, j]; z[j]) + """ + if axis == 1: + n_features = V.shape[1] + U = np.sort(V, axis=1)[:, ::-1] + z = np.ones(len(V)) * z + cssv = np.cumsum(U, axis=1) - z[:, np.newaxis] + ind = np.arange(n_features) + 1 + cond = U - cssv / ind > 0 + rho = np.count_nonzero(cond, axis=1) + theta = cssv[np.arange(len(V)), rho - 1] / rho + return np.maximum(V - theta[:, np.newaxis], 0) + + elif axis == 0: + return projection_simplex(V.T, z, axis=1).T + + else: + V = V.ravel().reshape(1, -1) + return projection_simplex(V, z, axis=1).ravel() + + +class Regularization(object): + + def __init__(self, gamma=1.0): + """ + + Parameters + ---------- + gamma: float + Regularization parameter. + We recover unregularized OT when gamma -> 0. + + """ + self.gamma = gamma + + def delta_Omega(X): + """ + Compute delta_Omega(X[:, j]) for each X[:, j]. + delta_Omega(x) = sup_{y >= 0} y^T x - Omega(y). + Parameters + ---------- + X: array, shape = len(a) x len(b) + Input array. + Returns + ------- + v: array, len(b) + Values: v[j] = delta_Omega(X[:, j]) + G: array, len(a) x len(b) + Gradients: G[:, j] = nabla delta_Omega(X[:, j]) + """ + raise NotImplementedError + + def max_Omega(X, b): + """ + Compute max_Omega_j(X[:, j]) for each X[:, j]. + max_Omega_j(x) = sup_{y >= 0, sum(y) = 1} y^T x - Omega(b[j] y) / b[j]. + + Parameters + ---------- + X: array, shape = len(a) x len(b) + Input array. + Returns + ------- + v: array, len(b) + Values: v[j] = max_Omega_j(X[:, j]) + G: array, len(a) x len(b) + Gradients: G[:, j] = nabla max_Omega_j(X[:, j]) + """ + raise NotImplementedError + + def Omega(T): + """ + Compute regularization term. + + Parameters + ---------- + T: array, shape = len(a) x len(b) + Input array. + Returns + ------- + value: float + Regularization term. + """ + raise NotImplementedError + + +class NegEntropy(Regularization): + + def delta_Omega(self, X): + G = np.exp(X / self.gamma - 1) + val = self.gamma * np.sum(G, axis=0) + return val, G + + def max_Omega(self, X, b): + max_X = np.max(X, axis=0) / self.gamma + exp_X = np.exp(X / self.gamma - max_X) + val = self.gamma * (np.log(np.sum(exp_X, axis=0)) + max_X) + val -= self.gamma * np.log(b) + G = exp_X / np.sum(exp_X, axis=0) + return val, G + + def Omega(self, T): + return self.gamma * np.sum(T * np.log(T)) + + +class SquaredL2(Regularization): + + def delta_Omega(self, X): + max_X = np.maximum(X, 0) + val = np.sum(max_X ** 2, axis=0) / (2 * self.gamma) + G = max_X / self.gamma + return val, G + + def max_Omega(self, X, b): + G = projection_simplex(X / (b * self.gamma), axis=0) + val = np.sum(X * G, axis=0) + val -= 0.5 * self.gamma * b * np.sum(G * G, axis=0) + return val, G + + def Omega(self, T): + return 0.5 * self.gamma * np.sum(T ** 2) + + +def dual_obj_grad(alpha, beta, a, b, C, regul): + """ + Compute objective value and gradients of dual objective. + + Parameters + ---------- + alpha: array, shape = len(a) + beta: array, shape = len(b) + Current iterate of dual potentials. + a: array, shape = len(a) + b: array, shape = len(b) + Input histograms (should be non-negative and sum to 1). + C: array, shape = len(a) x len(b) + Ground cost matrix. + regul: Regularization object + Should implement a delta_Omega(X) method. + + Returns + ------- + obj: float + Objective value (higher is better). + grad_alpha: array, shape = len(a) + Gradient w.r.t. alpha. + grad_beta: array, shape = len(b) + Gradient w.r.t. beta. + """ + obj = np.dot(alpha, a) + np.dot(beta, b) + grad_alpha = a.copy() + grad_beta = b.copy() + + # X[:, j] = alpha + beta[j] - C[:, j] + X = alpha[:, np.newaxis] + beta - C + + # val.shape = len(b) + # G.shape = len(a) x len(b) + val, G = regul.delta_Omega(X) + + obj -= np.sum(val) + grad_alpha -= G.sum(axis=1) + grad_beta -= G.sum(axis=0) + + return obj, grad_alpha, grad_beta + + +def solve_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): + """ + Solve the "smoothed" dual objective. + + Parameters + ---------- + a: array, shape = len(a) + b: array, shape = len(b) + Input histograms (should be non-negative and sum to 1). + C: array, shape = len(a) x len(b) + Ground cost matrix. + regul: Regularization object + Should implement a delta_Omega(X) method. + method: str + Solver to be used (passed to `scipy.optimize.minimize`). + tol: float + Tolerance parameter. + max_iter: int + Maximum number of iterations. + Returns + ------- + alpha: array, shape = len(a) + beta: array, shape = len(b) + Dual potentials. + """ + + def _func(params): + # Unpack alpha and beta. + alpha = params[:len(a)] + beta = params[len(a):] + + obj, grad_alpha, grad_beta = dual_obj_grad(alpha, beta, a, b, C, regul) + + # Pack grad_alpha and grad_beta. + grad = np.concatenate((grad_alpha, grad_beta)) + + # We need to maximize the dual. + return -obj, -grad + + # Unfortunately, `minimize` only supports functions whose argument is a + # vector. So, we need to concatenate alpha and beta. + alpha_init = np.zeros(len(a)) + beta_init = np.zeros(len(b)) + params_init = np.concatenate((alpha_init, beta_init)) + + res = minimize(_func, params_init, method=method, jac=True, + tol=tol, options=dict(maxiter=max_iter, disp=False)) + + alpha = res.x[:len(a)] + beta = res.x[len(a):] + + return alpha, beta, res + + +def semi_dual_obj_grad(alpha, a, b, C, regul): + """ + Compute objective value and gradient of semi-dual objective. + Parameters + ---------- + alpha: array, shape = len(a) + Current iterate of semi-dual potentials. + a: array, shape = len(a) + b: array, shape = len(b) + Input histograms (should be non-negative and sum to 1). + C: array, shape = len(a) x len(b) + Ground cost matrix. + regul: Regularization object + Should implement a max_Omega(X) method. + Returns + ------- + obj: float + Objective value (higher is better). + grad: array, shape = len(a) + Gradient w.r.t. alpha. + """ + obj = np.dot(alpha, a) + grad = a.copy() + + # X[:, j] = alpha - C[:, j] + X = alpha[:, np.newaxis] - C + + # val.shape = len(b) + # G.shape = len(a) x len(b) + val, G = regul.max_Omega(X, b) + + obj -= np.dot(b, val) + grad -= np.dot(G, b) + + return obj, grad + + +def solve_semi_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): + """ + Solve the "smoothed" semi-dual objective. + Parameters + ---------- + a: array, shape = len(a) + b: array, shape = len(b) + Input histograms (should be non-negative and sum to 1). + C: array, shape = len(a) x len(b) + Ground cost matrix. + regul: Regularization object + Should implement a max_Omega(X) method. + method: str + Solver to be used (passed to `scipy.optimize.minimize`). + tol: float + Tolerance parameter. + max_iter: int + Maximum number of iterations. + Returns + ------- + alpha: array, shape = len(a) + Semi-dual potentials. + """ + + def _func(alpha): + obj, grad = semi_dual_obj_grad(alpha, a, b, C, regul) + # We need to maximize the semi-dual. + return -obj, -grad + + alpha_init = np.zeros(len(a)) + + res = minimize(_func, alpha_init, method=method, jac=True, + tol=tol, options=dict(maxiter=max_iter, disp=False)) + + return res.x, res + + +def get_plan_from_dual(alpha, beta, C, regul): + """ + Retrieve optimal transportation plan from optimal dual potentials. + Parameters + ---------- + alpha: array, shape = len(a) + beta: array, shape = len(b) + Optimal dual potentials. + C: array, shape = len(a) x len(b) + Ground cost matrix. + regul: Regularization object + Should implement a delta_Omega(X) method. + Returns + ------- + T: array, shape = len(a) x len(b) + Optimal transportation plan. + """ + X = alpha[:, np.newaxis] + beta - C + return regul.delta_Omega(X)[1] + + +def get_plan_from_semi_dual(alpha, b, C, regul): + """ + Retrieve optimal transportation plan from optimal semi-dual potentials. + Parameters + ---------- + alpha: array, shape = len(a) + Optimal semi-dual potentials. + b: array, shape = len(b) + Second input histogram (should be non-negative and sum to 1). + C: array, shape = len(a) x len(b) + Ground cost matrix. + regul: Regularization object + Should implement a delta_Omega(X) method. + Returns + ------- + T: array, shape = len(a) x len(b) + Optimal transportation plan. + """ + X = alpha[:, np.newaxis] - C + return regul.max_Omega(X, b)[1] * b + + +def smooth_ot_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr=1e-9, + numItermax=500, log=False): + + if reg_type.lower()=='l2': + regul = SquaredL2(gamma=reg) + elif reg_type.lower() in ['entropic','negentropy','kl']: + regul = NegEntropy(gamma=reg) + + alpha, beta, res = solve_dual(a, b, M, regul, max_iter=numItermax,tol=stopThr) + G = get_plan_from_dual(alpha, beta, M, regul) + + if log: + log={'alpha':alpha,beta:'beta','res':res} + return G, log + else: + return G + + + -- cgit v1.2.3 From af99d3ff66062811a81454ea03e6b831a1292ae4 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Thu, 31 May 2018 11:26:32 +0200 Subject: add smooth.py + first tests --- ot/smooth.py | 13 +++++++++++++ test/test_smooth.py | 4 ++++ 2 files changed, 17 insertions(+) (limited to 'ot/smooth.py') diff --git a/ot/smooth.py b/ot/smooth.py index b8c0a9a..f8bb20a 100644 --- a/ot/smooth.py +++ b/ot/smooth.py @@ -85,10 +85,12 @@ class Regularization(object): """ Compute delta_Omega(X[:, j]) for each X[:, j]. delta_Omega(x) = sup_{y >= 0} y^T x - Omega(y). + Parameters ---------- X: array, shape = len(a) x len(b) Input array. + Returns ------- v: array, len(b) @@ -107,6 +109,7 @@ class Regularization(object): ---------- X: array, shape = len(a) x len(b) Input array. + Returns ------- v: array, len(b) @@ -124,6 +127,7 @@ class Regularization(object): ---------- T: array, shape = len(a) x len(b) Input array. + Returns ------- value: float @@ -232,6 +236,7 @@ def solve_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): Tolerance parameter. max_iter: int Maximum number of iterations. + Returns ------- alpha: array, shape = len(a) @@ -270,6 +275,7 @@ def solve_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): def semi_dual_obj_grad(alpha, a, b, C, regul): """ Compute objective value and gradient of semi-dual objective. + Parameters ---------- alpha: array, shape = len(a) @@ -281,6 +287,7 @@ def semi_dual_obj_grad(alpha, a, b, C, regul): Ground cost matrix. regul: Regularization object Should implement a max_Omega(X) method. + Returns ------- obj: float @@ -307,6 +314,7 @@ def semi_dual_obj_grad(alpha, a, b, C, regul): def solve_semi_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): """ Solve the "smoothed" semi-dual objective. + Parameters ---------- a: array, shape = len(a) @@ -322,6 +330,7 @@ def solve_semi_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): Tolerance parameter. max_iter: int Maximum number of iterations. + Returns ------- alpha: array, shape = len(a) @@ -344,6 +353,7 @@ def solve_semi_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): def get_plan_from_dual(alpha, beta, C, regul): """ Retrieve optimal transportation plan from optimal dual potentials. + Parameters ---------- alpha: array, shape = len(a) @@ -353,6 +363,7 @@ def get_plan_from_dual(alpha, beta, C, regul): Ground cost matrix. regul: Regularization object Should implement a delta_Omega(X) method. + Returns ------- T: array, shape = len(a) x len(b) @@ -365,6 +376,7 @@ def get_plan_from_dual(alpha, beta, C, regul): def get_plan_from_semi_dual(alpha, b, C, regul): """ Retrieve optimal transportation plan from optimal semi-dual potentials. + Parameters ---------- alpha: array, shape = len(a) @@ -375,6 +387,7 @@ def get_plan_from_semi_dual(alpha, b, C, regul): Ground cost matrix. regul: Regularization object Should implement a delta_Omega(X) method. + Returns ------- T: array, shape = len(a) x len(b) diff --git a/test/test_smooth.py b/test/test_smooth.py index f951bf9..4ca44f8 100644 --- a/test/test_smooth.py +++ b/test/test_smooth.py @@ -30,4 +30,8 @@ def test_smooth_ot_dual(): u, G.sum(1), atol=1e-05) # cf convergence sinkhorn np.testing.assert_allclose( u, G.sum(0), atol=1e-05) # cf convergence sinkhorn + + + G2 = ot.sinkhorn(u, u, M, 1, stopThr=1e-10) + np.testing.assert_allclose( G, G2 , atol=1e-05) \ No newline at end of file -- cgit v1.2.3 From e585d64dd08e5367350e70f23e81f9fd2d676a6b Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Thu, 31 May 2018 11:33:08 +0200 Subject: pep8 --- ot/smooth.py | 67 +++++++++++++++++++++++++---------------------------- test/test_smooth.py | 26 ++++++++++++--------- 2 files changed, 47 insertions(+), 46 deletions(-) (limited to 'ot/smooth.py') diff --git a/ot/smooth.py b/ot/smooth.py index f8bb20a..f4f4306 100644 --- a/ot/smooth.py +++ b/ot/smooth.py @@ -1,4 +1,4 @@ -#Copyright (c) 2018, Mathieu Blondel +#Copyright (c) 2018, Mathieu Blondel #All rights reserved. # #Redistribution and use in source and binary forms, with or without @@ -22,7 +22,7 @@ #OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF #THE POSSIBILITY OF SUCH DAMAGE. -# Author: Mathieu Blondel +# Author: Mathieu Blondel # Remi Flamary """ @@ -39,7 +39,7 @@ from scipy.optimize import minimize def projection_simplex(V, z=1, axis=None): """ Projection of x onto the simplex, scaled by z - + P(x; z) = argmin_{y >= 0, sum(y) = z} ||y - x||^2 z: float or array If array, len(z) must be compatible with V @@ -65,19 +65,19 @@ def projection_simplex(V, z=1, axis=None): else: V = V.ravel().reshape(1, -1) return projection_simplex(V, z, axis=1).ravel() - + class Regularization(object): def __init__(self, gamma=1.0): """ - + Parameters ---------- gamma: float Regularization parameter. We recover unregularized OT when gamma -> 0. - + """ self.gamma = gamma @@ -85,12 +85,12 @@ class Regularization(object): """ Compute delta_Omega(X[:, j]) for each X[:, j]. delta_Omega(x) = sup_{y >= 0} y^T x - Omega(y). - + Parameters ---------- X: array, shape = len(a) x len(b) Input array. - + Returns ------- v: array, len(b) @@ -104,12 +104,12 @@ class Regularization(object): """ Compute max_Omega_j(X[:, j]) for each X[:, j]. max_Omega_j(x) = sup_{y >= 0, sum(y) = 1} y^T x - Omega(b[j] y) / b[j]. - + Parameters ---------- X: array, shape = len(a) x len(b) Input array. - + Returns ------- v: array, len(b) @@ -122,12 +122,12 @@ class Regularization(object): def Omega(T): """ Compute regularization term. - + Parameters ---------- T: array, shape = len(a) x len(b) Input array. - + Returns ------- value: float @@ -176,7 +176,7 @@ class SquaredL2(Regularization): def dual_obj_grad(alpha, beta, a, b, C, regul): """ Compute objective value and gradients of dual objective. - + Parameters ---------- alpha: array, shape = len(a) @@ -189,7 +189,7 @@ def dual_obj_grad(alpha, beta, a, b, C, regul): Ground cost matrix. regul: Regularization object Should implement a delta_Omega(X) method. - + Returns ------- obj: float @@ -220,7 +220,7 @@ def dual_obj_grad(alpha, beta, a, b, C, regul): def solve_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): """ Solve the "smoothed" dual objective. - + Parameters ---------- a: array, shape = len(a) @@ -236,7 +236,7 @@ def solve_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): Tolerance parameter. max_iter: int Maximum number of iterations. - + Returns ------- alpha: array, shape = len(a) @@ -275,7 +275,7 @@ def solve_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): def semi_dual_obj_grad(alpha, a, b, C, regul): """ Compute objective value and gradient of semi-dual objective. - + Parameters ---------- alpha: array, shape = len(a) @@ -287,7 +287,7 @@ def semi_dual_obj_grad(alpha, a, b, C, regul): Ground cost matrix. regul: Regularization object Should implement a max_Omega(X) method. - + Returns ------- obj: float @@ -314,7 +314,7 @@ def semi_dual_obj_grad(alpha, a, b, C, regul): def solve_semi_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): """ Solve the "smoothed" semi-dual objective. - + Parameters ---------- a: array, shape = len(a) @@ -330,7 +330,7 @@ def solve_semi_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): Tolerance parameter. max_iter: int Maximum number of iterations. - + Returns ------- alpha: array, shape = len(a) @@ -353,7 +353,7 @@ def solve_semi_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): def get_plan_from_dual(alpha, beta, C, regul): """ Retrieve optimal transportation plan from optimal dual potentials. - + Parameters ---------- alpha: array, shape = len(a) @@ -363,7 +363,7 @@ def get_plan_from_dual(alpha, beta, C, regul): Ground cost matrix. regul: Regularization object Should implement a delta_Omega(X) method. - + Returns ------- T: array, shape = len(a) x len(b) @@ -376,7 +376,7 @@ def get_plan_from_dual(alpha, beta, C, regul): def get_plan_from_semi_dual(alpha, b, C, regul): """ Retrieve optimal transportation plan from optimal semi-dual potentials. - + Parameters ---------- alpha: array, shape = len(a) @@ -387,7 +387,7 @@ def get_plan_from_semi_dual(alpha, b, C, regul): Ground cost matrix. regul: Regularization object Should implement a delta_Omega(X) method. - + Returns ------- T: array, shape = len(a) x len(b) @@ -399,20 +399,17 @@ def get_plan_from_semi_dual(alpha, b, C, regul): def smooth_ot_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr=1e-9, numItermax=500, log=False): - - if reg_type.lower()=='l2': + + if reg_type.lower() == 'l2': regul = SquaredL2(gamma=reg) - elif reg_type.lower() in ['entropic','negentropy','kl']: - regul = NegEntropy(gamma=reg) - - alpha, beta, res = solve_dual(a, b, M, regul, max_iter=numItermax,tol=stopThr) + elif reg_type.lower() in ['entropic', 'negentropy', 'kl']: + regul = NegEntropy(gamma=reg) + + alpha, beta, res = solve_dual(a, b, M, regul, max_iter=numItermax, tol=stopThr) G = get_plan_from_dual(alpha, beta, M, regul) - + if log: - log={'alpha':alpha,beta:'beta','res':res} + log = {'alpha': alpha, beta: 'beta', 'res': res} return G, log else: return G - - - diff --git a/test/test_smooth.py b/test/test_smooth.py index 4ca44f8..e95b3fe 100644 --- a/test/test_smooth.py +++ b/test/test_smooth.py @@ -4,17 +4,14 @@ # # License: MIT License -import warnings - import numpy as np import ot -from ot.datasets import get_1D_gauss as gauss -import pytest def test_smooth_ot_dual(): - # test sinkhorn + + # get data n = 100 rng = np.random.RandomState(0) @@ -23,15 +20,22 @@ def test_smooth_ot_dual(): M = ot.dist(x, x) - G = ot.smooth.smooth_ot_dual(u, u, M, 1, stopThr=1e-10) + Gl2 = ot.smooth.smooth_ot_dual(u, u, M, 1, reg_type='l2', stopThr=1e-10) + + # check constratints + np.testing.assert_allclose( + u, Gl2.sum(1), atol=1e-05) # cf convergence sinkhorn + np.testing.assert_allclose( + u, Gl2.sum(0), atol=1e-05) # cf convergence sinkhorn + + # kl regyularisation + G = ot.smooth.smooth_ot_dual(u, u, M, 1, reg_type='kl', stopThr=1e-10) # check constratints np.testing.assert_allclose( u, G.sum(1), atol=1e-05) # cf convergence sinkhorn np.testing.assert_allclose( - u, G.sum(0), atol=1e-05) # cf convergence sinkhorn - - + u, G.sum(0), atol=1e-05) # cf convergence sinkhorn + G2 = ot.sinkhorn(u, u, M, 1, stopThr=1e-10) - np.testing.assert_allclose( G, G2 , atol=1e-05) - \ No newline at end of file + np.testing.assert_allclose(G, G2, atol=1e-05) -- cgit v1.2.3 From b188381a8e8190958e4e9aca64b7501cf0321406 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Thu, 31 May 2018 12:02:13 +0200 Subject: add semidual --- ot/smooth.py | 37 ++++++++++++++++++++++++++++++++++--- test/test_smooth.py | 34 +++++++++++++++++++++++++++++++++- 2 files changed, 67 insertions(+), 4 deletions(-) (limited to 'ot/smooth.py') diff --git a/ot/smooth.py b/ot/smooth.py index f4f4306..1a9972b 100644 --- a/ot/smooth.py +++ b/ot/smooth.py @@ -137,6 +137,7 @@ class Regularization(object): class NegEntropy(Regularization): + """ NegEntropy regularization """ def delta_Omega(self, X): G = np.exp(X / self.gamma - 1) @@ -156,6 +157,7 @@ class NegEntropy(Regularization): class SquaredL2(Regularization): + """ Squared L2 regularization """ def delta_Omega(self, X): max_X = np.maximum(X, 0) @@ -311,7 +313,8 @@ def semi_dual_obj_grad(alpha, a, b, C, regul): return obj, grad -def solve_semi_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): +def solve_semi_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500, + ): """ Solve the "smoothed" semi-dual objective. @@ -400,16 +403,44 @@ def get_plan_from_semi_dual(alpha, b, C, regul): def smooth_ot_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr=1e-9, numItermax=500, log=False): - if reg_type.lower() == 'l2': + if reg_type.lower() in ['l2', 'squaredl2']: regul = SquaredL2(gamma=reg) elif reg_type.lower() in ['entropic', 'negentropy', 'kl']: regul = NegEntropy(gamma=reg) + else: + raise NotImplementedError('Unknown regularization') + # solve dual alpha, beta, res = solve_dual(a, b, M, regul, max_iter=numItermax, tol=stopThr) + + # reconstruct transport matrix G = get_plan_from_dual(alpha, beta, M, regul) if log: - log = {'alpha': alpha, beta: 'beta', 'res': res} + log = {'alpha': alpha, 'beta': beta, 'res': res} + return G, log + else: + return G + + +def smooth_ot_semi_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr=1e-9, + numItermax=500, log=False): + + if reg_type.lower() in ['l2', 'squaredl2']: + regul = SquaredL2(gamma=reg) + elif reg_type.lower() in ['entropic', 'negentropy', 'kl']: + regul = NegEntropy(gamma=reg) + else: + raise NotImplementedError('Unknown regularization') + + # solve dual + alpha, res = solve_semi_dual(a, b, M, regul, max_iter=numItermax, tol=stopThr) + + # reconstruct transport matrix + G = get_plan_from_semi_dual(alpha, b, M, regul) + + if log: + log = {'alpha': alpha, 'res': res} return G, log else: return G diff --git a/test/test_smooth.py b/test/test_smooth.py index e95b3fe..37cc66e 100644 --- a/test/test_smooth.py +++ b/test/test_smooth.py @@ -20,7 +20,7 @@ def test_smooth_ot_dual(): M = ot.dist(x, x) - Gl2 = ot.smooth.smooth_ot_dual(u, u, M, 1, reg_type='l2', stopThr=1e-10) + Gl2, log = ot.smooth.smooth_ot_dual(u, u, M, 1, reg_type='l2', log=True, stopThr=1e-10) # check constratints np.testing.assert_allclose( @@ -39,3 +39,35 @@ def test_smooth_ot_dual(): G2 = ot.sinkhorn(u, u, M, 1, stopThr=1e-10) np.testing.assert_allclose(G, G2, atol=1e-05) + + +def test_smooth_ot_semi_dual(): + + # get data + n = 100 + rng = np.random.RandomState(0) + + x = rng.randn(n, 2) + u = ot.utils.unif(n) + + M = ot.dist(x, x) + + Gl2, log = ot.smooth.smooth_ot_semi_dual(u, u, M, 1, reg_type='l2', log=True, stopThr=1e-10) + + # check constratints + np.testing.assert_allclose( + u, Gl2.sum(1), atol=1e-05) # cf convergence sinkhorn + np.testing.assert_allclose( + u, Gl2.sum(0), atol=1e-05) # cf convergence sinkhorn + + # kl regyularisation + G = ot.smooth.smooth_ot_semi_dual(u, u, M, 1, reg_type='kl', stopThr=1e-10) + + # check constratints + np.testing.assert_allclose( + u, G.sum(1), atol=1e-05) # cf convergence sinkhorn + np.testing.assert_allclose( + u, G.sum(0), atol=1e-05) # cf convergence sinkhorn + + G2 = ot.sinkhorn(u, u, M, 1, stopThr=1e-10) + np.testing.assert_allclose(G, G2, atol=1e-05) -- cgit v1.2.3 From ed0d4171c6291a15360bdb8a955b0783585da749 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Thu, 31 May 2018 13:03:40 +0200 Subject: update readme --- README.md | 3 +++ examples/plot_OT_1D_smooth.py | 6 +++--- ot/smooth.py | 9 +++++++++ 3 files changed, 15 insertions(+), 3 deletions(-) (limited to 'ot/smooth.py') diff --git a/README.md b/README.md index 466c09c..b3068ee 100644 --- a/README.md +++ b/README.md @@ -15,6 +15,7 @@ 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). +* Smooth optimal transport solvers (dual and semi-dual) for KL and squared L2 regularization [17]. * Non regularized Wasserstein barycenters [16] with LP solver. * Bregman projections for Wasserstein barycenter [3] and unmixing [4]. * Optimal transport for domain adaptation with group lasso regularization [5] @@ -213,3 +214,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). diff --git a/examples/plot_OT_1D_smooth.py b/examples/plot_OT_1D_smooth.py index 4927617..dec84db 100644 --- a/examples/plot_OT_1D_smooth.py +++ b/examples/plot_OT_1D_smooth.py @@ -1,8 +1,8 @@ # -*- coding: utf-8 -*- """ -==================== -1D optimal transport -==================== +=========================== +1D smooth optimal transport +=========================== This example illustrates the computation of EMD, Sinkhorn and smooth OT plans and their visualization. diff --git a/ot/smooth.py b/ot/smooth.py index 1a9972b..b3649e9 100644 --- a/ot/smooth.py +++ b/ot/smooth.py @@ -22,6 +22,7 @@ #OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF #THE POSSIBILITY OF SUCH DAMAGE. + # Author: Mathieu Blondel # Remi Flamary @@ -31,6 +32,13 @@ Smooth and Sparse Optimal Transport. Mathieu Blondel, Vivien Seguy, Antoine Rolet. In Proc. of AISTATS 2018. https://arxiv.org/abs/1710.06276 + +[17] Blondel, M., Seguy, V., & Rolet, A. (2018). Smooth and Sparse Optimal +Transport. Proceedings of the Twenty-First International Conference on +Artificial Intelligence and Statistics (AISTATS). + +Original code from https://github.com/mblondel/smooth-ot/ + """ import numpy as np @@ -402,6 +410,7 @@ def get_plan_from_semi_dual(alpha, b, C, regul): def smooth_ot_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr=1e-9, numItermax=500, log=False): + if reg_type.lower() in ['l2', 'squaredl2']: regul = SquaredL2(gamma=reg) -- cgit v1.2.3 From 724984d3e13aa9390ef17a389fbd6ccf1f277d69 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Thu, 31 May 2018 13:05:37 +0200 Subject: pep8 --- examples/plot_OT_1D_smooth.py | 2 +- ot/smooth.py | 6 ++---- 2 files changed, 3 insertions(+), 5 deletions(-) (limited to 'ot/smooth.py') diff --git a/examples/plot_OT_1D_smooth.py b/examples/plot_OT_1D_smooth.py index dec84db..ec43279 100644 --- a/examples/plot_OT_1D_smooth.py +++ b/examples/plot_OT_1D_smooth.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- """ =========================== -1D smooth optimal transport +1D smooth optimal transport =========================== This example illustrates the computation of EMD, Sinkhorn and smooth OT plans diff --git a/ot/smooth.py b/ot/smooth.py index b3649e9..237b67a 100644 --- a/ot/smooth.py +++ b/ot/smooth.py @@ -22,7 +22,6 @@ #OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF #THE POSSIBILITY OF SUCH DAMAGE. - # Author: Mathieu Blondel # Remi Flamary @@ -33,8 +32,8 @@ Mathieu Blondel, Vivien Seguy, Antoine Rolet. In Proc. of AISTATS 2018. https://arxiv.org/abs/1710.06276 -[17] Blondel, M., Seguy, V., & Rolet, A. (2018). Smooth and Sparse Optimal -Transport. Proceedings of the Twenty-First International Conference on +[17] Blondel, M., Seguy, V., & Rolet, A. (2018). Smooth and Sparse Optimal +Transport. Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics (AISTATS). Original code from https://github.com/mblondel/smooth-ot/ @@ -410,7 +409,6 @@ def get_plan_from_semi_dual(alpha, b, C, regul): def smooth_ot_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr=1e-9, numItermax=500, log=False): - if reg_type.lower() in ['l2', 'squaredl2']: regul = SquaredL2(gamma=reg) -- cgit v1.2.3 From fb883fcddff31cc4e21c921ebeb9a550eaa48ff0 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Thu, 31 May 2018 13:36:12 +0200 Subject: proper documentation --- ot/smooth.py | 155 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 147 insertions(+), 8 deletions(-) (limited to 'ot/smooth.py') diff --git a/ot/smooth.py b/ot/smooth.py index 237b67a..4c46e73 100644 --- a/ot/smooth.py +++ b/ot/smooth.py @@ -226,7 +226,8 @@ def dual_obj_grad(alpha, beta, a, b, C, regul): return obj, grad_alpha, grad_beta -def solve_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): +def solve_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500, + verbose=False): """ Solve the "smoothed" dual objective. @@ -273,7 +274,7 @@ def solve_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500): params_init = np.concatenate((alpha_init, beta_init)) res = minimize(_func, params_init, method=method, jac=True, - tol=tol, options=dict(maxiter=max_iter, disp=False)) + tol=tol, options=dict(maxiter=max_iter, disp=verbose)) alpha = res.x[:len(a)] beta = res.x[len(a):] @@ -321,7 +322,7 @@ def semi_dual_obj_grad(alpha, a, b, C, regul): def solve_semi_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500, - ): + verbose=False): """ Solve the "smoothed" semi-dual objective. @@ -355,7 +356,7 @@ def solve_semi_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500, alpha_init = np.zeros(len(a)) res = minimize(_func, alpha_init, method=method, jac=True, - tol=tol, options=dict(maxiter=max_iter, disp=False)) + tol=tol, options=dict(maxiter=max_iter, disp=verbose)) return res.x, res @@ -408,7 +409,75 @@ def get_plan_from_semi_dual(alpha, b, C, regul): def smooth_ot_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr=1e-9, - numItermax=500, log=False): + numItermax=500, verbose=False, log=False): + r""" + Solve the regularized OT problem in the dual and return the OT matrix + + The function solves the smooth relaxed dual formulation (7) in [17]_ : + + .. math:: + \max_{\alpha,\beta}\quad a^T\alpha+b^T\beta-\sum_j\delta_\Omega(\alpha+\beta_j-\mathbf{m}_j) + + where : + + - :math:`\mathbf{m}_j` is the jth column of the cost matrix + - :math:`\delta_\Omega` is the convex conjugate of the regularization term :math:`\Omega` + - a and b are source and target weights (sum to 1) + + The OT matrix can is reconstructed from the gradient of :math:`\delta_\Omega` + (See [17]_ Proposition 1). + The optimization algorithm is using gradient decent (L-BFGS by default). + + + Parameters + ---------- + a : np.ndarray (ns,) + samples weights in the source domain + b : np.ndarray (nt,) or np.ndarray (nt,nbb) + samples in the target domain, compute sinkhorn with multiple targets + and fixed M if b is a matrix (return OT loss + dual variables in log) + M : np.ndarray (ns,nt) + loss matrix + reg : float + Regularization term >0 + reg_type : str + Regularization type, can be the following (default ='l2'): + - 'kl' : Kullback Leibler (~ Neg-entropy used in sinkhorn [2]_) + - 'l2' : Squared Euclidean regularization + method : str + Solver to use for scipy.optimize.minimize + 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 + ------- + gamma : (ns x nt) ndarray + Optimal transportation matrix for the given parameters + log : dict + log dictionary return only if log==True in parameters + + + References + ---------- + + .. [2] M. Cuturi, Sinkhorn Distances : Lightspeed Computation of Optimal Transport, Advances in Neural Information Processing Systems (NIPS) 26, 2013 + + .. [17] Blondel, M., Seguy, V., & Rolet, A. (2018). Smooth and Sparse Optimal Transport. Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics (AISTATS). + + See Also + -------- + ot.lp.emd : Unregularized OT + ot.sinhorn : Entropic regularized OT + ot.optim.cg : General regularized OT + + """ if reg_type.lower() in ['l2', 'squaredl2']: regul = SquaredL2(gamma=reg) @@ -418,7 +487,8 @@ def smooth_ot_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr=1e-9, raise NotImplementedError('Unknown regularization') # solve dual - alpha, beta, res = solve_dual(a, b, M, regul, max_iter=numItermax, tol=stopThr) + alpha, beta, res = solve_dual(a, b, M, regul, max_iter=numItermax, + tol=stopThr, verbose=verbose) # reconstruct transport matrix G = get_plan_from_dual(alpha, beta, M, regul) @@ -431,8 +501,77 @@ def smooth_ot_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr=1e-9, def smooth_ot_semi_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr=1e-9, - numItermax=500, log=False): + numItermax=500, verbose=False, log=False): + r""" + Solve the regularized OT problem in the semi-dual and return the OT matrix + The function solves the smooth relaxed dual formulation (10) in [17]_ : + + .. math:: + \max_{\alpha}\quad a^T\alpha-OT_\Omega^*(\alpha,b) + + where : + + .. math:: + OT_\Omega^*(\alpha,b)=\sum_j b_j + + - :math:`\mathbf{m}_j` is the jth column of the cost matrix + - :math:`OT_\Omega^*(\alpha,b)` is defined in Eq. (9) in [17] + - a and b are source and target weights (sum to 1) + + The OT matrix can is reconstructed using [17]_ Proposition 2. + The optimization algorithm is using gradient decent (L-BFGS by default). + + + Parameters + ---------- + a : np.ndarray (ns,) + samples weights in the source domain + b : np.ndarray (nt,) or np.ndarray (nt,nbb) + samples in the target domain, compute sinkhorn with multiple targets + and fixed M if b is a matrix (return OT loss + dual variables in log) + M : np.ndarray (ns,nt) + loss matrix + reg : float + Regularization term >0 + reg_type : str + Regularization type, can be the following (default ='l2'): + - 'kl' : Kullback Leibler (~ Neg-entropy used in sinkhorn [2]_) + - 'l2' : Squared Euclidean regularization + method : str + Solver to use for scipy.optimize.minimize + 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 + ------- + gamma : (ns x nt) ndarray + Optimal transportation matrix for the given parameters + log : dict + log dictionary return only if log==True in parameters + + + References + ---------- + + .. [2] M. Cuturi, Sinkhorn Distances : Lightspeed Computation of Optimal Transport, Advances in Neural Information Processing Systems (NIPS) 26, 2013 + + .. [17] Blondel, M., Seguy, V., & Rolet, A. (2018). Smooth and Sparse Optimal Transport. Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics (AISTATS). + + See Also + -------- + ot.lp.emd : Unregularized OT + ot.sinhorn : Entropic regularized OT + ot.optim.cg : General regularized OT + + """ if reg_type.lower() in ['l2', 'squaredl2']: regul = SquaredL2(gamma=reg) elif reg_type.lower() in ['entropic', 'negentropy', 'kl']: @@ -444,7 +583,7 @@ def smooth_ot_semi_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr= alpha, res = solve_semi_dual(a, b, M, regul, max_iter=numItermax, tol=stopThr) # reconstruct transport matrix - G = get_plan_from_semi_dual(alpha, b, M, regul) + G = get_plan_from_semi_dual(alpha, b, M, regul, verbose=verbose) if log: log = {'alpha': alpha, 'res': res} -- cgit v1.2.3 From d370f18aaac4ed6903e3b4251c8b9c4685c53cc3 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Thu, 31 May 2018 13:40:45 +0200 Subject: bug verbose semi-dual --- ot/smooth.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) (limited to 'ot/smooth.py') diff --git a/ot/smooth.py b/ot/smooth.py index 4c46e73..dd734b3 100644 --- a/ot/smooth.py +++ b/ot/smooth.py @@ -580,10 +580,11 @@ def smooth_ot_semi_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr= raise NotImplementedError('Unknown regularization') # solve dual - alpha, res = solve_semi_dual(a, b, M, regul, max_iter=numItermax, tol=stopThr) + alpha, res = solve_semi_dual(a, b, M, regul, max_iter=numItermax, + tol=stopThr, verbose=verbose) # reconstruct transport matrix - G = get_plan_from_semi_dual(alpha, b, M, regul, verbose=verbose) + G = get_plan_from_semi_dual(alpha, b, M, regul) if log: log = {'alpha': alpha, 'res': res} -- cgit v1.2.3 From 8046b8c424d7b80f520e212e2bf8de41cb624aab Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Thu, 31 May 2018 13:45:11 +0200 Subject: pep8 --- README.md | 2 +- ot/smooth.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) (limited to 'ot/smooth.py') diff --git a/README.md b/README.md index b3068ee..6d3286b 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ 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). * Smooth optimal transport solvers (dual and semi-dual) for KL and squared L2 regularization [17]. -* Non regularized Wasserstein barycenters [16] with LP solver. +* 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]. diff --git a/ot/smooth.py b/ot/smooth.py index dd734b3..87d61ef 100644 --- a/ot/smooth.py +++ b/ot/smooth.py @@ -580,7 +580,7 @@ def smooth_ot_semi_dual(a, b, M, reg, reg_type='l2', method="L-BFGS-B", stopThr= raise NotImplementedError('Unknown regularization') # solve dual - alpha, res = solve_semi_dual(a, b, M, regul, max_iter=numItermax, + alpha, res = solve_semi_dual(a, b, M, regul, max_iter=numItermax, tol=stopThr, verbose=verbose) # reconstruct transport matrix -- cgit v1.2.3 From ecb093b95ddbeaeb800b443d3a5c9d5c24c5897c Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Mon, 11 Jun 2018 13:00:57 +0200 Subject: ad documentation class Regularization --- ot/smooth.py | 7 +++++++ 1 file changed, 7 insertions(+) (limited to 'ot/smooth.py') diff --git a/ot/smooth.py b/ot/smooth.py index 87d61ef..5a8e4b5 100644 --- a/ot/smooth.py +++ b/ot/smooth.py @@ -75,6 +75,13 @@ def projection_simplex(V, z=1, axis=None): class Regularization(object): + """Base class for Regularization objects + + Notes + ----- + This class is not intended for direct use but as aparent for true + regularizatiojn implementation. + """ def __init__(self, gamma=1.0): """ -- cgit v1.2.3