From 94eb4a24fe09c8a193933d8c7f48a657da4e6c52 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Fri, 11 May 2018 16:07:03 +0200 Subject: update documentation in bregman --- ot/bregman.py | 2 ++ 1 file changed, 2 insertions(+) (limited to 'ot') diff --git a/ot/bregman.py b/ot/bregman.py index 07b8660..9c84aed 100644 --- a/ot/bregman.py +++ b/ot/bregman.py @@ -844,6 +844,8 @@ def barycenter(A, M, reg, weights=None, numItermax=1000, loss matrix for OT reg : float Regularization term >0 + weights : np.ndarray (n,) + Weights of each histogram i_i on the simplex numItermax : int, optional Max number of iterations stopThr : float, optional -- cgit v1.2.3 From 060d9046b291c76244deab2d78ee8356a294e91f Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Fri, 11 May 2018 16:56:47 +0200 Subject: add cvx barycenter solver --- ot/lp/cvx.py | 138 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 138 insertions(+) create mode 100644 ot/lp/cvx.py (limited to 'ot') diff --git a/ot/lp/cvx.py b/ot/lp/cvx.py new file mode 100644 index 0000000..4d08916 --- /dev/null +++ b/ot/lp/cvx.py @@ -0,0 +1,138 @@ +# -*- coding: utf-8 -*- +""" +LP solvers for optimal transport using cvxopt +""" + +# Author: Remi Flamary +# +# License: MIT License + +import numpy as np +import scipy as sp +import scipy.sparse as sps + +try: + import cvxopt + from cvxopt import solvers, matrix, sparse, spmatrix +except ImportError: + cvxopt=False + +def scipy_sparse_to_spmatrix(A): + """Efficient conversion from scipy sparse matrix to cvxopt sparse matrix""" + coo = A.tocoo() + SP = spmatrix(coo.data.tolist(), coo.row.tolist(), coo.col.tolist(), size=A.shape) + return SP + +def barycenter(A, M, weights=None, verbose=False, log=False,solver='interior-point'): + """Compute the entropic regularized wasserstein barycenter of distributions A + + The function solves the following optimization problem [16]: + + .. math:: + \mathbf{a} = arg\min_\mathbf{a} \sum_i W_{1}(\mathbf{a},\mathbf{a}_i) + + where : + + - :math:`W_1(\cdot,\cdot)` is the Wasserstein distance (see ot.emd.sinkhorn) + - :math:`\mathbf{a}_i` are training distributions in the columns of matrix :math:`\mathbf{A}` + + The linear program is solved using the default cvxopt solver if installed. + If cvxopt is not installed it uses the lp solver from scipy.optimize. + + Parameters + ---------- + A : np.ndarray (d,n) + n training distributions of size d + M : np.ndarray (d,d) + loss matrix for OT + reg : float + Regularization term >0 + weights : np.ndarray (n,) + Weights of each histogram i_i on the simplex + verbose : bool, optional + Print information along iterations + log : bool, optional + record log if True + solver : string, optional + the solver used, default 'interior-point' use the lp solver from + scipy.optimize. None, or 'glpk' or 'mosek' use the solver from cvxopt. + + Returns + ------- + a : (d,) ndarray + Wasserstein barycenter + log : dict + log dictionary return only if log==True in parameters + + + References + ---------- + + .. [16] Agueh, M., & Carlier, G. (2011). Barycenters in the Wasserstein space. SIAM Journal on Mathematical Analysis, 43(2), 904-924. + + + + """ + + if weights is None: + weights = np.ones(A.shape[1]) / A.shape[1] + else: + assert(len(weights) == A.shape[1]) + + n_distributions=A.shape[1] + n=A.shape[0] + + n2=n*n + c=np.zeros((0)) + b_eq1=np.zeros((0)) + for i in range(n_distributions): + c=np.concatenate((c,M.ravel()*weights[i])) + b_eq1=np.concatenate((b_eq1,A[:,i])) + c=np.concatenate((c,np.zeros(n))) + + lst_idiag1=[sps.kron(sps.eye(n),np.ones((1,n))) for i in range(n_distributions)] + # row constraints + A_eq1=sps.hstack((sps.block_diag(lst_idiag1),sps.coo_matrix((n_distributions*n,n)))) + + # columns constraints + lst_idiag2=[] + lst_eye=[] + for i in range(n_distributions): + if i==0: + lst_idiag2.append(sps.kron(np.ones((1,n)),sps.eye(n))) + lst_eye.append(-sps.eye(n)) + else: + lst_idiag2.append(sps.kron(np.ones((1,n)),sps.eye(n-1,n))) + lst_eye.append(-sps.eye(n-1,n)) + + A_eq2=sps.hstack((sps.block_diag(lst_idiag2),sps.vstack(lst_eye))) + b_eq2=np.zeros((A_eq2.shape[0])) + + # full problem + A_eq=sps.vstack((A_eq1,A_eq2)) + b_eq=np.concatenate((b_eq1,b_eq2)) + + if not cvxopt or solver in ['interior-point']: # cvxopt not installed or simplex/interior point + + if solver is None: + solver='interior-point' + + options={'sparse':True,'disp': verbose} + sol=sp.optimize.linprog(c,A_eq=A_eq,b_eq=b_eq,method=solver,options=options) + x=sol.x + b=x[-n:] + + else: + + h=np.zeros((n_distributions*n2+n)) + G=-sps.eye(n_distributions*n2+n) + + sol=solvers.lp(matrix(c),scipy_sparse_to_spmatrix(G),matrix(h),A=scipy_sparse_to_spmatrix(A_eq),b=matrix(b_eq),solver=solver) + + x=np.array(sol['x']) + b=x[-n:].ravel() + + if log: + return b, sol + else: + return b \ No newline at end of file -- cgit v1.2.3 From 3aee908ad42d65897f1916de6eab84921ac94a10 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Fri, 11 May 2018 16:58:06 +0200 Subject: pep8 --- ot/bregman.py | 2 +- ot/lp/cvx.py | 104 ++++++++++++++++++++++++++++++---------------------------- 2 files changed, 54 insertions(+), 52 deletions(-) (limited to 'ot') diff --git a/ot/bregman.py b/ot/bregman.py index 9c84aed..e788ef5 100644 --- a/ot/bregman.py +++ b/ot/bregman.py @@ -845,7 +845,7 @@ def barycenter(A, M, reg, weights=None, numItermax=1000, reg : float Regularization term >0 weights : np.ndarray (n,) - Weights of each histogram i_i on the simplex + Weights of each histogram i_i on the simplex numItermax : int, optional Max number of iterations stopThr : float, optional diff --git a/ot/lp/cvx.py b/ot/lp/cvx.py index 4d08916..93097d1 100644 --- a/ot/lp/cvx.py +++ b/ot/lp/cvx.py @@ -15,7 +15,8 @@ try: import cvxopt from cvxopt import solvers, matrix, sparse, spmatrix except ImportError: - cvxopt=False + cvxopt = False + def scipy_sparse_to_spmatrix(A): """Efficient conversion from scipy sparse matrix to cvxopt sparse matrix""" @@ -23,7 +24,8 @@ def scipy_sparse_to_spmatrix(A): SP = spmatrix(coo.data.tolist(), coo.row.tolist(), coo.col.tolist(), size=A.shape) return SP -def barycenter(A, M, weights=None, verbose=False, log=False,solver='interior-point'): + +def barycenter(A, M, weights=None, verbose=False, log=False, solver='interior-point'): """Compute the entropic regularized wasserstein barycenter of distributions A The function solves the following optimization problem [16]: @@ -36,7 +38,7 @@ def barycenter(A, M, weights=None, verbose=False, log=False,solver='interior-poi - :math:`W_1(\cdot,\cdot)` is the Wasserstein distance (see ot.emd.sinkhorn) - :math:`\mathbf{a}_i` are training distributions in the columns of matrix :math:`\mathbf{A}` - The linear program is solved using the default cvxopt solver if installed. + The linear program is solved using the default cvxopt solver if installed. If cvxopt is not installed it uses the lp solver from scipy.optimize. Parameters @@ -48,13 +50,13 @@ def barycenter(A, M, weights=None, verbose=False, log=False,solver='interior-poi reg : float Regularization term >0 weights : np.ndarray (n,) - Weights of each histogram i_i on the simplex + Weights of each histogram i_i on the simplex verbose : bool, optional Print information along iterations log : bool, optional record log if True solver : string, optional - the solver used, default 'interior-point' use the lp solver from + the solver used, default 'interior-point' use the lp solver from scipy.optimize. None, or 'glpk' or 'mosek' use the solver from cvxopt. Returns @@ -78,61 +80,61 @@ def barycenter(A, M, weights=None, verbose=False, log=False,solver='interior-poi weights = np.ones(A.shape[1]) / A.shape[1] else: assert(len(weights) == A.shape[1]) - - n_distributions=A.shape[1] - n=A.shape[0] - - n2=n*n - c=np.zeros((0)) - b_eq1=np.zeros((0)) + + n_distributions = A.shape[1] + n = A.shape[0] + + n2 = n * n + c = np.zeros((0)) + b_eq1 = np.zeros((0)) for i in range(n_distributions): - c=np.concatenate((c,M.ravel()*weights[i])) - b_eq1=np.concatenate((b_eq1,A[:,i])) - c=np.concatenate((c,np.zeros(n))) - - lst_idiag1=[sps.kron(sps.eye(n),np.ones((1,n))) for i in range(n_distributions)] + c = np.concatenate((c, M.ravel() * weights[i])) + b_eq1 = np.concatenate((b_eq1, A[:, i])) + c = np.concatenate((c, np.zeros(n))) + + lst_idiag1 = [sps.kron(sps.eye(n), np.ones((1, n))) for i in range(n_distributions)] # row constraints - A_eq1=sps.hstack((sps.block_diag(lst_idiag1),sps.coo_matrix((n_distributions*n,n)))) - + A_eq1 = sps.hstack((sps.block_diag(lst_idiag1), sps.coo_matrix((n_distributions * n, n)))) + # columns constraints - lst_idiag2=[] - lst_eye=[] + lst_idiag2 = [] + lst_eye = [] for i in range(n_distributions): - if i==0: - lst_idiag2.append(sps.kron(np.ones((1,n)),sps.eye(n))) + if i == 0: + lst_idiag2.append(sps.kron(np.ones((1, n)), sps.eye(n))) lst_eye.append(-sps.eye(n)) else: - lst_idiag2.append(sps.kron(np.ones((1,n)),sps.eye(n-1,n))) - lst_eye.append(-sps.eye(n-1,n)) - - A_eq2=sps.hstack((sps.block_diag(lst_idiag2),sps.vstack(lst_eye))) - b_eq2=np.zeros((A_eq2.shape[0])) - + lst_idiag2.append(sps.kron(np.ones((1, n)), sps.eye(n - 1, n))) + lst_eye.append(-sps.eye(n - 1, n)) + + A_eq2 = sps.hstack((sps.block_diag(lst_idiag2), sps.vstack(lst_eye))) + b_eq2 = np.zeros((A_eq2.shape[0])) + # full problem - A_eq=sps.vstack((A_eq1,A_eq2)) - b_eq=np.concatenate((b_eq1,b_eq2)) - - if not cvxopt or solver in ['interior-point']: # cvxopt not installed or simplex/interior point - + A_eq = sps.vstack((A_eq1, A_eq2)) + b_eq = np.concatenate((b_eq1, b_eq2)) + + if not cvxopt or solver in ['interior-point']: # cvxopt not installed or simplex/interior point + if solver is None: - solver='interior-point' - - options={'sparse':True,'disp': verbose} - sol=sp.optimize.linprog(c,A_eq=A_eq,b_eq=b_eq,method=solver,options=options) - x=sol.x - b=x[-n:] - + solver = 'interior-point' + + options = {'sparse': True, 'disp': verbose} + sol = sp.optimize.linprog(c, A_eq=A_eq, b_eq=b_eq, method=solver, options=options) + x = sol.x + b = x[-n:] + else: - - h=np.zeros((n_distributions*n2+n)) - G=-sps.eye(n_distributions*n2+n) - - sol=solvers.lp(matrix(c),scipy_sparse_to_spmatrix(G),matrix(h),A=scipy_sparse_to_spmatrix(A_eq),b=matrix(b_eq),solver=solver) - - x=np.array(sol['x']) - b=x[-n:].ravel() - + + h = np.zeros((n_distributions * n2 + n)) + G = -sps.eye(n_distributions * n2 + n) + + sol = solvers.lp(matrix(c), scipy_sparse_to_spmatrix(G), matrix(h), A=scipy_sparse_to_spmatrix(A_eq), b=matrix(b_eq), solver=solver) + + x = np.array(sol['x']) + b = x[-n:].ravel() + if log: return b, sol else: - return b \ No newline at end of file + return b -- cgit v1.2.3 From 4285cf64f8a2ec481586a190dd545d2a8946e134 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Fri, 11 May 2018 17:01:23 +0200 Subject: remove unused sparse --- ot/lp/cvx.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) (limited to 'ot') diff --git a/ot/lp/cvx.py b/ot/lp/cvx.py index 93097d1..193c0f5 100644 --- a/ot/lp/cvx.py +++ b/ot/lp/cvx.py @@ -13,7 +13,7 @@ import scipy.sparse as sps try: import cvxopt - from cvxopt import solvers, matrix, sparse, spmatrix + from cvxopt import solvers, matrix, spmatrix except ImportError: cvxopt = False @@ -114,13 +114,15 @@ def barycenter(A, M, weights=None, verbose=False, log=False, solver='interior-po A_eq = sps.vstack((A_eq1, A_eq2)) b_eq = np.concatenate((b_eq1, b_eq2)) - if not cvxopt or solver in ['interior-point']: # cvxopt not installed or simplex/interior point + if not cvxopt or solver in ['interior-point']: + # cvxopt not installed or interior point if solver is None: solver = 'interior-point' options = {'sparse': True, 'disp': verbose} - sol = sp.optimize.linprog(c, A_eq=A_eq, b_eq=b_eq, method=solver, options=options) + sol = sp.optimize.linprog(c, A_eq=A_eq, b_eq=b_eq, method=solver, + options=options) x = sol.x b = x[-n:] @@ -129,7 +131,9 @@ def barycenter(A, M, weights=None, verbose=False, log=False, solver='interior-po h = np.zeros((n_distributions * n2 + n)) G = -sps.eye(n_distributions * n2 + n) - sol = solvers.lp(matrix(c), scipy_sparse_to_spmatrix(G), matrix(h), A=scipy_sparse_to_spmatrix(A_eq), b=matrix(b_eq), solver=solver) + sol = solvers.lp(matrix(c), scipy_sparse_to_spmatrix(G), matrix(h), + A=scipy_sparse_to_spmatrix(A_eq), b=matrix(b_eq), + solver=solver) x = np.array(sol['x']) b = x[-n:].ravel() -- cgit v1.2.3 From 36f4f7ed2116841d7fe9514ee250bbf16e77b72d Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Fri, 11 May 2018 17:07:56 +0200 Subject: better documentation --- ot/lp/__init__.py | 3 +++ ot/lp/cvx.py | 4 ++-- 2 files changed, 5 insertions(+), 2 deletions(-) (limited to 'ot') diff --git a/ot/lp/__init__.py b/ot/lp/__init__.py index 6371feb..5dda82a 100644 --- a/ot/lp/__init__.py +++ b/ot/lp/__init__.py @@ -11,9 +11,12 @@ import multiprocessing import numpy as np +from .import cvx + # import compiled emd from .emd_wrap import emd_c, check_result from ..utils import parmap +from .cvx import barycenter def emd(a, b, M, numItermax=100000, log=False): diff --git a/ot/lp/cvx.py b/ot/lp/cvx.py index 193c0f5..c62da6a 100644 --- a/ot/lp/cvx.py +++ b/ot/lp/cvx.py @@ -38,8 +38,8 @@ def barycenter(A, M, weights=None, verbose=False, log=False, solver='interior-po - :math:`W_1(\cdot,\cdot)` is the Wasserstein distance (see ot.emd.sinkhorn) - :math:`\mathbf{a}_i` are training distributions in the columns of matrix :math:`\mathbf{A}` - The linear program is solved using the default cvxopt solver if installed. - If cvxopt is not installed it uses the lp solver from scipy.optimize. + The linear program is solved using the interior point solver from scipy.optimize. + If cvxopt solver if installed it can use cvxopt. Parameters ---------- -- cgit v1.2.3 From fdb2f3af19d04872bafa0d9ec5563732e1d6209b Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Fri, 11 May 2018 17:24:09 +0200 Subject: add test for barycenter --- ot/lp/cvx.py | 12 +++++++----- test/test_gpu.py | 2 +- test/test_ot.py | 35 ++++++++++++++++++++++++++++++++++- 3 files changed, 42 insertions(+), 7 deletions(-) (limited to 'ot') diff --git a/ot/lp/cvx.py b/ot/lp/cvx.py index c62da6a..fe9ac76 100644 --- a/ot/lp/cvx.py +++ b/ot/lp/cvx.py @@ -39,7 +39,9 @@ def barycenter(A, M, weights=None, verbose=False, log=False, solver='interior-po - :math:`\mathbf{a}_i` are training distributions in the columns of matrix :math:`\mathbf{A}` The linear program is solved using the interior point solver from scipy.optimize. - If cvxopt solver if installed it can use cvxopt. + If cvxopt solver if installed it can use cvxopt + + Note that this problem do not scale well (both in memory and computational time). Parameters ---------- @@ -114,14 +116,14 @@ def barycenter(A, M, weights=None, verbose=False, log=False, solver='interior-po A_eq = sps.vstack((A_eq1, A_eq2)) b_eq = np.concatenate((b_eq1, b_eq2)) - if not cvxopt or solver in ['interior-point']: + if not cvxopt or solver in ['interior-point']: # cvxopt not installed or interior point if solver is None: solver = 'interior-point' options = {'sparse': True, 'disp': verbose} - sol = sp.optimize.linprog(c, A_eq=A_eq, b_eq=b_eq, method=solver, + sol = sp.optimize.linprog(c, A_eq=A_eq, b_eq=b_eq, method=solver, options=options) x = sol.x b = x[-n:] @@ -131,8 +133,8 @@ def barycenter(A, M, weights=None, verbose=False, log=False, solver='interior-po h = np.zeros((n_distributions * n2 + n)) G = -sps.eye(n_distributions * n2 + n) - sol = solvers.lp(matrix(c), scipy_sparse_to_spmatrix(G), matrix(h), - A=scipy_sparse_to_spmatrix(A_eq), b=matrix(b_eq), + sol = solvers.lp(matrix(c), scipy_sparse_to_spmatrix(G), matrix(h), + A=scipy_sparse_to_spmatrix(A_eq), b=matrix(b_eq), solver=solver) x = np.array(sol['x']) diff --git a/test/test_gpu.py b/test/test_gpu.py index 615c2a7..1e97c45 100644 --- a/test/test_gpu.py +++ b/test/test_gpu.py @@ -76,4 +76,4 @@ def test_gpu_sinkhorn_lpl1(): time3 - time2)) describe_res(G2) - np.testing.assert_allclose(G1, G2, rtol=1e-5, atol=1e-5) + np.testing.assert_allclose(G1, G2, rtol=1e-3, atol=1e-3) diff --git a/test/test_ot.py b/test/test_ot.py index ea6d9dc..bf23e8c 100644 --- a/test/test_ot.py +++ b/test/test_ot.py @@ -10,7 +10,7 @@ import numpy as np import ot from ot.datasets import get_1D_gauss as gauss - +import pytest def test_doctest(): import doctest @@ -117,6 +117,39 @@ def test_emd2_multi(): np.testing.assert_allclose(emd1, emdn) +def test_lp_barycenter(): + + a1 = np.array([1.0, 0, 0])[:, None] + a2 = np.array([0, 0, 1.0])[:, None] + + A = np.hstack((a1, a2)) + M = np.array([[0, 1.0, 4.0], [1.0, 0, 1.0], [4.0, 1.0, 0]]) + + # obvious barycenter between two diracs + bary0 = np.array([0, 1.0, 0]) + + bary = ot.lp.barycenter(A, M, [.5, .5]) + + np.testing.assert_allclose(bary, bary0, rtol=1e-5, atol=1e-7) + np.testing.assert_allclose(bary.sum(), 1) + +@pytest.mark.skipif(not ot.lp.cvx.cvxopt, reason="No cvxopt available") +def test_lp_barycenter_cvxopt(): + + a1 = np.array([1.0, 0, 0])[:, None] + a2 = np.array([0, 0, 1.0])[:, None] + + A = np.hstack((a1, a2)) + M = np.array([[0, 1.0, 4.0], [1.0, 0, 1.0], [4.0, 1.0, 0]]) + + # obvious barycenter between two diracs + bary0 = np.array([0, 1.0, 0]) + + bary = ot.lp.barycenter(A, M, [.5, .5],solver=None) + + np.testing.assert_allclose(bary, bary0, rtol=1e-5, atol=1e-7) + np.testing.assert_allclose(bary.sum(), 1) + def test_warnings(): n = 100 # nb bins m = 100 # nb bins -- cgit v1.2.3 From 54f0b47e55c966d5492e4ce19ec4e704ef3278d6 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Tue, 29 May 2018 16:08:33 +0200 Subject: update documentation for barycenter function --- ot/bregman.py | 4 ++-- ot/lp/cvx.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) (limited to 'ot') diff --git a/ot/bregman.py b/ot/bregman.py index e788ef5..b017c1a 100644 --- a/ot/bregman.py +++ b/ot/bregman.py @@ -839,13 +839,13 @@ def barycenter(A, M, reg, weights=None, numItermax=1000, Parameters ---------- A : np.ndarray (d,n) - n training distributions of size d + n training distributions a_i of size d M : np.ndarray (d,d) loss matrix for OT reg : float Regularization term >0 weights : np.ndarray (n,) - Weights of each histogram i_i on the simplex + Weights of each histogram a_i on the simplex (barycentric coodinates) numItermax : int, optional Max number of iterations stopThr : float, optional diff --git a/ot/lp/cvx.py b/ot/lp/cvx.py index fe9ac76..c8c75bc 100644 --- a/ot/lp/cvx.py +++ b/ot/lp/cvx.py @@ -46,13 +46,13 @@ def barycenter(A, M, weights=None, verbose=False, log=False, solver='interior-po Parameters ---------- A : np.ndarray (d,n) - n training distributions of size d + n training distributions a_i of size d M : np.ndarray (d,d) loss matrix for OT reg : float Regularization term >0 weights : np.ndarray (n,) - Weights of each histogram i_i on the simplex + Weights of each histogram a_i on the simplex (barycentric coodinates) verbose : bool, optional Print information along iterations log : bool, optional -- cgit v1.2.3