ot.lp

Solvers for the original linear program OT problem

Functions

ot.lp.center_ot_dual(alpha0, beta0, a=None, b=None)[source]

Center dual OT potentials w.r.t. their weights

The main idea of this function is to find unique dual potentials that ensure some kind of centering/fairness. The main idea is to find dual potentials that lead to the same final objective value for both source and targets (see below for more details). It will help having stability when multiple calling of the OT solver with small changes.

Basically we add another constraint to the potential that will not change the objective value but will ensure unicity. The constraint is the following:

\[\alpha^T \mathbf{a} = \beta^T \mathbf{b}\]

in addition to the OT problem constraints.

since \(\sum_i a_i=\sum_j b_j\) this can be solved by adding/removing a constant from both \(\alpha_0\) and \(\beta_0\).

\[ \begin{align}\begin{aligned}c &= \frac{\beta_0^T \mathbf{b} - \alpha_0^T \mathbf{a}}{\mathbf{1}^T \mathbf{b} + \mathbf{1}^T \mathbf{a}}\\\alpha &= \alpha_0 + c\\\beta &= \beta_0 + c\end{aligned}\end{align} \]
Parameters
  • alpha0 ((ns,) numpy.ndarray, float64) – Source dual potential

  • beta0 ((nt,) numpy.ndarray, float64) – Target dual potential

  • a ((ns,) numpy.ndarray, float64) – Source histogram (uniform weight if empty list)

  • b ((nt,) numpy.ndarray, float64) – Target histogram (uniform weight if empty list)

Returns

  • alpha ((ns,) numpy.ndarray, float64) – Source centered dual potential

  • beta ((nt,) numpy.ndarray, float64) – Target centered dual potential

ot.lp.check_number_threads(numThreads)[source]

Checks whether or not the requested number of threads has a valid value.

Parameters

numThreads (int or str) – The requested number of threads, should either be a strictly positive integer or “max” or None

Returns

numThreads – Corrected number of threads

Return type

int

ot.lp.emd(a, b, M, numItermax=100000, log=False, center_dual=True, numThreads=1)[source]

Solves the Earth Movers distance problem and returns the OT matrix

\[ \begin{align}\begin{aligned}\gamma = \mathop{\arg \min}_\gamma \quad \langle \gamma, \mathbf{M} \rangle_F\\s.t. \ \gamma \mathbf{1} = \mathbf{a}\\ \gamma^T \mathbf{1} = \mathbf{b}\\ \gamma \geq 0\end{aligned}\end{align} \]

where :

  • \(\mathbf{M}\) is the metric cost matrix

  • \(\mathbf{a}\) and \(\mathbf{b}\) are the sample weights

Warning

Note that the \(\mathbf{M}\) matrix in numpy needs to be a C-order numpy.array in float64 format. It will be converted if not in this format

Note

This function is backend-compatible and will work on arrays from all compatible backends. But the algorithm uses the C++ CPU backend which can lead to copy overhead on GPU arrays.

Note

This function will cast the computed transport plan to the data type of the provided input with the following priority: \(\mathbf{a}\), then \(\mathbf{b}\), then \(\mathbf{M}\) if marginals are not provided. Casting to an integer tensor might result in a loss of precision. If this behaviour is unwanted, please make sure to provide a floating point input.

Uses the algorithm proposed in [1].

Parameters
  • a ((ns,) array-like, float) – Source histogram (uniform weight if empty list)

  • b ((nt,) array-like, float) – Target histogram (uniform weight if empty list)

  • M ((ns,nt) array-like, float) – Loss matrix (c-order array in numpy with type float64)

  • numItermax (int, optional (default=100000)) – The maximum number of iterations before stopping the optimization algorithm if it has not converged.

  • log (bool, optional (default=False)) – If True, returns a dictionary containing the cost and dual variables. Otherwise returns only the optimal transportation matrix.

  • center_dual (boolean, optional (default=True)) – If True, centers the dual potential using function center_ot_dual.

  • numThreads (int or "max", optional (default=1, i.e. OpenMP is not used)) – If compiled with OpenMP, chooses the number of threads to parallelize. “max” selects the highest number possible.

