diff options
| -rw-r--r-- | examples/plot_otda_laplacian.py | 127 | ||||
| -rw-r--r-- | ot/da.py | 257 | ||||
| -rw-r--r-- | ot/utils.py | 6 | ||||
| -rw-r--r-- | test/test_da.py | 64 |
4 files changed, 453 insertions, 1 deletions
diff --git a/examples/plot_otda_laplacian.py b/examples/plot_otda_laplacian.py new file mode 100644 index 0000000..67c8f67 --- /dev/null +++ b/examples/plot_otda_laplacian.py @@ -0,0 +1,127 @@ +# -*- coding: utf-8 -*- +""" +====================================================== +OT with Laplacian regularization for domain adaptation +====================================================== + +This example introduces a domain adaptation in a 2D setting and OTDA +approach with Laplacian regularization. + +""" + +# Authors: Ievgen Redko <ievgen.redko@univ-st-etienne.fr> + +# License: MIT License + +import matplotlib.pylab as pl +import ot + +############################################################################## +# Generate data +# ------------- + +n_source_samples = 150 +n_target_samples = 150 + +Xs, ys = ot.datasets.make_data_classif('3gauss', n_source_samples) +Xt, yt = ot.datasets.make_data_classif('3gauss2', n_target_samples) + + +############################################################################## +# Instantiate the different transport algorithms and fit them +# ----------------------------------------------------------- + +# EMD Transport +ot_emd = ot.da.EMDTransport() +ot_emd.fit(Xs=Xs, Xt=Xt) + +# Sinkhorn Transport +ot_sinkhorn = ot.da.SinkhornTransport(reg_e=.01) +ot_sinkhorn.fit(Xs=Xs, Xt=Xt) + +# EMD Transport with Laplacian regularization +ot_emd_laplace = ot.da.EMDLaplaceTransport(reg_lap=100, reg_src=1) +ot_emd_laplace.fit(Xs=Xs, Xt=Xt) + +# transport source samples onto target samples +transp_Xs_emd = ot_emd.transform(Xs=Xs) +transp_Xs_sinkhorn = ot_sinkhorn.transform(Xs=Xs) +transp_Xs_emd_laplace = ot_emd_laplace.transform(Xs=Xs) + +############################################################################## +# Fig 1 : plots source and target samples +# --------------------------------------- + +pl.figure(1, figsize=(10, 5)) +pl.subplot(1, 2, 1) +pl.scatter(Xs[:, 0], Xs[:, 1], c=ys, marker='+', label='Source samples') +pl.xticks([]) +pl.yticks([]) +pl.legend(loc=0) +pl.title('Source samples') + +pl.subplot(1, 2, 2) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', label='Target samples') +pl.xticks([]) +pl.yticks([]) +pl.legend(loc=0) +pl.title('Target samples') +pl.tight_layout() + + +############################################################################## +# Fig 2 : plot optimal couplings and transported samples +# ------------------------------------------------------ + +param_img = {'interpolation': 'nearest'} + +pl.figure(2, figsize=(15, 8)) +pl.subplot(2, 3, 1) +pl.imshow(ot_emd.coupling_, **param_img) +pl.xticks([]) +pl.yticks([]) +pl.title('Optimal coupling\nEMDTransport') + +pl.figure(2, figsize=(15, 8)) +pl.subplot(2, 3, 2) +pl.imshow(ot_sinkhorn.coupling_, **param_img) +pl.xticks([]) +pl.yticks([]) +pl.title('Optimal coupling\nSinkhornTransport') + +pl.subplot(2, 3, 3) +pl.imshow(ot_emd_laplace.coupling_, **param_img) +pl.xticks([]) +pl.yticks([]) +pl.title('Optimal coupling\nEMDLaplaceTransport') + +pl.subplot(2, 3, 4) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=0.3) +pl.scatter(transp_Xs_emd[:, 0], transp_Xs_emd[:, 1], c=ys, + marker='+', label='Transp samples', s=30) +pl.xticks([]) +pl.yticks([]) +pl.title('Transported samples\nEmdTransport') +pl.legend(loc="lower left") + +pl.subplot(2, 3, 5) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=0.3) +pl.scatter(transp_Xs_sinkhorn[:, 0], transp_Xs_sinkhorn[:, 1], c=ys, + marker='+', label='Transp samples', s=30) +pl.xticks([]) +pl.yticks([]) +pl.title('Transported samples\nSinkhornTransport') + +pl.subplot(2, 3, 6) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=0.3) +pl.scatter(transp_Xs_emd_laplace[:, 0], transp_Xs_emd_laplace[:, 1], c=ys, + marker='+', label='Transp samples', s=30) +pl.xticks([]) +pl.yticks([]) +pl.title('Transported samples\nEMDLaplaceTransport') +pl.tight_layout() + +pl.show() @@ -16,7 +16,7 @@ import scipy.linalg as linalg from .bregman import sinkhorn, jcpot_barycenter from .lp import emd -from .utils import unif, dist, kernel, cost_normalization, label_normalization +from .utils import unif, dist, kernel, cost_normalization, label_normalization, laplacian, dots from .utils import check_params, BaseEstimator from .unbalanced import sinkhorn_unbalanced from .optim import cg @@ -748,6 +748,139 @@ def OT_mapping_linear(xs, xt, reg=1e-6, ws=None, return A, b +def emd_laplace(a, b, xs, xt, M, reg, eta, alpha, + numItermax, stopThr, numInnerItermax, + stopInnerThr, log=False, verbose=False, **kwargs): + r"""Solve the optimal transport problem (OT) with Laplacian regularization + + .. math:: + \gamma = arg\min_\gamma <\gamma,M>_F + eta\Omega_\alpha(\gamma) + + s.t.\ \gamma 1 = a + + \gamma^T 1= b + + \gamma\geq 0 + + where: + + - a and b are source and target weights (sum to 1) + - xs and xt are source and target samples + - M is the (ns,nt) metric cost matrix + - :math:`\Omega_\alpha` is the Laplacian regularization term + :math:`\Omega_\alpha = (1-\alpha)/n_s^2\sum_{i,j}S^s_{i,j}\|T(\mathbf{x}^s_i)-T(\mathbf{x}^s_j)\|^2+\alpha/n_t^2\sum_{i,j}S^t_{i,j}^'\|T(\mathbf{x}^t_i)-T(\mathbf{x}^t_j)\|^2` + with :math:`S^s_{i,j}, S^t_{i,j}` denoting source and target similarity matrices and :math:`T(\cdot)` being a barycentric mapping + + The algorithm used for solving the problem is the conditional gradient algorithm as proposed in [5]. + + Parameters + ---------- + a : np.ndarray (ns,) + samples weights in the source domain + b : np.ndarray (nt,) + samples weights in the target domain + xs : np.ndarray (ns,d) + samples in the source domain + xt : np.ndarray (nt,d) + samples in the target domain + M : np.ndarray (ns,nt) + loss matrix + reg : string + Type of Laplacian regularization + eta : float + Regularization term for Laplacian regularization + alpha : float + Regularization term for source domain's importance in regularization + numItermax : int, optional + Max number of iterations + stopThr : float, optional + Stop threshold on error (inner emd solver) (>0) + numInnerItermax : int, optional + Max number of iterations (inner CG solver) + stopInnerThr : float, optional + Stop threshold on error (inner CG solver) (>0) + verbose : bool, optional + Print information along iterations + log : bool, optional + record log if True + kwargs : dict + Dictionary with attributes 'sim' ('knn' or 'gauss') and + 'param' (int, float or None) for similarity type and its parameter to be used. + If 'param' is None, it is computed as mean pairwise Euclidean distance over the data set + or set to 3 when sim is 'gauss' or 'knn', respectively. + + Returns + ------- + gamma : (ns x nt) ndarray + Optimal transportation matrix for the given parameters + log : dict + log dictionary return only if log==True in parameters + + + References + ---------- + + .. [5] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, + "Optimal Transport for Domain Adaptation," in IEEE + Transactions on Pattern Analysis and Machine Intelligence , + vol.PP, no.99, pp.1-1 + .. [28] R. Flamary, N. Courty, D. Tuia, A. Rakotomamonjy, + "Optimal transport with Laplacian regularization: Applications to domain adaptation and shape matching," + in NIPS Workshop on Optimal Transport and Machine Learning OTML, 2014. + + See Also + -------- + ot.lp.emd : Unregularized OT + ot.optim.cg : General regularized OT + + """ + if not isinstance(kwargs['param'], (int, float, type(None))): + raise ValueError( + 'Similarity parameter should be an int or a float. Got {type} instead.'.format(type=type(kwargs['param']))) + + if kwargs['sim'] == 'gauss': + if kwargs['param'] is None: + kwargs['param'] = 1 / (2 * (np.mean(dist(xs, xs, 'sqeuclidean')) ** 2)) + sS = kernel(xs, xs, method=kwargs['sim'], sigma=kwargs['param']) + sT = kernel(xt, xt, method=kwargs['sim'], sigma=kwargs['param']) + + elif kwargs['sim'] == 'knn': + if kwargs['param'] is None: + kwargs['param'] = 3 + + from sklearn.neighbors import kneighbors_graph + + sS = kneighbors_graph(X=xs, n_neighbors=int(kwargs['param'])).toarray() + sS = (sS + sS.T) / 2 + sT = kneighbors_graph(xt, n_neighbors=int(kwargs['param'])).toarray() + sT = (sT + sT.T) / 2 + else: + raise ValueError('Unknown similarity type {sim}. Currently supported similarity types are "knn" and "gauss".'.format(sim=kwargs['sim'])) + + lS = laplacian(sS) + lT = laplacian(sT) + + def f(G): + return alpha * np.trace(np.dot(xt.T, np.dot(G.T, np.dot(lS, np.dot(G, xt))))) \ + + (1 - alpha) * np.trace(np.dot(xs.T, np.dot(G, np.dot(lT, np.dot(G.T, xs))))) + + ls2 = lS + lS.T + lt2 = lT + lT.T + xt2 = np.dot(xt, xt.T) + + if reg == 'disp': + Cs = -eta * alpha / xs.shape[0] * dots(ls2, xs, xt.T) + Ct = -eta * (1 - alpha) / xt.shape[0] * dots(xs, xt.T, lt2) + M = M + Cs + Ct + + def df(G): + return alpha * np.dot(ls2, np.dot(G, xt2))\ + + (1 - alpha) * np.dot(xs, np.dot(xs.T, np.dot(G, lt2))) + + return cg(a, b, M, reg=eta, f=f, df=df, G0=None, numItermax=numItermax, numItermaxEmd=numInnerItermax, + stopThr=stopThr, stopThr2=stopInnerThr, verbose=verbose, log=log) + + def distribution_estimation_uniform(X): """estimates a uniform distribution from an array of samples X @@ -1576,6 +1709,128 @@ class SinkhornLpl1Transport(BaseTransport): return self +class EMDLaplaceTransport(BaseTransport): + + """Domain Adapatation OT method based on Earth Mover's Distance with Laplacian regularization + + Parameters + ---------- + reg_type : string optional (default='pos') + Type of the regularization term: 'pos' and 'disp' for + regularization term defined in [2] and [6], respectively. + reg_lap : float, optional (default=1) + Laplacian regularization parameter + reg_src : float, optional (default=0.5) + Source relative importance in regularization + metric : string, optional (default="sqeuclidean") + The ground metric for the Wasserstein problem + norm : string, optional (default=None) + If given, normalize the ground metric to avoid numerical errors that + can occur with large metric values. + similarity : string, optional (default="knn") + The similarity to use either knn or gaussian + similarity_param : int or float, optional (default=3) + Parameter for the similarity: number of nearest neighbors or bandwidth + if similarity="knn" or "gaussian", respectively. + max_iter : int, optional (default=100) + Max number of BCD iterations + tol : float, optional (default=1e-5) + Stop threshold on relative loss decrease (>0) + max_inner_iter : int, optional (default=10) + Max number of iterations (inner CG solver) + inner_tol : float, optional (default=1e-6) + Stop threshold on error (inner CG solver) (>0) + log : int, optional (default=False) + Controls the logs of the optimization algorithm + distribution_estimation : callable, optional (defaults to the uniform) + The kind of distribution estimation to employ + out_of_sample_map : string, optional (default="ferradans") + The kind of out of sample mapping to apply to transport samples + from a domain into another one. Currently the only possible option is + "ferradans" which uses the method proposed in [6]. + + Attributes + ---------- + coupling_ : array-like, shape (n_source_samples, n_target_samples) + The optimal coupling + + References + ---------- + .. [1] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, + "Optimal Transport for Domain Adaptation," in IEEE Transactions + on Pattern Analysis and Machine Intelligence , vol.PP, no.99, pp.1-1 + .. [2] R. Flamary, N. Courty, D. Tuia, A. Rakotomamonjy, + "Optimal transport with Laplacian regularization: Applications to domain adaptation and shape matching," + in NIPS Workshop on Optimal Transport and Machine Learning OTML, 2014. + """ + + def __init__(self, reg_type='pos', reg_lap=1., reg_src=1., alpha=0.5, + metric="sqeuclidean", norm=None, similarity="knn", similarity_param=None, max_iter=100, tol=1e-9, + max_inner_iter=100000, inner_tol=1e-9, log=False, verbose=False, + distribution_estimation=distribution_estimation_uniform, + out_of_sample_map='ferradans'): + self.reg = reg_type + self.reg_lap = reg_lap + self.reg_src = reg_src + self.alpha = alpha + self.metric = metric + self.norm = norm + self.similarity = similarity + self.sim_param = similarity_param + self.max_iter = max_iter + self.tol = tol + self.max_inner_iter = max_inner_iter + self.inner_tol = inner_tol + self.log = log + self.verbose = verbose + self.distribution_estimation = distribution_estimation + self.out_of_sample_map = out_of_sample_map + + def fit(self, Xs, ys=None, Xt=None, yt=None): + """Build a coupling matrix from source and target sets of samples + (Xs, ys) and (Xt, yt) + + Parameters + ---------- + Xs : array-like, shape (n_source_samples, n_features) + The training input samples. + ys : array-like, shape (n_source_samples,) + The class labels + Xt : array-like, shape (n_target_samples, n_features) + The training input samples. + yt : array-like, shape (n_target_samples,) + The class labels. If some target samples are unlabeled, fill the + yt's elements with -1. + + Warning: Note that, due to this convention -1 cannot be used as a + class label + + Returns + ------- + self : object + Returns self. + """ + + super(EMDLaplaceTransport, self).fit(Xs, ys, Xt, yt) + + kwargs = dict() + kwargs["sim"] = self.similarity + kwargs["param"] = self.sim_param + + returned_ = emd_laplace(a=self.mu_s, b=self.mu_t, xs=self.xs_, + xt=self.xt_, M=self.cost_, reg=self.reg, eta=self.reg_lap, alpha=self.reg_src, + numItermax=self.max_iter, stopThr=self.tol, numInnerItermax=self.max_inner_iter, + stopInnerThr=self.inner_tol, log=self.log, verbose=self.verbose, **kwargs) + + # coupling estimation + if self.log: + self.coupling_, self.log_ = returned_ + else: + self.coupling_ = returned_ + self.log_ = dict() + return self + + class SinkhornL1l2Transport(BaseTransport): """Domain Adapatation OT method based on sinkhorn algorithm + diff --git a/ot/utils.py b/ot/utils.py index c154f99..f9911a1 100644 --- a/ot/utils.py +++ b/ot/utils.py @@ -49,6 +49,12 @@ def kernel(x1, x2, method='gaussian', sigma=1, **kwargs): return K +def laplacian(x): + """Compute Laplacian matrix""" + L = np.diag(np.sum(x, axis=0)) - x + return L + + def unif(n): """ return a uniform histogram of length n (simplex) diff --git a/test/test_da.py b/test/test_da.py index 7d0fdda..3b28119 100644 --- a/test/test_da.py +++ b/test/test_da.py @@ -689,3 +689,67 @@ def test_jcpot_barycenter(): numItermax=10000, stopThr=1e-9, verbose=False, log=False) np.testing.assert_allclose(prop, [1 - pt, pt], rtol=1e-3, atol=1e-3) + + +def test_emd_laplace_class(): + """test_emd_laplace_transport + """ + ns = 150 + nt = 200 + + Xs, ys = make_data_classif('3gauss', ns) + Xt, yt = make_data_classif('3gauss2', nt) + + otda = ot.da.EMDLaplaceTransport(reg_lap=0.01, max_iter=1000, tol=1e-9, verbose=False, log=True) + + # test its computed + otda.fit(Xs=Xs, ys=ys, Xt=Xt) + + assert hasattr(otda, "coupling_") + assert hasattr(otda, "log_") + + # test dimensions of coupling + assert_equal(otda.coupling_.shape, ((Xs.shape[0], Xt.shape[0]))) + + # test all margin constraints + mu_s = unif(ns) + mu_t = unif(nt) + + assert_allclose( + np.sum(otda.coupling_, axis=0), mu_t, rtol=1e-3, atol=1e-3) + assert_allclose( + np.sum(otda.coupling_, axis=1), mu_s, rtol=1e-3, atol=1e-3) + + # test transform + transp_Xs = otda.transform(Xs=Xs) + [assert_equal(x.shape, y.shape) for x, y in zip(transp_Xs, Xs)] + + Xs_new, _ = make_data_classif('3gauss', ns + 1) + transp_Xs_new = otda.transform(Xs_new) + + # check that the oos method is working + assert_equal(transp_Xs_new.shape, Xs_new.shape) + + # test inverse transform + transp_Xt = otda.inverse_transform(Xt=Xt) + assert_equal(transp_Xt.shape, Xt.shape) + + Xt_new, _ = make_data_classif('3gauss2', nt + 1) + transp_Xt_new = otda.inverse_transform(Xt=Xt_new) + + # check that the oos method is working + assert_equal(transp_Xt_new.shape, Xt_new.shape) + + # test fit_transform + transp_Xs = otda.fit_transform(Xs=Xs, Xt=Xt) + assert_equal(transp_Xs.shape, Xs.shape) + + # check label propagation + transp_yt = otda.transform_labels(ys) + assert_equal(transp_yt.shape[0], yt.shape[0]) + assert_equal(transp_yt.shape[1], len(np.unique(ys))) + + # check inverse label propagation + transp_ys = otda.inverse_transform_labels(yt) + assert_equal(transp_ys.shape[0], ys.shape[0]) + assert_equal(transp_ys.shape[1], len(np.unique(yt))) |