From 9f63ee92e281427ab3d520f75bb9c3406b547365 Mon Sep 17 00:00:00 2001 From: Laetitia Chapel Date: Thu, 9 Apr 2020 13:55:27 +0200 Subject: initial commit partial Wass and GW --- ot/partial.py | 1015 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1015 insertions(+) create mode 100755 ot/partial.py (limited to 'ot/partial.py') diff --git a/ot/partial.py b/ot/partial.py new file mode 100755 index 0000000..746f337 --- /dev/null +++ b/ot/partial.py @@ -0,0 +1,1015 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Partial OT +""" + +# Author: Laetitia Chapel +# License: MIT License + +import numpy as np + +from ot.lp import emd + + +def partial_wasserstein_lagrange(a, b, M, reg_m=None, nb_dummies=1, log=False, + **kwargs): + + r""" + Solves the partial optimal transport problem for the quadratic cost + and returns the OT plan + + The function considers the following problem: + + .. math:: + \gamma = \arg\min_\gamma <\gamma,(M-\lambda)>_F + + s.t. + \gamma\geq 0 \\ + \gamma 1 \leq a\\ + \gamma^T 1 \leq b\\ + 1^T \gamma^T 1 = m \leq \min\{\|a\|_1, \|b\|_1\} + + + or equivalently: + + .. math:: + \gamma = \arg\min_\gamma <\gamma,M>_F + \sqrt(\lambda/2) + (\|\gamma 1 - a\|_1 + \|\gamma^T 1 - b\|_1) + + s.t. + \gamma\geq 0 \\ + + + where : + + - M is the metric cost matrix + - a and b are source and target unbalanced distributions + - :math:`\lambda` is the lagragian cost. Tuning its value allows attaining + a given mass to be transported m + + The formulation of the problem has been proposed in [26]_ + + Parameters + ---------- + a : np.ndarray (dim_a,) + Unnormalized histogram of dimension dim_a + b : np.ndarray (dim_b,) + Unnormalized histograms of dimension dim_b + M : np.ndarray (dim_a, dim_b) + cost matrix for the quadratic cost + reg_m : float, optional + Lagragian cost + log : bool, optional + record log if True + + + Returns + ------- + gamma : (dim_a x dim_b) ndarray + Optimal transportation matrix for the given parameters + log : dict + log dictionary returned only if `log` is `True` + + + Examples + -------- + + >>> import ot + >>> a = [.1, .2] + >>> b = [.1, .1] + >>> M = [[0., 1.], [2., 3.]] + >>> np.round(partial_wasserstein_lagrange(a,b,M), 2) + array([[0.1, 0. ], + [0. , 0.1]]) + >>> np.round(partial_wasserstein_lagrange(a,b,M,reg_m=2), 2) + array([[0.1, 0. ], + [0. , 0. ]]) + + References + ---------- + + .. [26] Caffarelli, L. A., & McCann, R. J. (2010) Free boundaries in + optimal transport and Monge-Ampere obstacle problems. Annals of + mathematics, 673-730. + + See Also + -------- + ot.partial.partial_wasserstein : Partial Wasserstein with fixed mass + """ + + if np.sum(a) > 1 or np.sum(b) > 1: + raise ValueError("Problem infeasible. Check that a and b are in the " + "simplex") + + if reg_m is None: + reg_m = np.max(M) + 1 + if reg_m < -np.max(M): + return np.zeros((len(a), len(b))) + + eps = 1e-20 + M = np.asarray(M, dtype=np.float64) + b = np.asarray(b, dtype=np.float64) + a = np.asarray(a, dtype=np.float64) + + M_star = M - reg_m # modified cost matrix + + # trick to fasten the computation: select only the subset of columns/lines + # that can have marginals greater than 0 (that is to say M < 0) + idx_x = np.where(np.min(M_star, axis=1) < eps)[0] + idx_y = np.where(np.min(M_star, axis=0) < eps)[0] + + # extend a, b, M with "reservoir" or "dummy" points + M_extended = np.zeros((len(idx_x) + nb_dummies, len(idx_y) + nb_dummies)) + M_extended[:len(idx_x), :len(idx_y)] = M_star[np.ix_(idx_x, idx_y)] + + a_extended = np.append(a[idx_x], [(np.sum(a) - np.sum(a[idx_x]) + + np.sum(b)) / nb_dummies] * nb_dummies) + b_extended = np.append(b[idx_y], [(np.sum(b) - np.sum(b[idx_y]) + + np.sum(a)) / nb_dummies] * nb_dummies) + + gamma_extended, log_emd = emd(a_extended, b_extended, M_extended, log=True, + **kwargs) + gamma = np.zeros((len(a), len(b))) + gamma[np.ix_(idx_x, idx_y)] = gamma_extended[:-nb_dummies, :-nb_dummies] + + if log_emd['warning'] is not None: + raise ValueError("Error in the EMD resolution: try to increase the" + " number of dummy points") + log_emd['cost'] = np.sum(gamma*M) + if log: + return gamma, log_emd + else: + return gamma + + +def partial_wasserstein(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): + r""" + Solves the partial optimal transport problem for the quadratic cost + and returns the OT plan + + The function considers the following problem: + + .. math:: + \gamma = \arg\min_\gamma <\gamma,M>_F + + s.t. + \gamma\geq 0 \\ + \gamma 1 \leq a\\ + \gamma^T 1 \leq b\\ + 1^T \gamma^T 1 = m \leq \min\{\|a\|_1, \|b\|_1\} + + + where : + + - M is the metric cost matrix + - a and b are source and target unbalanced distributions + - m is the amount of mass to be transported + + Parameters + ---------- + a : np.ndarray (dim_a,) + Unnormalized histogram of dimension dim_a + b : np.ndarray (dim_b,) + Unnormalized histograms of dimension dim_b + M : np.ndarray (dim_a, dim_b) + cost matrix for the quadratic cost + m : float, optional + amount of mass to be transported + nb_dummies : int, optional, default:1 + number of reservoir points to be added (to avoid numerical + instabilities, increase its value if an error is raised) + log : bool, optional + record log if True + + + Returns + ------- + :math:`gamma` : (dim_a x dim_b) ndarray + Optimal transportation matrix for the given parameters + log : dict + log dictionary returned only if `log` is `True` + + + Examples + -------- + + >>> import ot + >>> a = [.1, .2] + >>> b = [.1, .1] + >>> M = [[0., 1.], [2., 3.]] + >>> np.round(partial_wasserstein(a,b,M), 2) + array([[0.1, 0. ], + [0. , 0.1]]) + >>> np.round(partial_wasserstein(a,b,M,m=0.1), 2) + array([[0.1, 0. ], + [0. , 0. ]]) + + References + ---------- + .. [26] Caffarelli, L. A., & McCann, R. J. (2010) Free boundaries in + optimal transport and Monge-Ampere obstacle problems. Annals of + mathematics, 673-730. + .. [27] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + Wasserstein with Applications on Positive-Unlabeled Learning". + arXiv preprint arXiv:2002.08276. + + See Also + -------- + ot.partial.partial_wasserstein_lagrange: Partial Wasserstein with + regularization on the marginals + ot.partial.entropic_partial_wasserstein: Partial Wasserstein with a + entropic regularization parameter + """ + + if m is None: + return partial_wasserstein_lagrange(a, b, M, log=log, **kwargs) + elif m < 0: + raise ValueError("Problem infeasible. Parameter m should be greater" + " than 0.") + elif m > np.min((np.sum(a), np.sum(b))): + raise ValueError("Problem infeasible. Parameter m should lower or" + " equal than min(|a|_1, |b|_1).") + + b_extended = np.append(b, [(np.sum(a) - m) / nb_dummies] * nb_dummies) + a_extended = np.append(a, [(np.sum(b) - m) / nb_dummies] * nb_dummies) + M_extended = np.ones((len(a_extended), len(b_extended))) * np.max(M) * 1e2 + M_extended[-1, -1] = np.max(M) * 1e5 + M_extended[:len(a), :len(b)] = M + + gamma, log_emd = emd(a_extended, b_extended, M_extended, log=True, + **kwargs) + if log_emd['warning'] is not None: + raise ValueError("Error in the EMD resolution: try to increase the" + " number of dummy points") + log_emd['partial_w_dist'] = np.sum(M * gamma[:len(a), :len(b)]) + + if log: + return gamma[:len(a), :len(b)], log_emd + else: + return gamma[:len(a), :len(b)] + + +def partial_wasserstein2(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): + r""" + Solves the partial optimal transport problem for the quadratic cost + and returns the partial GW discrepancy + + The function considers the following problem: + + .. math:: + \gamma = \arg\min_\gamma <\gamma,M>_F + + s.t. + \gamma\geq 0 \\ + \gamma 1 \leq a\\ + \gamma^T 1 \leq b\\ + 1^T \gamma^T 1 = m \leq \min\{\|a\|_1, \|b\|_1\} + + + where : + + - M is the metric cost matrix + - a and b are source and target unbalanced distributions + - m is the amount of mass to be transported + + Parameters + ---------- + a : np.ndarray (dim_a,) + Unnormalized histogram of dimension dim_a + b : np.ndarray (dim_b,) + Unnormalized histograms of dimension dim_b + M : np.ndarray (dim_a, dim_b) + cost matrix for the quadratic cost + m : float, optional + amount of mass to be transported + nb_dummies : int, optional, default:1 + number of reservoir points to be added (to avoid numerical + instabilities, increase its value if an error is raised) + log : bool, optional + record log if True + + + Returns + ------- + :math:`gamma` : (dim_a x dim_b) ndarray + Optimal transportation matrix for the given parameters + log : dict + log dictionary returned only if `log` is `True` + + + Examples + -------- + + >>> import ot + >>> a=[.1, .2] + >>> b=[.1, .1] + >>> M=[[0., 1.], [2., 3.]] + >>> np.round(partial_wasserstein2(a, b, M), 1) + 0.3 + >>> np.round(partial_wasserstein2(a,b,M,m=0.1), 1) + 0.0 + + References + ---------- + .. [26] Caffarelli, L. A., & McCann, R. J. (2010) Free boundaries in + optimal transport and Monge-Ampere obstacle problems. Annals of + mathematics, 673-730. + .. [27] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + Wasserstein with Applications on Positive-Unlabeled Learning". + arXiv preprint arXiv:2002.08276. + """ + + partial_gw, log_w = partial_wasserstein(a, b, M, m, nb_dummies, log=True, + **kwargs) + + log_w['T'] = partial_gw + + if log: + return np.sum(partial_gw * M), log_w + else: + return np.sum(partial_gw * M) + + +def gwgrad_partial(C1, C2, T): + """Compute the GW gradient. Note: we can not use the trick in [12]_ as + the marginals may not sum to 1. + + Parameters + ---------- + C1: array of shape (n_p,n_p) + intra-source (P) cost matrix + + C2: array of shape (n_u,n_u) + intra-target (U) cost matrix + + T : array of shape(n_p+nb_dummies, n_u) (default: None) + Transport matrix + + Returns + ------- + numpy.array of shape (n_p+nb_dummies, n_u) + gradient + + References + ---------- + .. [12] Peyré, Gabriel, Marco Cuturi, and Justin Solomon, + "Gromov-Wasserstein averaging of kernel and distance matrices." + International Conference on Machine Learning (ICML). 2016. + """ + cC1 = np.dot(C1 ** 2 / 2, np.dot(T, np.ones(C2.shape[0]).reshape(-1, 1))) + cC2 = np.dot(np.dot(np.ones(C1.shape[0]).reshape(1, -1), T), C2 ** 2 / 2) + constC = cC1 + cC2 + A = -np.dot(C1, T).dot(C2.T) + tens = constC + A + return tens * 2 + + +def gwloss_partial(C1, C2, T): + """Compute the GW loss. + + Parameters + ---------- + C1: array of shape (n_p,n_p) + intra-source (P) cost matrix + + C2: array of shape (n_u,n_u) + intra-target (U) cost matrix + + T : array of shape(n_p+nb_dummies, n_u) (default: None) + Transport matrix + + Returns + ------- + GW loss + """ + g = gwgrad_partial(C1, C2, T) * 0.5 + return np.sum(g * T) + + +def partial_gromov_wasserstein(C1, C2, p, q, m=None, nb_dummies=1, G0=None, + thres=1, numItermax=1000, tol=1e-7, + log=False, verbose=False, **kwargs): + r""" + Solves the partial optimal transport problem + and returns the OT plan + + The function considers the following problem: + + .. math:: + \gamma = arg\min_\gamma <\gamma,M>_F + + s.t. \gamma 1 \leq a \\ + \gamma^T 1 \leq b \\ + \gamma\geq 0 \\ + 1^T \gamma^T 1 = m \leq \min\{\|a\|_1, \|b\|_1\} \\ + + where : + + - M is the metric cost matrix + - :math:`\Omega` is the entropic regularization term :math:`\Omega(\gamma) + =\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` + - a and b are the sample weights + - m is the amount of mass to be transported + + The formulation of the problem has been proposed in [27]_ + + + Parameters + ---------- + C1 : ndarray, shape (ns, ns) + Metric cost matrix in the source space + C2 : ndarray, shape (nt, nt) + Metric costfr matrix in the target space + p : ndarray, shape (ns,) + Distribution in the source space + q : ndarray, shape (nt,) + Distribution in the target space + m : float, optional + Amount of mass to be transported (default: min (|p|_1, |q|_1)) + nb_dummies : int, optional + Number of dummy points to add (avoid instabilities in the EMD solver) + G0 : ndarray, shape (ns, nt), optional + Initialisation of the transportation matrix + thres : float, optional + quantile of the gradient matrix to populate the cost matrix when 0 + (default: 1) + numItermax : int, optional + Max number of iterations + log : bool, optional + return log if True + verbose : bool, optional + Print information along iterations + armijo : bool, optional + If True the steps of the line-search is found via an armijo research. Else closed form is used. + If there is convergence issues use False. + **kwargs : dict + parameters can be directly passed to the emd solver + + + Returns + ------- + gamma : (dim_a x dim_b) ndarray + Optimal transportation matrix for the given parameters + log : dict + log dictionary returned only if `log` is `True` + + + Examples + -------- + >>> import ot + >>> import scipy as sp + >>> a = np.array([0.25] * 4) + >>> b = np.array([0.25] * 4) + >>> x = np.array([1,2,100,200]).reshape((-1,1)) + >>> y = np.array([3,2,98,199]).reshape((-1,1)) + >>> C1 = sp.spatial.distance.cdist(x, x) + >>> C2 = sp.spatial.distance.cdist(y, y) + >>> np.round(partial_gromov_wasserstein(C1, C2, a, b),2) + array([[0. , 0.25, 0. , 0. ], + [0.25, 0. , 0. , 0. ], + [0. , 0. , 0.25, 0. ], + [0. , 0. , 0. , 0.25]]) + >>> np.round(partial_gromov_wasserstein(C1, C2, a, b, m=0.25),2) + array([[0. , 0. , 0. , 0. ], + [0. , 0. , 0. , 0. ], + [0. , 0. , 0. , 0. ], + [0. , 0. , 0. , 0.25]]) + + References + ---------- + .. [27] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + Wasserstein with Applications on Positive-Unlabeled Learning". + arXiv preprint arXiv:2002.08276. + + """ + + if m is None: + m = np.min((np.sum(p), np.sum(q))) + elif m < 0: + raise ValueError("Problem infeasible. Parameter m should be greater" + " than 0.") + elif m > np.min((np.sum(p), np.sum(q))): + raise ValueError("Problem infeasible. Parameter m should lower or" + " equal than min(|a|_1, |b|_1).") + + if G0 is None: + G0 = np.outer(p, q) + + dim_G_extended = (len(p) + nb_dummies, len(q) + nb_dummies) + q_extended = np.append(q, [(np.sum(p) - m) / nb_dummies] * nb_dummies) + p_extended = np.append(p, [(np.sum(q) - m) / nb_dummies] * nb_dummies) + + cpt = 0 + err = 1 + eps = 1e-20 + if log: + log = {'err': []} + + while (err > tol and cpt < numItermax): + + Gprev = G0 + + M = gwgrad_partial(C1, C2, G0) + M[M < eps] = np.quantile(M[M > eps], thres) + + M_emd = np.ones(dim_G_extended) * np.max(M) * 1e2 + M_emd[:len(p), :len(q)] = M + M_emd[-nb_dummies:, -nb_dummies:] = np.max(M) * 1e5 + M_emd = np.asarray(M_emd, dtype=np.float64) + + Gc, logemd = emd(p_extended, q_extended, M_emd, log=True, **kwargs) + + if logemd['warning'] is not None: + raise ValueError("Error in the EMD resolution: try to increase the" + " number of dummy points") + + G0 = Gc[:len(p), :len(q)] + + if cpt % 10 == 0: # to speed up the computations + err = np.linalg.norm(G0 - Gprev) + if log: + log['err'].append(err) + if verbose: + if cpt % 200 == 0: + print('{:5s}|{:12s}|{:12s}'.format( + 'It.', 'Err', 'Loss') + '\n' + '-' * 31) + print('{:5d}|{:8e}|{:8e}'.format(cpt, err, + gwloss_partial(C1, C2, G0))) + + cpt += 1 + + if log: + log['partial_gw_dist'] = gwloss_partial(C1, C2, G0) + return G0[:len(p), :len(q)], log + else: + return G0[:len(p), :len(q)] + + +def partial_gromov_wasserstein2(C1, C2, p, q, m=None, nb_dummies=1, G0=None, + thres=0.75, numItermax=1000, tol=1e-7, + log=False, verbose=False, **kwargs): + r""" + Solves the partial optimal transport problem + and returns the partial Gromov-Wasserstein discrepancy + + The function considers the following problem: + + .. math:: + \gamma = arg\min_\gamma <\gamma,M>_F + + s.t. \gamma 1 \leq a \\ + \gamma^T 1 \leq b \\ + \gamma\geq 0 \\ + 1^T \gamma^T 1 = m \leq \min\{\|a\|_1, \|b\|_1\} \\ + + where : + + - M is the metric cost matrix + - :math:`\Omega` is the entropic regularization term + :math:`\Omega=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` + - a and b are the sample weights + - m is the amount of mass to be transported + + The formulation of the problem has been proposed in [27]_ + + + Parameters + ---------- + C1 : ndarray, shape (ns, ns) + Metric cost matrix in the source space + C2 : ndarray, shape (nt, nt) + Metric costfr matrix in the target space + p : ndarray, shape (ns,) + Distribution in the source space + q : ndarray, shape (nt,) + Distribution in the target space + m : float, optional + Amount of mass to be transported (default: min (|p|_1, |q|_1)) + nb_dummies : int, optional + Number of dummy points to add (avoid instabilities in the EMD solver) + G0 : ndarray, shape (ns, nt), optional + Initialisation of the transportation matrix + thres : float, optional + quantile of the gradient matrix to populate the cost matrix when 0 + (default: 1) + numItermax : int, optional + Max number of iterations + log : bool, optional + return log if True + verbose : bool, optional + Print information along iterations + **kwargs : dict + parameters can be directly passed to the emd solver + + + Returns + ------- + partial_gw_dist : (dim_a x dim_b) ndarray + partial GW discrepancy + log : dict + log dictionary returned only if `log` is `True` + + + Examples + -------- + >>> import ot + >>> import scipy as sp + >>> a = np.array([0.25] * 4) + >>> b = np.array([0.25] * 4) + >>> x = np.array([1,2,100,200]).reshape((-1,1)) + >>> y = np.array([3,2,98,199]).reshape((-1,1)) + >>> C1 = sp.spatial.distance.cdist(x, x) + >>> C2 = sp.spatial.distance.cdist(y, y) + >>> np.round(partial_gromov_wasserstein2(C1, C2, a, b),2) + 1.69 + >>> np.round(partial_gromov_wasserstein2(C1, C2, a, b, m=0.25),2) + 0.0 + + References + ---------- + .. [27] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + Wasserstein with Applications on Positive-Unlabeled Learning". + arXiv preprint arXiv:2002.08276. + + """ + + partial_gw, log_gw = partial_gromov_wasserstein(C1, C2, p, q, m, + nb_dummies, G0, thres, + numItermax, tol, True, + verbose, **kwargs) + + log_gw['T'] = partial_gw + + if log: + return log_gw['partial_gw_dist'], log_gw + else: + return log_gw['partial_gw_dist'] + + +def entropic_partial_wasserstein(a, b, M, reg, m=None, numItermax=1000, + stopThr=1e-100, verbose=False, log=False): + r""" + Solves the partial optimal transport problem + and returns the OT plan + + The function considers the following problem: + + .. math:: + \gamma = arg\min_\gamma <\gamma,M>_F + reg\cdot\Omega(\gamma) + + s.t. \gamma 1 \leq a \\ + \gamma^T 1 \leq b \\ + \gamma\geq 0 \\ + 1^T \gamma^T 1 = m \leq \min\{\|a\|_1, \|b\|_1\} \\ + + where : + + - M is the metric cost matrix + - :math:`\Omega` is the entropic regularization term + :math:`\Omega=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` + - a and b are the sample weights + - m is the amount of mass to be transported + + The formulation of the problem has been proposed in [3]_ + + + Parameters + ---------- + a : np.ndarray (dim_a,) + Unnormalized histogram of dimension dim_a + b : np.ndarray (dim_b,) + Unnormalized histograms of dimension dim_b + M : np.ndarray (dim_a, dim_b) + cost matrix + reg : float + Regularization term > 0 + m : float, optional + Amount of mass to be transported + numItermax : int, optional + Max number of iterations + stopThr : float, optional + Stop threshold on error (>0) + verbose : bool, optional + Print information along iterations + log : bool, optional + record log if True + + + Returns + ------- + gamma : (dim_a x dim_b) ndarray + Optimal transportation matrix for the given parameters + log : dict + log dictionary returned only if `log` is `True` + + + Examples + -------- + >>> import ot + >>> a = [.1, .2] + >>> b = [.1, .1] + >>> M = [[0., 1.], [2., 3.]] + >>> np.round(entropic_partial_wasserstein(a, b, M, 1, 0.1), 2) + array([[0.06, 0.02], + [0.01, 0. ]]) + + + References + ---------- + .. [3] Benamou, J. D., Carlier, G., Cuturi, M., Nenna, L., & Peyré, G. + (2015). Iterative Bregman projections for regularized transportation + problems. SIAM Journal on Scientific Computing, 37(2), A1111-A1138. + + See Also + -------- + ot.partial.partial_wasserstein: exact Partial Wasserstein + """ + + a = np.asarray(a, dtype=np.float64) + b = np.asarray(b, dtype=np.float64) + M = np.asarray(M, dtype=np.float64) + + dim_a, dim_b = M.shape + dx = np.ones(dim_a) + dy = np.ones(dim_b) + + if len(a) == 0: + a = np.ones(dim_a, dtype=np.float64) / dim_a + if len(b) == 0: + b = np.ones(dim_b, dtype=np.float64) / dim_b + + if m is None: + m = np.min((np.sum(a), np.sum(b))) + if m < 0: + raise ValueError("Problem infeasible. Parameter m should be greater" + " than 0.") + if m > np.min((np.sum(a), np.sum(b))): + raise ValueError("Problem infeasible. Parameter m should lower or" + " equal than min(|a|_1, |b|_1).") + + log_e = {'err': []} + + # Next 3 lines equivalent to K=np.exp(-M/reg), but faster to compute + K = np.empty(M.shape, dtype=M.dtype) + np.divide(M, -reg, out=K) + np.exp(K, out=K) + np.multiply(K, m/np.sum(K), out=K) + + err, cpt = 1, 0 + + while (err > stopThr and cpt < numItermax): + Kprev = K + K1 = np.dot(np.diag(np.minimum(a / np.sum(K, axis=1), dx)), K) + K2 = np.dot(K1, np.diag(np.minimum(b / np.sum(K1, axis=0), dy))) + K = K2 * (m / np.sum(K2)) + + if np.any(np.isnan(K)) or np.any(np.isinf(K)): + print('Warning: numerical errors at iteration', cpt) + break + if cpt % 10 == 0: + err = np.linalg.norm(Kprev - K) + if log: + log_e['err'].append(err) + if verbose: + if cpt % 200 == 0: + print( + '{:5s}|{:12s}'.format('It.', 'Err') + '\n' + '-' * 11) + print('{:5d}|{:8e}|'.format(cpt, err)) + + cpt = cpt + 1 + log_e['partial_w_dist'] = np.sum(M * K) + if log: + return K, log_e + else: + return K + + +def entropic_partial_gromov_wasserstein(C1, C2, p, q, reg, m=None, G0=None, + numItermax=1000, tol=1e-7, log=False, + verbose=False): + r""" + Returns the partial Gromov-Wasserstein transport between (C1,p) and (C2,q) + + The function solves the following optimization problem: + + .. math:: + GW = \arg\min_{\gamma} \sum_{i,j,k,l} L(C1_{i,k},C2_{j,l})\cdot + \gamma_{i,j}\cdot\gamma_{k,l} + reg\cdot\Omega(\gamma) + + s.t. + \gamma\geq 0 \\ + \gamma 1 \leq a\\ + \gamma^T 1 \leq b\\ + 1^T \gamma^T 1 = m \leq \min\{\|a\|_1, \|b\|_1\} + + where : + + - C1 is the metric cost matrix in the source space + - C2 is the metric cost matrix in the target space + - p and q are the sample weights + - L : quadratic loss function + - :math:`\Omega` is the entropic regularization term + :math:`\Omega=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` + - m is the amount of mass to be transported + + The formulation of the problem has been proposed in [12]. + + Parameters + ---------- + C1 : ndarray, shape (ns, ns) + Metric cost matrix in the source space + C2 : ndarray, shape (nt, nt) + Metric costfr matrix in the target space + p : ndarray, shape (ns,) + Distribution in the source space + q : ndarray, shape (nt,) + Distribution in the target space + reg: float + entropic regularization parameter + m : float, optional + Amount of mass to be transported (default: min (|p|_1, |q|_1)) + G0 : ndarray, shape (ns, nt), optional + Initialisation of the transportation matrix + numItermax : int, optional + Max number of iterations + log : bool, optional + return log if True + verbose : bool, optional + Print information along iterations + + Examples + -------- + >>> import ot + >>> import scipy as sp + >>> a = np.array([0.25] * 4) + >>> b = np.array([0.25] * 4) + >>> x = np.array([1,2,100,200]).reshape((-1,1)) + >>> y = np.array([3,2,98,199]).reshape((-1,1)) + >>> C1 = sp.spatial.distance.cdist(x, x) + >>> C2 = sp.spatial.distance.cdist(y, y) + >>> np.round(entropic_partial_gromov_wasserstein(C1, C2, a, b,50), 2) + array([[0.12, 0.13, 0. , 0. ], + [0.13, 0.12, 0. , 0. ], + [0. , 0. , 0.25, 0. ], + [0. , 0. , 0. , 0.25]]) + >>> np.round(entropic_partial_gromov_wasserstein(C1, C2, a, b, 50, m=0.25) + , 2) + array([[0.02, 0.03, 0. , 0.03], + [0.03, 0.03, 0. , 0.03], + [0. , 0. , 0.03, 0. ], + [0.02, 0.02, 0. , 0.03]]) + + Returns + ------- + :math: `gamma` : (dim_a x dim_b) ndarray + Optimal transportation matrix for the given parameters + log : dict + log dictionary returned only if `log` is `True` + + References + ---------- + .. [12] Peyré, Gabriel, Marco Cuturi, and Justin Solomon, + "Gromov-Wasserstein averaging of kernel and distance matrices." + International Conference on Machine Learning (ICML). 2016. + + See Also + -------- + ot.partial.partial_gromov_wasserstein: exact Partial Gromov-Wasserstein + """ + + if G0 is None: + G0 = np.outer(p, q) + + if m is None: + m = np.min((np.sum(p), np.sum(q))) + elif m < 0: + raise ValueError("Problem infeasible. Parameter m should be greater" + " than 0.") + elif m > np.min((np.sum(p), np.sum(q))): + raise ValueError("Problem infeasible. Parameter m should lower or" + " equal than min(|a|_1, |b|_1).") + + cpt = 0 + err = 1 + + loge = {'err': []} + + while (err > tol and cpt < numItermax): + Gprev = G0 + M_entr = gwgrad_partial(C1, C2, G0) + G0 = entropic_partial_wasserstein(p, q, M_entr, reg, m) + if cpt % 10 == 0: # to speed up the computations + err = np.linalg.norm(G0 - Gprev) + if log: + loge['err'].append(err) + if verbose: + if cpt % 200 == 0: + print('{:5s}|{:12s}|{:12s}'.format( + 'It.', 'Err', 'Loss') + '\n' + '-' * 31) + print('{:5d}|{:8e}|{:8e}'.format(cpt, err, + gwloss_partial(C1, C2, G0))) + + cpt += 1 + + if log: + loge['partial_gw_dist'] = gwloss_partial(C1, C2, G0) + return G0, loge + else: + return G0 + + +def entropic_partial_gromov_wasserstein2(C1, C2, p, q, reg, m=None, G0=None, + numItermax=1000, tol=1e-7, log=False, + verbose=False): + r""" + Returns the partial Gromov-Wasserstein discrepancy between (C1,p) and + (C2,q) + + The function solves the following optimization problem: + + .. math:: + GW = \arg\min_{\gamma} \sum_{i,j,k,l} L(C1_{i,k},C2_{j,l})\cdot + \gamma_{i,j}\cdot\gamma_{k,l} + reg\cdot\Omega(\gamma) + + s.t. + \gamma\geq 0 \\ + \gamma 1 \leq a\\ + \gamma^T 1 \leq b\\ + 1^T \gamma^T 1 = m \leq \min\{\|a\|_1, \|b\|_1\} + + where : + + - C1 is the metric cost matrix in the source space + - C2 is the metric cost matrix in the target space + - p and q are the sample weights + - L : quadratic loss function + - :math:`\Omega` is the entropic regularization term + :math:`\Omega=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` + - m is the amount of mass to be transported + + The formulation of the problem has been proposed in [12]. + + + Parameters + ---------- + C1 : ndarray, shape (ns, ns) + Metric cost matrix in the source space + C2 : ndarray, shape (nt, nt) + Metric costfr matrix in the target space + p : ndarray, shape (ns,) + Distribution in the source space + q : ndarray, shape (nt,) + Distribution in the target space + reg: float + entropic regularization parameter + m : float, optional + Amount of mass to be transported (default: min (|p|_1, |q|_1)) + G0 : ndarray, shape (ns, nt), optional + Initialisation of the transportation matrix + numItermax : int, optional + Max number of iterations + log : bool, optional + return log if True + verbose : bool, optional + Print information along iterations + + + Returns + ------- + partial_gw_dist: float + Gromov-Wasserstein distance + log : dict + log dictionary returned only if `log` is `True` + + Examples + -------- + >>> import ot + >>> import scipy as sp + >>> a = np.array([0.25] * 4) + >>> b = np.array([0.25] * 4) + >>> x = np.array([1,2,100,200]).reshape((-1,1)) + >>> y = np.array([3,2,98,199]).reshape((-1,1)) + >>> C1 = sp.spatial.distance.cdist(x, x) + >>> C2 = sp.spatial.distance.cdist(y, y) + >>> np.round(entropic_partial_gromov_wasserstein2(C1, C2, a, b,50), 2) + 1.87 + + References + ---------- + .. [12] Peyré, Gabriel, Marco Cuturi, and Justin Solomon, + "Gromov-Wasserstein averaging of kernel and distance matrices." + International Conference on Machine Learning (ICML). 2016. + """ + + partial_gw, log_gw = entropic_partial_gromov_wasserstein(C1, C2, p, q, reg, + m, G0, numItermax, + tol, True, + verbose) + + log_gw['T'] = partial_gw + + if log: + return log_gw['partial_gw_dist'], log_gw + else: + return log_gw['partial_gw_dist'] -- cgit v1.2.3 From 8c724ad3579959e9d369c0b7fbaa22ea19ced614 Mon Sep 17 00:00:00 2001 From: Laetitia Chapel Date: Wed, 15 Apr 2020 15:32:56 +0200 Subject: partial with tests --- examples/plot_partial_wass_and_gromov.py | 18 +++---- ot/__init__.py | 81 -------------------------------- ot/partial.py | 23 +++++---- ot/unbalanced.py | 66 ++++++++------------------ 4 files changed, 39 insertions(+), 149 deletions(-) delete mode 100644 ot/__init__.py (limited to 'ot/partial.py') diff --git a/examples/plot_partial_wass_and_gromov.py b/examples/plot_partial_wass_and_gromov.py index 2ddeb68..30b3fc0 100755 --- a/examples/plot_partial_wass_and_gromov.py +++ b/examples/plot_partial_wass_and_gromov.py @@ -33,9 +33,9 @@ mu = np.array([0, 0]) cov = np.array([[1, 0], [0, 2]]) xs = ot.datasets.make_2D_samples_gauss(n_samples, mu, cov) -xs = np.append(xs, (np.random.rand(n_noise, 2)+1)*4).reshape((-1, 2)) +xs = np.append(xs, (np.random.rand(n_noise, 2) + 1) * 4).reshape((-1, 2)) xt = ot.datasets.make_2D_samples_gauss(n_samples, mu, cov) -xt = np.append(xt, (np.random.rand(n_noise, 2)+1)*-3).reshape((-1, 2)) +xt = np.append(xt, (np.random.rand(n_noise, 2) + 1) * -3).reshape((-1, 2)) M = sp.spatial.distance.cdist(xs, xt) @@ -62,7 +62,7 @@ w, log = ot.partial.entropic_partial_wasserstein(p, q, M, reg=0.1, m=0.5, log=True) print('Partial Wasserstein distance (m = 0.5): ' + str(log0['partial_w_dist'])) -print('Entropic partial Wasserstein distance (m = 0.5): ' + \ +print('Entropic partial Wasserstein distance (m = 0.5): ' + str(log['partial_w_dist'])) pl.figure(1, (10, 5)) @@ -98,10 +98,10 @@ cov_t = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) xs = ot.datasets.make_2D_samples_gauss(n_samples, mu_s, cov_s) -xs = np.concatenate((xs, ((np.random.rand(n_noise, 2)+1)*4)), axis=0) +xs = np.concatenate((xs, ((np.random.rand(n_noise, 2) + 1) * 4)), axis=0) P = sp.linalg.sqrtm(cov_t) xt = np.random.randn(n_samples, 3).dot(P) + mu_t -xt = np.concatenate((xt, ((np.random.rand(n_noise, 3)+1)*10)), axis=0) +xt = np.concatenate((xt, ((np.random.rand(n_noise, 3) + 1) * 10)), axis=0) fig = pl.figure() ax1 = fig.add_subplot(121) @@ -128,7 +128,7 @@ res, log = ot.partial.entropic_partial_gromov_wasserstein(C1, C2, p, q, 10, m=m, log=True) print('Partial Wasserstein distance (m = 1): ' + str(log0['partial_gw_dist'])) -print('Entropic partial Wasserstein distance (m = 1): ' + \ +print('Entropic partial Wasserstein distance (m = 1): ' + str(log['partial_gw_dist'])) pl.figure(1, (10, 5)) @@ -142,14 +142,14 @@ pl.title('Entropic partial Wasserstein') pl.show() print('-----m = 2/3') -m = 2/3 +m = 2 / 3 res0, log0 = ot.partial.partial_gromov_wasserstein(C1, C2, p, q, m=m, log=True) res, log = ot.partial.entropic_partial_gromov_wasserstein(C1, C2, p, q, 10, m=m, log=True) -print('Partial Wasserstein distance (m = 2/3): ' + \ +print('Partial Wasserstein distance (m = 2/3): ' + str(log0['partial_gw_dist'])) -print('Entropic partial Wasserstein distance (m = 2/3): ' + \ +print('Entropic partial Wasserstein distance (m = 2/3): ' + str(log['partial_gw_dist'])) pl.figure(1, (10, 5)) diff --git a/ot/__init__.py b/ot/__init__.py deleted file mode 100644 index 89c7936..0000000 --- a/ot/__init__.py +++ /dev/null @@ -1,81 +0,0 @@ -""" - -This is the main module of the POT toolbox. It provides easy access to -a number of sub-modules and functions described below. - -.. note:: - - - Here is a list of the submodules and short description of what they contain. - - - :any:`ot.lp` contains OT solvers for the exact (Linear Program) OT problems. - - :any:`ot.bregman` contains OT solvers for the entropic OT problems using - Bregman projections. - - :any:`ot.lp` contains OT solvers for the exact (Linear Program) OT problems. - - :any:`ot.smooth` contains OT solvers for the regularized (l2 and kl) smooth OT - problems. - - :any:`ot.gromov` contains solvers for Gromov-Wasserstein and Fused Gromov - Wasserstein problems. - - :any:`ot.optim` contains generic solvers OT based optimization problems - - :any:`ot.da` contains classes and function related to Monge mapping - estimation and Domain Adaptation (DA). - - :any:`ot.gpu` contains GPU (cupy) implementation of some OT solvers - - :any:`ot.dr` contains Dimension Reduction (DR) methods such as Wasserstein - Discriminant Analysis. - - :any:`ot.utils` contains utility functions such as distance computation and - timing. - - :any:`ot.datasets` contains toy dataset generation functions. - - :any:`ot.plot` contains visualization functions - - :any:`ot.stochastic` contains stochastic solvers for regularized OT. - - :any:`ot.unbalanced` contains solvers for regularized unbalanced OT. - -.. warning:: - The list of automatically imported sub-modules is as follows: - :py:mod:`ot.lp`, :py:mod:`ot.bregman`, :py:mod:`ot.optim` - :py:mod:`ot.utils`, :py:mod:`ot.datasets`, - :py:mod:`ot.gromov`, :py:mod:`ot.smooth` - :py:mod:`ot.stochastic` - - The following sub-modules are not imported due to additional dependencies: - - - :any:`ot.dr` : depends on :code:`pymanopt` and :code:`autograd`. - - :any:`ot.gpu` : depends on :code:`cupy` and a CUDA GPU. - - :any:`ot.plot` : depends on :code:`matplotlib` - -""" - -# Author: Remi Flamary -# Nicolas Courty -# -# License: MIT License - - -# All submodules and packages -from . import lp -from . import bregman -from . import optim -from . import utils -from . import datasets -from . import da -from . import gromov -from . import smooth -from . import stochastic -from . import unbalanced - -# OT functions -from .lp import emd, emd2, emd_1d, emd2_1d, wasserstein_1d -from .bregman import sinkhorn, sinkhorn2, barycenter -from .unbalanced import sinkhorn_unbalanced, barycenter_unbalanced, sinkhorn_unbalanced2 -from .da import sinkhorn_lpl1_mm - -# utils functions -from .utils import dist, unif, tic, toc, toq - -__version__ = "0.6.0" - -__all__ = ['emd', 'emd2', 'emd_1d', 'sinkhorn', 'sinkhorn2', 'utils', 'datasets', - 'bregman', 'lp', 'tic', 'toc', 'toq', 'gromov', - 'emd_1d', 'emd2_1d', 'wasserstein_1d', - 'dist', 'unif', 'barycenter', 'sinkhorn_lpl1_mm', 'da', 'optim', - 'sinkhorn_unbalanced', 'barycenter_unbalanced', - 'sinkhorn_unbalanced2'] diff --git a/ot/partial.py b/ot/partial.py index 746f337..3425acb 100755 --- a/ot/partial.py +++ b/ot/partial.py @@ -9,12 +9,11 @@ Partial OT import numpy as np -from ot.lp import emd +from .lp import emd def partial_wasserstein_lagrange(a, b, M, reg_m=None, nb_dummies=1, log=False, **kwargs): - r""" Solves the partial optimal transport problem for the quadratic cost and returns the OT plan @@ -136,7 +135,7 @@ def partial_wasserstein_lagrange(a, b, M, reg_m=None, nb_dummies=1, log=False, if log_emd['warning'] is not None: raise ValueError("Error in the EMD resolution: try to increase the" " number of dummy points") - log_emd['cost'] = np.sum(gamma*M) + log_emd['cost'] = np.sum(gamma * M) if log: return gamma, log_emd else: @@ -233,7 +232,7 @@ def partial_wasserstein(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): b_extended = np.append(b, [(np.sum(a) - m) / nb_dummies] * nb_dummies) a_extended = np.append(a, [(np.sum(b) - m) / nb_dummies] * nb_dummies) - M_extended = np.ones((len(a_extended), len(b_extended))) * np.max(M) * 1e2 + M_extended = np.ones((len(a_extended), len(b_extended))) * 0 M_extended[-1, -1] = np.max(M) * 1e5 M_extended[:len(a), :len(b)] = M @@ -381,7 +380,7 @@ def gwloss_partial(C1, C2, T): Returns ------- - GW loss + GW loss """ g = gwgrad_partial(C1, C2, T) * 0.5 return np.sum(g * T) @@ -432,7 +431,7 @@ def partial_gromov_wasserstein(C1, C2, p, q, m=None, nb_dummies=1, G0=None, G0 : ndarray, shape (ns, nt), optional Initialisation of the transportation matrix thres : float, optional - quantile of the gradient matrix to populate the cost matrix when 0 + quantile of the gradient matrix to populate the cost matrix when 0 (default: 1) numItermax : int, optional Max number of iterations @@ -566,7 +565,7 @@ def partial_gromov_wasserstein2(C1, C2, p, q, m=None, nb_dummies=1, G0=None, where : - M is the metric cost matrix - - :math:`\Omega` is the entropic regularization term + - :math:`\Omega` is the entropic regularization term :math:`\Omega=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` - a and b are the sample weights - m is the amount of mass to be transported @@ -591,7 +590,7 @@ def partial_gromov_wasserstein2(C1, C2, p, q, m=None, nb_dummies=1, G0=None, G0 : ndarray, shape (ns, nt), optional Initialisation of the transportation matrix thres : float, optional - quantile of the gradient matrix to populate the cost matrix when 0 + quantile of the gradient matrix to populate the cost matrix when 0 (default: 1) numItermax : int, optional Max number of iterations @@ -666,7 +665,7 @@ def entropic_partial_wasserstein(a, b, M, reg, m=None, numItermax=1000, where : - M is the metric cost matrix - - :math:`\Omega` is the entropic regularization term + - :math:`\Omega` is the entropic regularization term :math:`\Omega=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` - a and b are the sample weights - m is the amount of mass to be transported @@ -754,7 +753,7 @@ def entropic_partial_wasserstein(a, b, M, reg, m=None, numItermax=1000, K = np.empty(M.shape, dtype=M.dtype) np.divide(M, -reg, out=K) np.exp(K, out=K) - np.multiply(K, m/np.sum(K), out=K) + np.multiply(K, m / np.sum(K), out=K) err, cpt = 1, 0 @@ -809,7 +808,7 @@ def entropic_partial_gromov_wasserstein(C1, C2, p, q, reg, m=None, G0=None, - C2 is the metric cost matrix in the target space - p and q are the sample weights - L : quadratic loss function - - :math:`\Omega` is the entropic regularization term + - :math:`\Omega` is the entropic regularization term :math:`\Omega=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` - m is the amount of mass to be transported @@ -944,7 +943,7 @@ def entropic_partial_gromov_wasserstein2(C1, C2, p, q, reg, m=None, G0=None, - C2 is the metric cost matrix in the target space - p and q are the sample weights - L : quadratic loss function - - :math:`\Omega` is the entropic regularization term + - :math:`\Omega` is the entropic regularization term :math:`\Omega=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` - m is the amount of mass to be transported diff --git a/ot/unbalanced.py b/ot/unbalanced.py index 66a8830..23f6607 100644 --- a/ot/unbalanced.py +++ b/ot/unbalanced.py @@ -14,7 +14,7 @@ from scipy.special import logsumexp # from .utils import unif, dist -def sinkhorn_unbalanced(a, b, M, reg, reg_m, method='sinkhorn', div = "TV", numItermax=1000, +def sinkhorn_unbalanced(a, b, M, reg, reg_m, method='sinkhorn', numItermax=1000, stopThr=1e-6, verbose=False, log=False, **kwargs): r""" Solve the unbalanced entropic regularization optimal transport problem @@ -120,20 +120,20 @@ def sinkhorn_unbalanced(a, b, M, reg, reg_m, method='sinkhorn', div = "TV", numI """ if method.lower() == 'sinkhorn': - return sinkhorn_knopp_unbalanced(a, b, M, reg, reg_m, div, + return sinkhorn_knopp_unbalanced(a, b, M, reg, reg_m, numItermax=numItermax, stopThr=stopThr, verbose=verbose, log=log, **kwargs) elif method.lower() == 'sinkhorn_stabilized': - return sinkhorn_stabilized_unbalanced(a, b, M, reg, reg_m, div, + return sinkhorn_stabilized_unbalanced(a, b, M, reg, reg_m, numItermax=numItermax, stopThr=stopThr, verbose=verbose, log=log, **kwargs) elif method.lower() in ['sinkhorn_reg_scaling']: warnings.warn('Method not implemented yet. Using classic Sinkhorn Knopp') - return sinkhorn_knopp_unbalanced(a, b, M, reg, reg_m, reg, + return sinkhorn_knopp_unbalanced(a, b, M, reg, reg_m, numItermax=numItermax, stopThr=stopThr, verbose=verbose, log=log, **kwargs) @@ -261,8 +261,8 @@ def sinkhorn_unbalanced2(a, b, M, reg, reg_m, method='sinkhorn', else: raise ValueError('Unknown method %s.' % method) -# TODO: update the doc -def sinkhorn_knopp_unbalanced(a, b, M, reg, reg_m, div="KL", numItermax=1000, + +def sinkhorn_knopp_unbalanced(a, b, M, reg, reg_m, numItermax=1000, stopThr=1e-6, verbose=False, log=False, **kwargs): r""" Solve the entropic regularization unbalanced optimal transport problem and return the loss @@ -349,7 +349,6 @@ def sinkhorn_knopp_unbalanced(a, b, M, reg, reg_m, div="KL", numItermax=1000, """ a = np.asarray(a, dtype=np.float64) - print(a) b = np.asarray(b, dtype=np.float64) M = np.asarray(M, dtype=np.float64) @@ -377,39 +376,24 @@ def sinkhorn_knopp_unbalanced(a, b, M, reg, reg_m, div="KL", numItermax=1000, else: u = np.ones(dim_a) / dim_a v = np.ones(dim_b) / dim_b - u = np.ones(dim_a) - v = np.ones(dim_b) # Next 3 lines equivalent to K= np.exp(-M/reg), but faster to compute K = np.empty(M.shape, dtype=M.dtype) - np.true_divide(M, -reg, out=K) + np.divide(M, -reg, out=K) np.exp(K, out=K) - - if div == "KL": - fi = reg_m / (reg_m + reg) - elif div == "TV": - fi = reg_m / reg + + fi = reg_m / (reg_m + reg) err = 1. - - dx = np.ones(dim_a) / dim_a - dy = np.ones(dim_b) / dim_b - z = 1 for i in range(numItermax): uprev = u vprev = v - Kv = z*K.dot(v*dy) - u = scaling_iter_prox(Kv, a, fi, div) - #u = (a / Kv) ** fi - Ktu = z*K.T.dot(u*dx) - v = scaling_iter_prox(Ktu, b, fi, div) - #v = (b / Ktu) ** fi - #print(v*dy) - z = np.dot((u*dx).T, np.dot(K,v*dy))/0.35 - print(z) - + Kv = K.dot(v) + u = (a / Kv) ** fi + Ktu = K.T.dot(u) + v = (b / Ktu) ** fi if (np.any(Ktu == 0.) or np.any(np.isnan(u)) or np.any(np.isnan(v)) @@ -450,12 +434,12 @@ def sinkhorn_knopp_unbalanced(a, b, M, reg, reg_m, div="KL", numItermax=1000, if log: return u[:, None] * K * v[None, :], log else: - return z*u[:, None] * K * v[None, :] + return u[:, None] * K * v[None, :] + -# TODO: update the doc -def sinkhorn_stabilized_unbalanced(a, b, M, reg, reg_m, div = "KL", tau=1e5, - numItermax=1000, stopThr=1e-6, - verbose=False, log=False, **kwargs): +def sinkhorn_stabilized_unbalanced(a, b, M, reg, reg_m, tau=1e5, numItermax=1000, + stopThr=1e-6, verbose=False, log=False, + **kwargs): r""" Solve the entropic regularization unbalanced optimal transport problem and return the loss @@ -580,10 +564,7 @@ def sinkhorn_stabilized_unbalanced(a, b, M, reg, reg_m, div = "KL", tau=1e5, np.divide(M, -reg, out=K) np.exp(K, out=K) - if div == "KL": - fi = reg_m / (reg_m + reg) - elif div == "TV": - fi = reg_m / reg + fi = reg_m / (reg_m + reg) cpt = 0 err = 1. @@ -669,15 +650,6 @@ def sinkhorn_stabilized_unbalanced(a, b, M, reg, reg_m, div = "KL", tau=1e5, else: return ot_matrix -def scaling_iter_prox(s, p, fi, div): - if div == "KL": - return (p / s) ** fi - elif div == "TV": - return np.minimum(s*np.exp(fi), np.maximum(s*np.exp(-fi), p)) / s - else: - raise ValueError("Unknown divergence '%s'." % div) - - def barycenter_unbalanced_stabilized(A, M, reg, reg_m, weights=None, tau=1e3, numItermax=1000, stopThr=1e-6, -- cgit v1.2.3 From 18b64556aaa477b5499dc05110c96d32b04147ff Mon Sep 17 00:00:00 2001 From: Laetitia Chapel Date: Thu, 16 Apr 2020 13:56:53 +0200 Subject: partial with doctests --- ot/partial.py | 38 ++++++++++++++++++-------------------- 1 file changed, 18 insertions(+), 20 deletions(-) (limited to 'ot/partial.py') diff --git a/ot/partial.py b/ot/partial.py index 3425acb..8698d9d 100755 --- a/ot/partial.py +++ b/ot/partial.py @@ -80,10 +80,10 @@ def partial_wasserstein_lagrange(a, b, M, reg_m=None, nb_dummies=1, log=False, >>> M = [[0., 1.], [2., 3.]] >>> np.round(partial_wasserstein_lagrange(a,b,M), 2) array([[0.1, 0. ], - [0. , 0.1]]) + [0. , 0.1]]) >>> np.round(partial_wasserstein_lagrange(a,b,M,reg_m=2), 2) array([[0.1, 0. ], - [0. , 0. ]]) + [0. , 0. ]]) References ---------- @@ -199,10 +199,10 @@ def partial_wasserstein(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): >>> M = [[0., 1.], [2., 3.]] >>> np.round(partial_wasserstein(a,b,M), 2) array([[0.1, 0. ], - [0. , 0.1]]) + [0. , 0.1]]) >>> np.round(partial_wasserstein(a,b,M,m=0.1), 2) array([[0.1, 0. ], - [0. , 0. ]]) + [0. , 0. ]]) References ---------- @@ -466,14 +466,14 @@ def partial_gromov_wasserstein(C1, C2, p, q, m=None, nb_dummies=1, G0=None, >>> C2 = sp.spatial.distance.cdist(y, y) >>> np.round(partial_gromov_wasserstein(C1, C2, a, b),2) array([[0. , 0.25, 0. , 0. ], - [0.25, 0. , 0. , 0. ], - [0. , 0. , 0.25, 0. ], - [0. , 0. , 0. , 0.25]]) + [0.25, 0. , 0. , 0. ], + [0. , 0. , 0.25, 0. ], + [0. , 0. , 0. , 0.25]]) >>> np.round(partial_gromov_wasserstein(C1, C2, a, b, m=0.25),2) array([[0. , 0. , 0. , 0. ], - [0. , 0. , 0. , 0. ], - [0. , 0. , 0. , 0. ], - [0. , 0. , 0. , 0.25]]) + [0. , 0. , 0. , 0. ], + [0. , 0. , 0. , 0. ], + [0. , 0. , 0. , 0.25]]) References ---------- @@ -711,8 +711,7 @@ def entropic_partial_wasserstein(a, b, M, reg, m=None, numItermax=1000, >>> M = [[0., 1.], [2., 3.]] >>> np.round(entropic_partial_wasserstein(a, b, M, 1, 0.1), 2) array([[0.06, 0.02], - [0.01, 0. ]]) - + [0.01, 0. ]]) References ---------- @@ -849,15 +848,14 @@ def entropic_partial_gromov_wasserstein(C1, C2, p, q, reg, m=None, G0=None, >>> C2 = sp.spatial.distance.cdist(y, y) >>> np.round(entropic_partial_gromov_wasserstein(C1, C2, a, b,50), 2) array([[0.12, 0.13, 0. , 0. ], - [0.13, 0.12, 0. , 0. ], - [0. , 0. , 0.25, 0. ], - [0. , 0. , 0. , 0.25]]) - >>> np.round(entropic_partial_gromov_wasserstein(C1, C2, a, b, 50, m=0.25) - , 2) + [0.13, 0.12, 0. , 0. ], + [0. , 0. , 0.25, 0. ], + [0. , 0. , 0. , 0.25]]) + >>> np.round(entropic_partial_gromov_wasserstein(C1, C2, a, b, 50, m=0.25), 2) array([[0.02, 0.03, 0. , 0.03], - [0.03, 0.03, 0. , 0.03], - [0. , 0. , 0.03, 0. ], - [0.02, 0.02, 0. , 0.03]]) + [0.03, 0.03, 0. , 0.03], + [0. , 0. , 0.03, 0. ], + [0.02, 0.02, 0. , 0.03]]) Returns ------- -- cgit v1.2.3 From ef7c11a5df3cf6c82864472f0cfa65d6b2036f2f Mon Sep 17 00:00:00 2001 From: Laetitia Chapel Date: Thu, 16 Apr 2020 15:52:00 +0200 Subject: partial with python 3.8 --- .travis.yml | 2 +- ot/partial.py | 12 ++++++------ test/test_partial.py | 9 ++++----- 3 files changed, 11 insertions(+), 12 deletions(-) (limited to 'ot/partial.py') diff --git a/.travis.yml b/.travis.yml index 072bc55..7ff1b3c 100644 --- a/.travis.yml +++ b/.travis.yml @@ -16,7 +16,7 @@ matrix: python: 3.7 - os: linux sudo: required - python: 2.7 + python: 3.8 # - os: osx # sudo: required # language: generic diff --git a/ot/partial.py b/ot/partial.py index 8698d9d..726a590 100755 --- a/ot/partial.py +++ b/ot/partial.py @@ -232,7 +232,7 @@ def partial_wasserstein(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): b_extended = np.append(b, [(np.sum(a) - m) / nb_dummies] * nb_dummies) a_extended = np.append(a, [(np.sum(b) - m) / nb_dummies] * nb_dummies) - M_extended = np.ones((len(a_extended), len(b_extended))) * 0 + M_extended = np.zeros((len(a_extended), len(b_extended))) M_extended[-1, -1] = np.max(M) * 1e5 M_extended[:len(a), :len(b)] = M @@ -510,9 +510,9 @@ def partial_gromov_wasserstein(C1, C2, p, q, m=None, nb_dummies=1, G0=None, Gprev = G0 M = gwgrad_partial(C1, C2, G0) - M[M < eps] = np.quantile(M[M > eps], thres) + M[M < eps] = np.quantile(M, thres) - M_emd = np.ones(dim_G_extended) * np.max(M) * 1e2 + M_emd = np.zeros(dim_G_extended) M_emd[:len(p), :len(q)] = M M_emd[-nb_dummies:, -nb_dummies:] = np.max(M) * 1e5 M_emd = np.asarray(M_emd, dtype=np.float64) @@ -729,8 +729,8 @@ def entropic_partial_wasserstein(a, b, M, reg, m=None, numItermax=1000, M = np.asarray(M, dtype=np.float64) dim_a, dim_b = M.shape - dx = np.ones(dim_a) - dy = np.ones(dim_b) + dx = np.ones(dim_a, dtype=np.float64) + dy = np.ones(dim_b, dtype=np.float64) if len(a) == 0: a = np.ones(dim_a, dtype=np.float64) / dim_a @@ -738,7 +738,7 @@ def entropic_partial_wasserstein(a, b, M, reg, m=None, numItermax=1000, b = np.ones(dim_b, dtype=np.float64) / dim_b if m is None: - m = np.min((np.sum(a), np.sum(b))) + m = np.min((np.sum(a), np.sum(b))) * 1.0 if m < 0: raise ValueError("Problem infeasible. Parameter m should be greater" " than 0.") diff --git a/test/test_partial.py b/test/test_partial.py index fbcd3c2..1799fd4 100755 --- a/test/test_partial.py +++ b/test/test_partial.py @@ -93,10 +93,7 @@ def test_partial_gromov_wasserstein(): m = 2 / 3 res0, log0 = ot.partial.partial_gromov_wasserstein(C1, C3, p, q, m=m, log=True) - res, log = ot.partial.entropic_partial_gromov_wasserstein(C1, C3, p, q, 10, - m=m, log=True) np.testing.assert_allclose(res0, 0, atol=1e-1, rtol=1e-1) - np.testing.assert_allclose(res, 0, atol=1e-1, rtol=1e-1) C1 = sp.spatial.distance.cdist(xs, xs) C2 = sp.spatial.distance.cdist(xt, xt) @@ -123,8 +120,10 @@ def test_partial_gromov_wasserstein(): m = 2 / 3 res0, log0 = ot.partial.partial_gromov_wasserstein(C1, C2, p, q, m=m, log=True) - res, log = ot.partial.entropic_partial_gromov_wasserstein(C1, C2, p, q, 10, - m=m, log=True) + res, log = ot.partial.entropic_partial_gromov_wasserstein(C1, C2, p, q, + 100, m=m, + log=True) + # check constratints np.testing.assert_equal( res0.sum(1) <= p, [True] * len(p)) # cf convergence wasserstein -- cgit v1.2.3 From dc942ac386423870277ea69fae723f216ea61030 Mon Sep 17 00:00:00 2001 From: Laetitia Chapel Date: Fri, 17 Apr 2020 09:08:02 +0200 Subject: partial with readme updated --- README.md | 4 +++- ot/partial.py | 12 ++++++------ 2 files changed, 9 insertions(+), 7 deletions(-) (limited to 'ot/partial.py') diff --git a/README.md b/README.md index b6baf14..fecaa35 100644 --- a/README.md +++ b/README.md @@ -259,4 +259,6 @@ You can also post bug reports and feature requests in Github issues. Make sure t [26] Alaya M. Z., Bérar M., Gasso G., Rakotomamonjy A. (2019). [Screening Sinkhorn Algorithm for Regularized Optimal Transport](https://papers.nips.cc/paper/9386-screening-sinkhorn-algorithm-for-regularized-optimal-transport), Advances in Neural Information Processing Systems 33 (NeurIPS). -[27] Redko I., Courty N., Flamary R., Tuia D. (2019). [Optimal Transport for Multi-source Domain Adaptation under Target Shift](http://proceedings.mlr.press/v89/redko19a.html), Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics (AISTATS) 22, 2019. \ No newline at end of file +[27] Redko I., Courty N., Flamary R., Tuia D. (2019). [Optimal Transport for Multi-source Domain Adaptation under Target Shift](http://proceedings.mlr.press/v89/redko19a.html), Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics (AISTATS) 22, 2019. + +[28] Chapel, L., Alaya, M., Gasso, G. (2019). [Partial Gromov-Wasserstein with Applications on Positive-Unlabeled Learning"] (https://arxiv.org/abs/2002.08276), arXiv preprint arXiv:2002.08276. \ No newline at end of file diff --git a/ot/partial.py b/ot/partial.py index 726a590..d32e054 100755 --- a/ot/partial.py +++ b/ot/partial.py @@ -209,7 +209,7 @@ def partial_wasserstein(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): .. [26] Caffarelli, L. A., & McCann, R. J. (2010) Free boundaries in optimal transport and Monge-Ampere obstacle problems. Annals of mathematics, 673-730. - .. [27] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + .. [28] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- Wasserstein with Applications on Positive-Unlabeled Learning". arXiv preprint arXiv:2002.08276. @@ -314,7 +314,7 @@ def partial_wasserstein2(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): .. [26] Caffarelli, L. A., & McCann, R. J. (2010) Free boundaries in optimal transport and Monge-Ampere obstacle problems. Annals of mathematics, 673-730. - .. [27] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + .. [28] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- Wasserstein with Applications on Positive-Unlabeled Learning". arXiv preprint arXiv:2002.08276. """ @@ -411,7 +411,7 @@ def partial_gromov_wasserstein(C1, C2, p, q, m=None, nb_dummies=1, G0=None, - a and b are the sample weights - m is the amount of mass to be transported - The formulation of the problem has been proposed in [27]_ + The formulation of the problem has been proposed in [28]_ Parameters @@ -477,7 +477,7 @@ def partial_gromov_wasserstein(C1, C2, p, q, m=None, nb_dummies=1, G0=None, References ---------- - .. [27] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + .. [28] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- Wasserstein with Applications on Positive-Unlabeled Learning". arXiv preprint arXiv:2002.08276. @@ -570,7 +570,7 @@ def partial_gromov_wasserstein2(C1, C2, p, q, m=None, nb_dummies=1, G0=None, - a and b are the sample weights - m is the amount of mass to be transported - The formulation of the problem has been proposed in [27]_ + The formulation of the problem has been proposed in [28]_ Parameters @@ -627,7 +627,7 @@ def partial_gromov_wasserstein2(C1, C2, p, q, m=None, nb_dummies=1, G0=None, References ---------- - .. [27] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + .. [28] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- Wasserstein with Applications on Positive-Unlabeled Learning". arXiv preprint arXiv:2002.08276. -- cgit v1.2.3 From 429abe06d53e1ebdd2492b275f70ba1bfe751f0f Mon Sep 17 00:00:00 2001 From: Laetitia Chapel Date: Fri, 17 Apr 2020 11:49:28 +0200 Subject: partial added on quick start guide --- README.md | 4 +- docs/source/quickstart.rst | 64 ++++++++++++++++++++++++++++ docs/source/readme.rst | 104 +++++++++++++++++++++++++++------------------ ot/partial.py | 68 ++++++++++++++++++++++------- 4 files changed, 181 insertions(+), 59 deletions(-) (limited to 'ot/partial.py') diff --git a/README.md b/README.md index fecaa35..1f62c2c 100644 --- a/README.md +++ b/README.md @@ -261,4 +261,6 @@ You can also post bug reports and feature requests in Github issues. Make sure t [27] Redko I., Courty N., Flamary R., Tuia D. (2019). [Optimal Transport for Multi-source Domain Adaptation under Target Shift](http://proceedings.mlr.press/v89/redko19a.html), Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics (AISTATS) 22, 2019. -[28] Chapel, L., Alaya, M., Gasso, G. (2019). [Partial Gromov-Wasserstein with Applications on Positive-Unlabeled Learning"] (https://arxiv.org/abs/2002.08276), arXiv preprint arXiv:2002.08276. \ No newline at end of file +[28] Caffarelli, L. A., McCann, R. J. (2020). [Free boundaries in optimal transport and Monge-Ampere obstacle problems](http://www.math.toronto.edu/~mccann/papers/annals2010.pdf), Annals of mathematics, 673-730. + +[29] Chapel, L., Alaya, M., Gasso, G. (2019). [Partial Gromov-Wasserstein with Applications on Positive-Unlabeled Learning](https://arxiv.org/abs/2002.08276), arXiv preprint arXiv:2002.08276. \ No newline at end of file diff --git a/docs/source/quickstart.rst b/docs/source/quickstart.rst index 978eaff..d56f812 100644 --- a/docs/source/quickstart.rst +++ b/docs/source/quickstart.rst @@ -645,6 +645,53 @@ implemented the main function :any:`ot.barycenter_unbalanced`. - :any:`auto_examples/plot_UOT_barycenter_1D` +Partial optimal transport +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Partial OT is a variant of the optimal transport problem when only a fixed amount of mass m +is to be transported. The partial OT metric between two histograms a and b is defined as [28]_: + +.. math:: + \gamma = \arg\min_\gamma <\gamma,M>_F + + s.t. + \gamma\geq 0 \\ + \gamma 1 \leq a\\ + \gamma^T 1 \leq b\\ + 1^T \gamma^T 1 = m \leq \min\{\|a\|_1, \|b\|_1\} + + +Interestingly the problem can be casted into a regular OT problem by adding reservoir points +in which the surplus mass is sent [29]_. We provide a solver for partial OT +in :any:`ot.partial`. The exact resolution of the problem is computed in :any:`ot.partial.partial_wasserstein` +and :any:`ot.partial.partial_wasserstein2` that return respectively the OT matrix and the value of the +linear term. The entropic solution of the problem is computed in :any:`ot.partial.entropic_partial_wasserstein` +(see [3]_). + +The partial Gromov-Wasserstein formulation of the problem + +.. math:: + GW = \min_\gamma \sum_{i,j,k,l} L(C1_{i,k},C2_{j,l})*\gamma_{i,j}*\gamma_{k,l} + + s.t. + \gamma\geq 0 \\ + \gamma 1 \leq a\\ + \gamma^T 1 \leq b\\ + 1^T \gamma^T 1 = m \leq \min\{\|a\|_1, \|b\|_1\} + +is computed in :any:`ot.partial.partial_gromov_wasserstein` and in +:any:`ot.partial.entropic_partial_gromov_wasserstein` when considering the entropic +regularization of the problem. + + +.. hint:: + + Examples of the use of :any:`ot.partial` are available in : + + - :any:`auto_examples/plot_partial` + + + Gromov-Wasserstein ^^^^^^^^^^^^^^^^^^ @@ -921,3 +968,20 @@ References .. [25] Frogner C., Zhang C., Mobahi H., Araya-Polo M., Poggio T. : Learning with a Wasserstein Loss, Advances in Neural Information Processing Systems (NIPS) 2015 + +.. [26] Alaya M. Z., Bérar M., Gasso G., Rakotomamonjy A. (2019). Screening Sinkhorn + Algorithm for Regularized Optimal Transport , + Advances in Neural Information Processing Systems 33 (NeurIPS). + +.. [27] Redko I., Courty N., Flamary R., Tuia D. (2019). Optimal Transport for Multi-source + Domain Adaptation under Target Shift , + Proceedings of the Twenty-Second International Conference on Artificial Intelligence + and Statistics (AISTATS) 22, 2019. + +.. [28] Caffarelli, L. A., McCann, R. J. (2020). Free boundaries in optimal transport and + Monge-Ampere obstacle problems , + Annals of mathematics, 673-730. + +.. [29] Chapel, L., Alaya, M., Gasso, G. (2019). Partial Gromov-Wasserstein with + Applications on Positive-Unlabeled Learning , + arXiv preprint arXiv:2002.08276. diff --git a/docs/source/readme.rst b/docs/source/readme.rst index d5f2161..6d98dc5 100644 --- a/docs/source/readme.rst +++ b/docs/source/readme.rst @@ -36,6 +36,9 @@ It provides the following solvers: problem [18] and dual problem [19]) - Non regularized free support Wasserstein barycenters [20]. - Unbalanced OT with KL relaxation distance and barycenter [10, 25]. +- Screening Sinkhorn Algorithm for OT [26]. +- JCPOT algorithm for multi-source domain adaptation with target shift + [27]. Some demonstrations (both in Python and Jupyter Notebook format) are available in the examples folder. @@ -48,19 +51,19 @@ 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} - } + @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: +C++ compiler for building/installing the EMD solver and relies on the +following Python modules: - Numpy (>=1.11) - Scipy (>=1.0) @@ -75,19 +78,19 @@ be installed prior to installing POT. This can be done easily with :: - pip install numpy cython + pip install numpy cython You can install the toolbox through PyPI with: :: - pip install POT + pip install POT or get the very latest version by downloading it and then running: :: - python setup.py install --user # for user install (no root) + python setup.py install --user # for user install (no root) Anaconda installation with conda-forge ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -98,7 +101,7 @@ required dependencies: :: - conda install -c conda-forge pot + conda install -c conda-forge pot Post installation check ^^^^^^^^^^^^^^^^^^^^^^^ @@ -108,7 +111,7 @@ without errors: .. code:: python - import ot + import ot Note that for easier access the module is name ot instead of pot. @@ -121,9 +124,9 @@ below - **ot.dr** (Wasserstein dimensionality reduction) depends on autograd and pymanopt that can be installed with: - :: +:: - pip install pymanopt autograd + pip install pymanopt autograd - **ot.gpu** (GPU accelerated OT) depends on cupy that have to be installed following instructions on `this @@ -139,36 +142,36 @@ Short examples - Import the toolbox - .. code:: python +.. code:: python - import ot + import ot - Compute Wasserstein distances - .. code:: python +.. code:: python - # a,b are 1D histograms (sum to 1 and positive) - # M is the ground cost matrix - Wd=ot.emd2(a,b,M) # exact linear program - Wd_reg=ot.sinkhorn2(a,b,M,reg) # entropic regularized OT - # if b is a matrix compute all distances to a and return a vector + # a,b are 1D histograms (sum to 1 and positive) + # M is the ground cost matrix + Wd=ot.emd2(a,b,M) # exact linear program + Wd_reg=ot.sinkhorn2(a,b,M,reg) # entropic regularized OT + # if b is a matrix compute all distances to a and return a vector - Compute OT matrix - .. code:: python +.. code:: python - # a,b are 1D histograms (sum to 1 and positive) - # M is the ground cost matrix - T=ot.emd(a,b,M) # exact linear program - T_reg=ot.sinkhorn(a,b,M,reg) # entropic regularized OT + # a,b are 1D histograms (sum to 1 and positive) + # M is the ground cost matrix + T=ot.emd(a,b,M) # exact linear program + T_reg=ot.sinkhorn(a,b,M,reg) # entropic regularized OT - Compute Wasserstein barycenter - .. code:: python +.. code:: python - # A is a n*d matrix containing d 1D histograms - # M is the ground cost matrix - ba=ot.barycenter(A,M,reg) # reg is regularization parameter + # A is a n*d matrix containing d 1D histograms + # M is the ground cost matrix + ba=ot.barycenter(A,M,reg) # reg is regularization parameter Examples and Notebooks ~~~~~~~~~~~~~~~~~~~~~~ @@ -207,6 +210,10 @@ want a quick look: Wasserstein `__ - `Gromov Wasserstein Barycenter `__ +- `Fused Gromov + Wasserstein `__ +- `Fused Gromov Wasserstein + Barycenter `__ You can also see the notebooks with `Jupyter nbviewer `__. @@ -237,6 +244,7 @@ The contributors to this library are - `Vayer Titouan `__ - `Hicham Janati `__ (Unbalanced OT) - `Romain Tavenard `__ (1d Wasserstein) +- `Mokhtar Z. Alaya `__ (Screenkhorn) This toolbox benefit a lot from open source research and we would like to thank the following persons for providing some code (in various @@ -274,11 +282,11 @@ References [1] Bonneel, N., Van De Panne, M., Paris, S., & Heidrich, W. (2011, December). `Displacement interpolation using Lagrangian mass transport `__. -In ACM Transactions on Graphics (TOG) (Vol. 30, No. 6, p. 158). ACM. +In ACM Transactions on Graphics (TOG) (Vol. 30, No. 6, p. 158). ACM. [2] Cuturi, M. (2013). `Sinkhorn distances: Lightspeed computation of optimal transport `__. In Advances -in Neural Information Processing Systems (pp. 2292-2300). +in Neural Information Processing Systems (pp. 2292-2300). [3] Benamou, J. D., Carlier, G., Cuturi, M., Nenna, L., & Peyré, G. (2015). `Iterative Bregman projections for regularized transportation @@ -387,17 +395,29 @@ and Statistics, (AISTATS) 21, 2018 graphs `__ Proceedings of the 36th International Conference on Machine Learning (ICML). -[25] Frogner C., Zhang C., Mobahi H., Araya-Polo M., Poggio T. (2019). +[25] Frogner C., Zhang C., Mobahi H., Araya-Polo M., Poggio T. (2015). `Learning with a Wasserstein Loss `__ Advances in Neural Information Processing Systems (NIPS). -[26] Caffarelli, L. A., McCann, R. J. (2020). `Free boundaries in optimal transport and -Monge-Ampere obstacle problems `__, -Annals of mathematics, 673-730. - -[27] Chapel, L., Alaya, M., Gasso, G. (2019). `Partial Gromov-Wasserstein with Applications -on Positive-Unlabeled Learning `__. arXiv preprint -arXiv:2002.08276. +[26] Alaya M. Z., Bérar M., Gasso G., Rakotomamonjy A. (2019). +`Screening Sinkhorn Algorithm for Regularized Optimal +Transport `__, +Advances in Neural Information Processing Systems 33 (NeurIPS). + +[27] Redko I., Courty N., Flamary R., Tuia D. (2019). `Optimal Transport +for Multi-source Domain Adaptation under Target +Shift `__, Proceedings +of the Twenty-Second International Conference on Artificial Intelligence +and Statistics (AISTATS) 22, 2019. + +[28] Caffarelli, L. A., McCann, R. J. (2020). [Free boundaries in +optimal transport and Monge-Ampere obstacle problems] +(http://www.math.toronto.edu/~mccann/papers/annals2010.pdf), Annals of +mathematics, 673-730. + +[29] Chapel, L., Alaya, M., Gasso, G. (2019). [Partial +Gromov-Wasserstein with Applications on Positive-Unlabeled Learning"] +(https://arxiv.org/abs/2002.08276), arXiv preprint arXiv:2002.08276. .. |PyPI version| image:: https://badge.fury.io/py/POT.svg :target: https://badge.fury.io/py/POT diff --git a/ot/partial.py b/ot/partial.py index d32e054..f325d98 100755 --- a/ot/partial.py +++ b/ot/partial.py @@ -1,4 +1,3 @@ -#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Partial OT @@ -30,7 +29,9 @@ def partial_wasserstein_lagrange(a, b, M, reg_m=None, nb_dummies=1, log=False, 1^T \gamma^T 1 = m \leq \min\{\|a\|_1, \|b\|_1\} - or equivalently: + or equivalently (see Chizat, L., Peyré, G., Schmitzer, B., & Vialard, F. X. + (2018). An interpolating distance between optimal transport and Fisher–Rao + metrics. Foundations of Computational Mathematics, 18(1), 1-44.) .. math:: \gamma = \arg\min_\gamma <\gamma,M>_F + \sqrt(\lambda/2) @@ -47,7 +48,8 @@ def partial_wasserstein_lagrange(a, b, M, reg_m=None, nb_dummies=1, log=False, - :math:`\lambda` is the lagragian cost. Tuning its value allows attaining a given mass to be transported m - The formulation of the problem has been proposed in [26]_ + The formulation of the problem has been proposed in [28]_ + Parameters ---------- @@ -59,9 +61,17 @@ def partial_wasserstein_lagrange(a, b, M, reg_m=None, nb_dummies=1, log=False, cost matrix for the quadratic cost reg_m : float, optional Lagragian cost + nb_dummies : int, optional, default:1 + number of reservoir points to be added (to avoid numerical + instabilities, increase its value if an error is raised) log : bool, optional record log if True + .. warning:: + When dealing with a large number of points, the EMD solver may face + some instabilities, especially when the mass associated to the dummy + point is large. To avoid them, increase the number of dummy points + (allows a smoother repartition of the mass over the points). Returns ------- @@ -88,7 +98,7 @@ def partial_wasserstein_lagrange(a, b, M, reg_m=None, nb_dummies=1, log=False, References ---------- - .. [26] Caffarelli, L. A., & McCann, R. J. (2010) Free boundaries in + .. [28] Caffarelli, L. A., & McCann, R. J. (2010) Free boundaries in optimal transport and Monge-Ampere obstacle problems. Annals of mathematics, 673-730. @@ -182,6 +192,13 @@ def partial_wasserstein(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): record log if True + .. warning:: + When dealing with a large number of points, the EMD solver may face + some instabilities, especially when the mass associated to the dummy + point is large. To avoid them, increase the number of dummy points + (allows a smoother repartition of the mass over the points). + + Returns ------- :math:`gamma` : (dim_a x dim_b) ndarray @@ -206,10 +223,10 @@ def partial_wasserstein(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): References ---------- - .. [26] Caffarelli, L. A., & McCann, R. J. (2010) Free boundaries in + .. [28] Caffarelli, L. A., & McCann, R. J. (2010) Free boundaries in optimal transport and Monge-Ampere obstacle problems. Annals of mathematics, 673-730. - .. [28] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + .. [29] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- Wasserstein with Applications on Positive-Unlabeled Learning". arXiv preprint arXiv:2002.08276. @@ -289,6 +306,13 @@ def partial_wasserstein2(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): record log if True + .. warning:: + When dealing with a large number of points, the EMD solver may face + some instabilities, especially when the mass associated to the dummy + point is large. To avoid them, increase the number of dummy points + (allows a smoother repartition of the mass over the points). + + Returns ------- :math:`gamma` : (dim_a x dim_b) ndarray @@ -311,10 +335,10 @@ def partial_wasserstein2(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): References ---------- - .. [26] Caffarelli, L. A., & McCann, R. J. (2010) Free boundaries in + .. [28] Caffarelli, L. A., & McCann, R. J. (2010) Free boundaries in optimal transport and Monge-Ampere obstacle problems. Annals of mathematics, 673-730. - .. [28] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + .. [29] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- Wasserstein with Applications on Positive-Unlabeled Learning". arXiv preprint arXiv:2002.08276. """ @@ -411,7 +435,7 @@ def partial_gromov_wasserstein(C1, C2, p, q, m=None, nb_dummies=1, G0=None, - a and b are the sample weights - m is the amount of mass to be transported - The formulation of the problem has been proposed in [28]_ + The formulation of the problem has been proposed in [29]_ Parameters @@ -435,13 +459,12 @@ def partial_gromov_wasserstein(C1, C2, p, q, m=None, nb_dummies=1, G0=None, (default: 1) numItermax : int, optional Max number of iterations + tol : float, optional + tolerance for stopping iterations log : bool, optional return log if True verbose : bool, optional Print information along iterations - armijo : bool, optional - If True the steps of the line-search is found via an armijo research. Else closed form is used. - If there is convergence issues use False. **kwargs : dict parameters can be directly passed to the emd solver @@ -477,7 +500,7 @@ def partial_gromov_wasserstein(C1, C2, p, q, m=None, nb_dummies=1, G0=None, References ---------- - .. [28] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + .. [29] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- Wasserstein with Applications on Positive-Unlabeled Learning". arXiv preprint arXiv:2002.08276. @@ -546,7 +569,7 @@ def partial_gromov_wasserstein(C1, C2, p, q, m=None, nb_dummies=1, G0=None, def partial_gromov_wasserstein2(C1, C2, p, q, m=None, nb_dummies=1, G0=None, - thres=0.75, numItermax=1000, tol=1e-7, + thres=1, numItermax=1000, tol=1e-7, log=False, verbose=False, **kwargs): r""" Solves the partial optimal transport problem @@ -570,7 +593,7 @@ def partial_gromov_wasserstein2(C1, C2, p, q, m=None, nb_dummies=1, G0=None, - a and b are the sample weights - m is the amount of mass to be transported - The formulation of the problem has been proposed in [28]_ + The formulation of the problem has been proposed in [29]_ Parameters @@ -594,6 +617,8 @@ def partial_gromov_wasserstein2(C1, C2, p, q, m=None, nb_dummies=1, G0=None, (default: 1) numItermax : int, optional Max number of iterations + tol : float, optional + tolerance for stopping iterations log : bool, optional return log if True verbose : bool, optional @@ -602,6 +627,13 @@ def partial_gromov_wasserstein2(C1, C2, p, q, m=None, nb_dummies=1, G0=None, parameters can be directly passed to the emd solver + .. warning:: + When dealing with a large number of points, the EMD solver may face + some instabilities, especially when the mass associated to the dummy + point is large. To avoid them, increase the number of dummy points + (allows a smoother repartition of the mass over the points). + + Returns ------- partial_gw_dist : (dim_a x dim_b) ndarray @@ -627,7 +659,7 @@ def partial_gromov_wasserstein2(C1, C2, p, q, m=None, nb_dummies=1, G0=None, References ---------- - .. [28] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + .. [29] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- Wasserstein with Applications on Positive-Unlabeled Learning". arXiv preprint arXiv:2002.08276. @@ -831,6 +863,8 @@ def entropic_partial_gromov_wasserstein(C1, C2, p, q, reg, m=None, G0=None, Initialisation of the transportation matrix numItermax : int, optional Max number of iterations + tol : float, optional + Stop threshold on error (>0) log : bool, optional return log if True verbose : bool, optional @@ -966,6 +1000,8 @@ def entropic_partial_gromov_wasserstein2(C1, C2, p, q, reg, m=None, G0=None, Initialisation of the transportation matrix numItermax : int, optional Max number of iterations + tol : float, optional + Stop threshold on error (>0) log : bool, optional return log if True verbose : bool, optional -- cgit v1.2.3 From 5ed4a27c8054397aae51ca49ddfcc8fa01e64db7 Mon Sep 17 00:00:00 2001 From: Laetitia Chapel Date: Mon, 20 Apr 2020 09:09:57 +0200 Subject: partial update refs --- ot/partial.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) (limited to 'ot/partial.py') diff --git a/ot/partial.py b/ot/partial.py index f325d98..5f4b836 100755 --- a/ot/partial.py +++ b/ot/partial.py @@ -702,7 +702,7 @@ def entropic_partial_wasserstein(a, b, M, reg, m=None, numItermax=1000, - a and b are the sample weights - m is the amount of mass to be transported - The formulation of the problem has been proposed in [3]_ + The formulation of the problem has been proposed in [3]_ (prop. 5) Parameters @@ -843,7 +843,8 @@ def entropic_partial_gromov_wasserstein(C1, C2, p, q, reg, m=None, G0=None, :math:`\Omega=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` - m is the amount of mass to be transported - The formulation of the problem has been proposed in [12]. + The formulation of the GW problem has been proposed in [12]_ and the + partial GW in [29]_. Parameters ---------- @@ -903,6 +904,9 @@ def entropic_partial_gromov_wasserstein(C1, C2, p, q, reg, m=None, G0=None, .. [12] Peyré, Gabriel, Marco Cuturi, and Justin Solomon, "Gromov-Wasserstein averaging of kernel and distance matrices." International Conference on Machine Learning (ICML). 2016. + .. [29] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + Wasserstein with Applications on Positive-Unlabeled Learning". + arXiv preprint arXiv:2002.08276. See Also -------- @@ -979,7 +983,8 @@ def entropic_partial_gromov_wasserstein2(C1, C2, p, q, reg, m=None, G0=None, :math:`\Omega=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` - m is the amount of mass to be transported - The formulation of the problem has been proposed in [12]. + The formulation of the GW problem has been proposed in [12]_ and the + partial GW in [29]_. Parameters @@ -1033,6 +1038,9 @@ def entropic_partial_gromov_wasserstein2(C1, C2, p, q, reg, m=None, G0=None, .. [12] Peyré, Gabriel, Marco Cuturi, and Justin Solomon, "Gromov-Wasserstein averaging of kernel and distance matrices." International Conference on Machine Learning (ICML). 2016. + .. [29] Chapel, L., Alaya, M., Gasso, G. (2019). "Partial Gromov- + Wasserstein with Applications on Positive-Unlabeled Learning". + arXiv preprint arXiv:2002.08276. """ partial_gw, log_gw = entropic_partial_gromov_wasserstein(C1, C2, p, q, reg, -- cgit v1.2.3 From fcf23770c5ed5252dee233b578239118d9d8ba62 Mon Sep 17 00:00:00 2001 From: Laetitia Chapel Date: Mon, 20 Apr 2020 13:50:14 +0200 Subject: doc kwargs + fix rst issues --- examples/plot_partial_wass_and_gromov.py | 24 +++++++++++------------- ot/partial.py | 6 ++++++ 2 files changed, 17 insertions(+), 13 deletions(-) (limited to 'ot/partial.py') diff --git a/examples/plot_partial_wass_and_gromov.py b/examples/plot_partial_wass_and_gromov.py index 30b3fc0..a5af441 100755 --- a/examples/plot_partial_wass_and_gromov.py +++ b/examples/plot_partial_wass_and_gromov.py @@ -1,8 +1,8 @@ # -*- coding: utf-8 -*- """ -========================== +================================================== Partial Wasserstein and Gromov-Wasserstein example -========================== +================================================== This example is designed to show how to use the Partial (Gromov-)Wassertsein distance computation in POT. @@ -50,8 +50,7 @@ pl.show() ############################################################################# # -# Compute partial Wasserstein plans and distance, -# by transporting 50% of the mass +# Compute partial Wasserstein plans and distance # ---------------------------------------------- p = ot.unif(n_samples + n_noise) @@ -113,34 +112,33 @@ pl.show() ############################################################################# # -# Compute partial Gromov-Wasserstein plans and distance, -# by transporting 100% and 2/3 of the mass +# Compute partial Gromov-Wasserstein plans and distance # ----------------------------------------------------- C1 = sp.spatial.distance.cdist(xs, xs) C2 = sp.spatial.distance.cdist(xt, xt) +# transport 100% of the mass print('-----m = 1') m = 1 -res0, log0 = ot.partial.partial_gromov_wasserstein(C1, C2, p, q, m=m, - log=True) +res0, log0 = ot.partial.partial_gromov_wasserstein(C1, C2, p, q, m=m, log=True) res, log = ot.partial.entropic_partial_gromov_wasserstein(C1, C2, p, q, 10, m=m, log=True) -print('Partial Wasserstein distance (m = 1): ' + str(log0['partial_gw_dist'])) -print('Entropic partial Wasserstein distance (m = 1): ' + - str(log['partial_gw_dist'])) +print('Wasserstein distance (m = 1): ' + str(log0['partial_gw_dist'])) +print('Entropic Wasserstein distance (m = 1): ' + str(log['partial_gw_dist'])) pl.figure(1, (10, 5)) pl.title("mass to be transported m = 1") pl.subplot(1, 2, 1) pl.imshow(res0, cmap='jet') -pl.title('Partial Wasserstein') +pl.title('Wasserstein') pl.subplot(1, 2, 2) pl.imshow(res, cmap='jet') -pl.title('Entropic partial Wasserstein') +pl.title('Entropic Wasserstein') pl.show() +# transport 2/3 of the mass print('-----m = 2/3') m = 2 / 3 res0, log0 = ot.partial.partial_gromov_wasserstein(C1, C2, p, q, m=m, log=True) diff --git a/ot/partial.py b/ot/partial.py index 5f4b836..c03ec25 100755 --- a/ot/partial.py +++ b/ot/partial.py @@ -66,6 +66,8 @@ def partial_wasserstein_lagrange(a, b, M, reg_m=None, nb_dummies=1, log=False, instabilities, increase its value if an error is raised) log : bool, optional record log if True + **kwargs : dict + parameters can be directly passed to the emd solver .. warning:: When dealing with a large number of points, the EMD solver may face @@ -190,6 +192,8 @@ def partial_wasserstein(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): instabilities, increase its value if an error is raised) log : bool, optional record log if True + **kwargs : dict + parameters can be directly passed to the emd solver .. warning:: @@ -304,6 +308,8 @@ def partial_wasserstein2(a, b, M, m=None, nb_dummies=1, log=False, **kwargs): instabilities, increase its value if an error is raised) log : bool, optional record log if True + **kwargs : dict + parameters can be directly passed to the emd solver .. warning:: -- cgit v1.2.3 From ad7aa892b47f039366a30103c1cede804811fb46 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Fri, 24 Apr 2020 16:39:29 +0200 Subject: better doc per module --- ot/__init__.py | 30 ------------------------------ ot/bregman.py | 2 +- ot/datasets.py | 2 +- ot/dr.py | 4 ++-- ot/gpu/__init__.py | 5 +++-- ot/gromov.py | 2 +- ot/optim.py | 2 +- ot/partial.py | 2 +- ot/smooth.py | 4 +++- ot/unbalanced.py | 2 +- 10 files changed, 14 insertions(+), 41 deletions(-) (limited to 'ot/partial.py') diff --git a/ot/__init__.py b/ot/__init__.py index 1e57b78..2d23610 100644 --- a/ot/__init__.py +++ b/ot/__init__.py @@ -1,35 +1,5 @@ """ -This is the main module of the POT toolbox. It provides easy access to -a number of sub-modules and functions described below. - -.. note:: - - - Here is a list of the submodules and short description of what they contain. - - - :any:`ot.lp` contains OT solvers for the exact (Linear Program) OT problems. - - :any:`ot.bregman` contains OT solvers for the entropic OT problems using - Bregman projections. - - :any:`ot.lp` contains OT solvers for the exact (Linear Program) OT problems. - - :any:`ot.smooth` contains OT solvers for the regularized (l2 and kl) smooth OT - problems. - - :any:`ot.gromov` contains solvers for Gromov-Wasserstein and Fused Gromov - Wasserstein problems. - - :any:`ot.optim` contains generic solvers OT based optimization problems - - :any:`ot.da` contains classes and function related to Monge mapping - estimation and Domain Adaptation (DA). - - :any:`ot.gpu` contains GPU (cupy) implementation of some OT solvers - - :any:`ot.dr` contains Dimension Reduction (DR) methods such as Wasserstein - Discriminant Analysis. - - :any:`ot.utils` contains utility functions such as distance computation and - timing. - - :any:`ot.datasets` contains toy dataset generation functions. - - :any:`ot.plot` contains visualization functions - - :any:`ot.stochastic` contains stochastic solvers for regularized OT. - - :any:`ot.unbalanced` contains solvers for regularized unbalanced OT. - - :any:`ot.partial` contains solvers for partial OT. - .. warning:: The list of automatically imported sub-modules is as follows: :py:mod:`ot.lp`, :py:mod:`ot.bregman`, :py:mod:`ot.optim` diff --git a/ot/bregman.py b/ot/bregman.py index b4365d0..f1f8437 100644 --- a/ot/bregman.py +++ b/ot/bregman.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- """ -Bregman projections for regularized OT +Bregman projections solvers for entropic regularized OT """ # Author: Remi Flamary diff --git a/ot/datasets.py b/ot/datasets.py index daca1ae..bd39e13 100644 --- a/ot/datasets.py +++ b/ot/datasets.py @@ -1,5 +1,5 @@ """ -Simple example datasets for OT +Simple example datasets """ # Author: Remi Flamary diff --git a/ot/dr.py b/ot/dr.py index 680dabf..11d2e10 100644 --- a/ot/dr.py +++ b/ot/dr.py @@ -1,10 +1,10 @@ # -*- coding: utf-8 -*- """ -Dimension reduction with optimal transport +Dimension reduction with OT .. warning:: - Note that by default the module is not import in :mod:`ot`. In order to + Note that by default the module is not imported in :mod:`ot`. In order to use it you need to explicitely import :mod:`ot.dr` """ diff --git a/ot/gpu/__init__.py b/ot/gpu/__init__.py index 1ab95bb..7478fb9 100644 --- a/ot/gpu/__init__.py +++ b/ot/gpu/__init__.py @@ -1,8 +1,9 @@ # -*- coding: utf-8 -*- """ +GPU implementation for several OT solvers and utility +functions. -This module provides GPU implementation for several OT solvers and utility -functions. The GPU backend in handled by `cupy +The GPU backend in handled by `cupy `_. .. warning:: diff --git a/ot/gromov.py b/ot/gromov.py index 43780a4..a678722 100644 --- a/ot/gromov.py +++ b/ot/gromov.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- """ -Gromov-Wasserstein transport method +Gromov-Wasserstein and Fused-Gromov-Wasserstein solvers """ # Author: Erwan Vautier diff --git a/ot/optim.py b/ot/optim.py index 4012e0d..b9ca891 100644 --- a/ot/optim.py +++ b/ot/optim.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- """ -Optimization algorithms for OT +Generic solvers for regularized OT """ # Author: Remi Flamary diff --git a/ot/partial.py b/ot/partial.py index c03ec25..eb707d8 100755 --- a/ot/partial.py +++ b/ot/partial.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- """ -Partial OT +Partial OT solvers """ # Author: Laetitia Chapel diff --git a/ot/smooth.py b/ot/smooth.py index 5a8e4b5..81f6a3e 100644 --- a/ot/smooth.py +++ b/ot/smooth.py @@ -26,7 +26,9 @@ # Remi Flamary """ -Implementation of +Smooth and Sparse Optimal Transport solvers (KL an L2 reg.) + +Implementation of : Smooth and Sparse Optimal Transport. Mathieu Blondel, Vivien Seguy, Antoine Rolet. In Proc. of AISTATS 2018. diff --git a/ot/unbalanced.py b/ot/unbalanced.py index 23f6607..e37f10c 100644 --- a/ot/unbalanced.py +++ b/ot/unbalanced.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- """ -Regularized Unbalanced OT +Regularized Unbalanced OT solvers """ # Author: Hicham Janati -- cgit v1.2.3