Returns

  • gamma (array-like, shape (ns, nt)) – Optimal transportation matrix for the given parameters

  • log (dict, optional) – If input log is true, a dictionary containing the cost and dual variables and exit status

Examples

Simple example with obvious solution. The function emd accepts lists and perform automatic conversion to numpy arrays

>>> import ot
>>> a=[.5,.5]
>>> b=[.5,.5]
>>> M=[[0.,1.],[1.,0.]]
>>> ot.emd(a, b, M)
array([[0.5, 0. ],
       [0. , 0.5]])

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.

See also

ot.bregman.sinkhorn

Entropic regularized OT

ot.optim.cg

General regularized OT

ot.lp.emd2(a, b, M, processes=1, numItermax=100000, log=False, return_matrix=False, center_dual=True, numThreads=1)[source]

Solves the Earth Movers distance problem and returns the loss

\[ \begin{align}\begin{aligned}\min_\gamma \quad \langle \gamma, \mathbf{M} \rangle_F\\s.t. \ \gamma \mathbf{1} = \mathbf{a}\\ \gamma^T \mathbf{1} = \mathbf{b}\\ \gamma \geq 0\end{aligned}\end{align} \]

where :

  • \(\mathbf{M}\) is the metric cost matrix

  • \(\mathbf{a}\) and \(\mathbf{b}\) are the sample weights

Note

This function is backend-compatible and will work on arrays from all compatible backends. But the algorithm uses the C++ CPU backend which can lead to copy overhead on GPU arrays.

Note

This function will cast the computed transport plan and transportation loss to the data type of the provided input with the following priority: \(\mathbf{a}\), then \(\mathbf{b}\), then \(\mathbf{M}\) if marginals are not provided. Casting to an integer tensor might result in a loss of precision. If this behaviour is unwanted, please make sure to provide a floating point input.

Uses the algorithm proposed in [1].

Parameters
  • a ((ns,) array-like, float64) – Source histogram (uniform weight if empty list)

  • b ((nt,) array-like, float64) – Target histogram (uniform weight if empty list)

  • M ((ns,nt) array-like, float64) – Loss matrix (for numpy c-order array with type float64)

  • processes (int, optional (default=1)) – Nb of processes used for multiple emd computation (deprecated)

  • numItermax (int, optional (default=100000)) – The maximum number of iterations before stopping the optimization algorithm if it has not converged.

  • log (boolean, optional (default=False)) – If True, returns a dictionary containing dual variables. Otherwise returns only the optimal transportation cost.

  • return_matrix (boolean, optional (default=False)) – If True, returns the optimal transportation matrix in the log.

  • center_dual (boolean, optional (default=True)) – If True, centers the dual potential using function center_ot_dual.

  • numThreads (int or "max", optional (default=1, i.e. OpenMP is not used)) – If compiled with OpenMP, chooses the number of threads to parallelize. “max” selects the highest number possible.

Returns

  • W (float, array-like) – Optimal transportation loss for the given parameters

  • log (dict) – If input log is true, a dictionary containing dual variables and exit status

Examples

Simple example with obvious solution. The function emd accepts lists and perform automatic conversion to numpy arrays

>>> import ot
>>> a=[.5,.5]
>>> b=[.5,.5]
>>> M=[[0.,1.],[1.,0.]]
>>> ot.emd2(a,b,M)
0.0

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.

See also

ot.bregman.sinkhorn

Entropic regularized OT

ot.optim.cg

General regularized OT

ot.lp.estimate_dual_null_weights(alpha0, beta0, a, b, M)[source]

Estimate feasible values for 0-weighted dual potentials

The feasible values are computed efficiently but rather coarsely.

Warning

This function is necessary because the C++ solver in emd_c discards all samples in the distributions with zeros weights. This means that while the primal variable (transport matrix) is exact, the solver only returns feasible dual potentials on the samples with weights different from zero.

First we compute the constraints violations:

\[\mathbf{V} = \alpha + \beta^T - \mathbf{M}\]

Next we compute the max amount of violation per row (\(\alpha\)) and columns (\(beta\))

\[ \begin{align}\begin{aligned}\mathbf{v^a}_i = \max_j \mathbf{V}_{i,j}\\\mathbf{v^b}_j = \max_i \mathbf{V}_{i,j}\end{aligned}\end{align} \]

