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.
- 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.