Finally we update the dual potential with 0 weights if a constraint is violated

\[ \begin{align}\begin{aligned}\alpha_i = \alpha_i - \mathbf{v^a}_i \quad \text{ if } \mathbf{a}_i=0 \text{ and } \mathbf{v^a}_i>0\\\beta_j = \beta_j - \mathbf{v^b}_j \quad \text{ if } \mathbf{b}_j=0 \text{ and } \mathbf{v^b}_j > 0\end{aligned}\end{align} \]

In the end the dual potentials are centered using function ot.lp.center_ot_dual().

Note that all those updates do not change the objective value of the solution but provide dual potentials that do not violate the constraints.

Parameters
  • alpha0 ((ns,) numpy.ndarray, float64) – Source dual potential

  • beta0 ((nt,) numpy.ndarray, float64) – Target dual potential

  • alpha0 – Source dual potential

  • beta0 – Target dual potential

  • a ((ns,) numpy.ndarray, float64) – Source distribution (uniform weights if empty list)

  • b ((nt,) numpy.ndarray, float64) – Target distribution (uniform weights if empty list)

  • M ((ns,nt) numpy.ndarray, float64) – Loss matrix (c-order array with type float64)

Returns

  • alpha ((ns,) numpy.ndarray, float64) – Source corrected dual potential

  • beta ((nt,) numpy.ndarray, float64) – Target corrected dual potential

ot.lp.free_support_barycenter(measures_locations, measures_weights, X_init, b=None, weights=None, numItermax=100, stopThr=1e-07, verbose=False, log=None, numThreads=1)[source]

Solves the free support (locations of the barycenters are optimized, not the weights) Wasserstein barycenter problem (i.e. the weighted Frechet mean for the 2-Wasserstein distance), formally:

\[\min_\mathbf{X} \quad \sum_{i=1}^N w_i W_2^2(\mathbf{b}, \mathbf{X}, \mathbf{a}_i, \mathbf{X}_i)\]

where :

  • \(w \in \mathbb{(0, 1)}^{N}\)’s are the barycenter weights and sum to one

  • the \(\mathbf{a}_i \in \mathbb{R}^{k_i}\) are the empirical measures weights and sum to one for each \(i\)

  • the \(\mathbf{X}_i \in \mathbb{R}^{k_i, d}\) are the empirical measures atoms locations

  • \(\mathbf{b} \in \mathbb{R}^{k}\) is the desired weights vector of the barycenter

This problem is considered in [1] (Algorithm 2). There are two differences with the following codes:

  • we do not optimize over the weights

  • we do not do line search for the locations updates, we use i.e. \(\theta = 1\) in [1] (Algorithm 2). This can be seen as a discrete implementation of the fixed-point algorithm of [2] proposed in the continuous setting.

Parameters
  • measures_locations (list of N (k_i,d) array-like) – The discrete support of a measure supported on \(k_i\) locations of a d-dimensional space (\(k_i\) can be different for each element of the list)

  • measures_weights (list of N (k_i,) array-like) – Numpy arrays where each numpy array has \(k_i\) non-negatives values summing to one representing the weights of each discrete input measure

  • X_init ((k,d) array-like) – Initialization of the support locations (on k atoms) of the barycenter

  • b ((k,) array-like) – Initialization of the weights of the barycenter (non-negatives, sum to 1)

  • weights ((N,) array-like) – Initialization of the coefficients of the barycenter (non-negatives, sum to 1)

  • 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

  • numThreads (int or "max", optional (default=1, i.e. OpenMP is not used)) – If compiled with OpenMP, chooses the number of threads to parallelize. “max” selects the highest number possible.

Returns

X – Support locations (on k atoms) of the barycenter

Return type

(k,d) array-like

References

1

Cuturi, Marco, and Arnaud Doucet. “Fast computation of Wasserstein barycenters.” International Conference on Machine Learning. 2014.

2

Álvarez-Esteban, Pedro C., et al. “A fixed-point approach to barycenters in Wasserstein space.” Journal of Mathematical Analysis and Applications 441.2 (2016): 744-762.

Examples using ot.lp.free_support_barycenter