Regularized OT with generic solver

Illustrates the use of the generic solver for regularized OT with user-designed regularization term. It uses Conditional gradient as in [6] and generalized Conditional Gradient as proposed in [5,7].

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

[6] Ferradans, S., Papadakis, N., Peyré, G., & Aujol, J. F. (2014). Regularized discrete optimal transport. SIAM Journal on Imaging Sciences, 7(3), 1853-1882.

[7] Rakotomamonjy, A., Flamary, R., & Courty, N. (2015). Generalized conditional gradient: analysis of convergence and applications. arXiv preprint arXiv:1510.06567.

# sphinx_gallery_thumbnail_number = 5

import numpy as np
import matplotlib.pylab as pl
import ot
import ot.plot

Generate data

n = 100  # nb bins

# bin positions
x = np.arange(n, dtype=np.float64)

# Gaussian distributions
a = ot.datasets.make_1D_gauss(n, m=20, s=5)  # m= mean, s= std
b = ot.datasets.make_1D_gauss(n, m=60, s=10)

# loss matrix
M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1)))
M /= M.max()

Solve EMD

G0 = ot.emd(a, b, M)

pl.figure(1, figsize=(5, 5))
ot.plot.plot1D_mat(a, b, G0, "OT matrix G0")

(<Axes: >, <Axes: >, <Axes: >)

Solve EMD with Frobenius norm regularization

def f(G):
    return 0.5 * np.sum(G**2)

def df(G):
    return G

reg = 1e-1

Gl2 = ot.optim.cg(a, b, M, reg, f, df, verbose=True)

pl.figure(2)
ot.plot.plot1D_mat(a, b, Gl2, "OT matrix Frob. reg")

It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.760578e-01|0.000000e+00|0.000000e+00
1.669467e-01|5.457501e-02|9.111116e-03
1.665639e-01|2.298097e-03|3.827801e-04
1.664378e-01|7.573707e-04|1.260551e-04
1.664075e-01|1.822159e-04|3.032210e-05
1.663910e-01|9.895078e-05|1.646452e-05
1.663851e-01|3.549058e-05|5.905106e-06
1.663814e-01|2.252945e-05|3.748482e-06
1.663785e-01|1.757545e-05|2.924177e-06
1.663767e-01|1.060202e-05|1.763930e-06
1.663751e-01|9.640161e-06|1.603883e-06
1.663737e-01|8.501133e-06|1.414365e-06
1.663727e-01|5.654674e-06|9.407836e-07
1.663720e-01|4.674546e-06|7.777135e-07
1.663712e-01|4.600984e-06|7.654712e-07
1.663707e-01|3.087272e-06|5.136316e-07
1.663702e-01|3.227426e-06|5.369474e-07
1.663696e-01|3.398875e-06|5.654694e-07
1.663691e-01|2.877154e-06|4.786695e-07
1.663687e-01|2.733590e-06|4.547836e-07
It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.663683e-01|2.417235e-06|4.021512e-07
1.663678e-01|2.668394e-06|4.439349e-07
1.663675e-01|2.063226e-06|3.432536e-07
1.663671e-01|2.165567e-06|3.602791e-07
1.663668e-01|2.102621e-06|3.498062e-07
1.663664e-01|2.004991e-06|3.335631e-07
1.663661e-01|2.095035e-06|3.485427e-07
1.663658e-01|1.810346e-06|3.011796e-07
1.663655e-01|1.730306e-06|2.878632e-07
1.663652e-01|1.691085e-06|2.813377e-07
1.663649e-01|1.712194e-06|2.848491e-07
1.663647e-01|1.599027e-06|2.660216e-07
1.663644e-01|1.537822e-06|2.558388e-07
1.663641e-01|1.578779e-06|2.626522e-07
1.663639e-01|1.196930e-06|1.991260e-07
1.663637e-01|1.479310e-06|2.461035e-07
1.663635e-01|1.298958e-06|2.160992e-07
1.663633e-01|1.226052e-06|2.039700e-07
1.663630e-01|1.407332e-06|2.341280e-07
1.663628e-01|1.121523e-06|1.865798e-07
It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.663627e-01|1.025285e-06|1.705692e-07
1.663625e-01|9.811935e-07|1.632338e-07
1.663623e-01|1.061869e-06|1.766550e-07
1.663622e-01|1.088559e-06|1.810950e-07
1.663620e-01|9.551075e-07|1.588936e-07
1.663618e-01|1.014606e-06|1.687917e-07
1.663617e-01|1.015650e-06|1.689652e-07
1.663615e-01|8.803835e-07|1.464619e-07
1.663614e-01|8.544397e-07|1.421458e-07
1.663612e-01|7.635494e-07|1.270250e-07
1.663611e-01|7.780657e-07|1.294399e-07
1.663610e-01|7.944753e-07|1.321697e-07
1.663609e-01|7.184089e-07|1.195151e-07
1.663607e-01|8.663411e-07|1.441251e-07
1.663606e-01|7.813255e-07|1.299818e-07
1.663605e-01|7.969488e-07|1.325808e-07
1.663603e-01|7.436876e-07|1.237201e-07
1.663602e-01|7.143391e-07|1.188376e-07
1.663601e-01|6.864534e-07|1.141985e-07
1.663600e-01|6.857649e-07|1.140838e-07
It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.663599e-01|7.077854e-07|1.177471e-07
1.663598e-01|6.728464e-07|1.119346e-07
1.663597e-01|6.234005e-07|1.037087e-07
1.663596e-01|5.925432e-07|9.857522e-08
1.663595e-01|5.223502e-07|8.689791e-08
1.663594e-01|5.685675e-07|9.458653e-08
1.663593e-01|5.089579e-07|8.466987e-08
1.663592e-01|6.008686e-07|9.996001e-08
1.663591e-01|5.507201e-07|9.161730e-08
1.663590e-01|5.377172e-07|8.945410e-08
1.663589e-01|4.843623e-07|8.057800e-08
1.663588e-01|4.893034e-07|8.139994e-08
1.663588e-01|4.561640e-07|7.588689e-08
1.663587e-01|4.885435e-07|8.127345e-08
1.663586e-01|4.125727e-07|6.863502e-08
1.663585e-01|4.081210e-07|6.789442e-08
1.663585e-01|4.562190e-07|7.589589e-08
1.663584e-01|5.554839e-07|9.240940e-08
1.663583e-01|4.948121e-07|8.231610e-08
1.663582e-01|4.499924e-07|7.485994e-08
It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.663582e-01|4.292707e-07|7.141268e-08
1.663581e-01|4.076378e-07|6.781384e-08
1.663580e-01|3.436588e-07|5.717040e-08
1.663580e-01|4.070341e-07|6.771337e-08
1.663579e-01|3.680297e-07|6.122465e-08
1.663578e-01|3.546290e-07|5.899532e-08
1.663578e-01|4.509493e-07|7.501891e-08
1.663577e-01|3.720299e-07|6.189004e-08
1.663576e-01|4.068402e-07|6.768098e-08
1.663576e-01|3.074031e-07|5.113883e-08
1.663575e-01|3.748969e-07|6.236692e-08
1.663575e-01|3.451085e-07|5.741138e-08
1.663574e-01|3.413178e-07|5.678075e-08
1.663574e-01|3.359396e-07|5.588603e-08
1.663573e-01|3.291271e-07|5.475270e-08
1.663572e-01|3.293071e-07|5.478261e-08
1.663572e-01|3.094564e-07|5.148030e-08
1.663571e-01|2.921218e-07|4.859654e-08
1.663571e-01|2.648265e-07|4.405576e-08
1.663570e-01|3.304197e-07|5.496765e-08
It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.663570e-01|3.205966e-07|5.333349e-08
1.663569e-01|2.957885e-07|4.920646e-08
1.663569e-01|2.603218e-07|4.330633e-08
1.663569e-01|2.644096e-07|4.398636e-08
1.663568e-01|2.868374e-07|4.771735e-08
1.663568e-01|3.023438e-07|5.029694e-08
1.663567e-01|2.363974e-07|3.932629e-08
1.663567e-01|2.821861e-07|4.694353e-08
1.663566e-01|2.917226e-07|4.852999e-08
1.663566e-01|2.351657e-07|3.912135e-08
1.663565e-01|2.800888e-07|4.659460e-08
1.663565e-01|2.556844e-07|4.253477e-08
1.663565e-01|2.101852e-07|3.496566e-08
1.663564e-01|2.143122e-07|3.565221e-08
1.663564e-01|2.513780e-07|4.181833e-08
1.663563e-01|2.311899e-07|3.845991e-08
1.663563e-01|2.342579e-07|3.897027e-08
1.663563e-01|2.178142e-07|3.623475e-08
1.663562e-01|2.358231e-07|3.923065e-08
1.663562e-01|2.402907e-07|3.997384e-08
It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.663561e-01|2.565224e-07|4.267407e-08
1.663561e-01|2.200381e-07|3.660468e-08
1.663561e-01|2.081655e-07|3.462959e-08
1.663560e-01|2.437787e-07|4.055405e-08
1.663560e-01|1.918928e-07|3.192251e-08
1.663560e-01|2.028104e-07|3.373873e-08
1.663559e-01|1.620420e-07|2.695665e-08
1.663559e-01|1.769910e-07|2.944350e-08
1.663559e-01|1.901214e-07|3.162781e-08
1.663558e-01|2.248171e-07|3.739964e-08
1.663558e-01|1.947031e-07|3.238999e-08
1.663558e-01|1.973556e-07|3.283124e-08
1.663557e-01|1.827801e-07|3.040651e-08
1.663557e-01|2.029551e-07|3.376274e-08
1.663557e-01|1.749633e-07|2.910614e-08
1.663556e-01|2.044018e-07|3.400339e-08
1.663556e-01|1.773860e-07|2.950915e-08
1.663556e-01|1.810515e-07|3.011893e-08
1.663556e-01|1.726333e-07|2.871850e-08
1.663555e-01|1.629448e-07|2.710677e-08
It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.663555e-01|1.550954e-07|2.580097e-08
1.663555e-01|2.053401e-07|3.415945e-08
1.663554e-01|1.657553e-07|2.757430e-08
1.663554e-01|1.622977e-07|2.699911e-08
1.663554e-01|1.745864e-07|2.904339e-08
1.663554e-01|1.643929e-07|2.734764e-08
1.663553e-01|1.451168e-07|2.414095e-08
1.663553e-01|1.620844e-07|2.696361e-08
1.663553e-01|1.460663e-07|2.429890e-08
1.663553e-01|1.483641e-07|2.468114e-08
1.663552e-01|1.413107e-07|2.350778e-08
1.663552e-01|1.347709e-07|2.241984e-08
1.663552e-01|1.524089e-07|2.535402e-08
1.663552e-01|1.638024e-07|2.724937e-08
1.663551e-01|1.581713e-07|2.631261e-08
1.663551e-01|1.224170e-07|2.036470e-08
1.663551e-01|1.521732e-07|2.531479e-08
1.663551e-01|1.430895e-07|2.380367e-08
1.663550e-01|1.314782e-07|2.187206e-08
1.663550e-01|1.678995e-07|2.793092e-08
It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.663550e-01|1.240491e-07|2.063619e-08
1.663550e-01|1.419641e-07|2.361643e-08
1.663550e-01|1.385635e-07|2.305073e-08
1.663549e-01|1.442790e-07|2.400153e-08
1.663549e-01|1.388963e-07|2.310608e-08
1.663549e-01|1.226898e-07|2.041005e-08
1.663549e-01|1.168132e-07|1.943245e-08
1.663548e-01|1.193300e-07|1.985112e-08
1.663548e-01|1.228166e-07|2.043113e-08
1.663548e-01|1.267728e-07|2.108926e-08
1.663548e-01|1.377074e-07|2.290828e-08
1.663548e-01|1.246988e-07|2.074424e-08
1.663547e-01|1.178106e-07|1.959835e-08
1.663547e-01|1.290212e-07|2.146328e-08
1.663547e-01|1.210762e-07|2.014159e-08
1.663547e-01|1.305489e-07|2.171742e-08
1.663547e-01|1.080205e-07|1.796971e-08
1.663546e-01|1.077142e-07|1.791876e-08
1.663546e-01|1.231012e-07|2.047846e-08
1.663546e-01|1.154717e-07|1.920926e-08
It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.663546e-01|1.099821e-07|1.829603e-08
1.663546e-01|1.014430e-07|1.687551e-08
1.663545e-01|1.086045e-07|1.806686e-08
1.663545e-01|9.685350e-08|1.611202e-08
1.663545e-01|9.256526e-08|1.539865e-08
1.663545e-01|9.798911e-08|1.630093e-08
1.663545e-01|9.658834e-08|1.606790e-08
1.663545e-01|1.198947e-07|1.994502e-08
1.663544e-01|1.006155e-07|1.673784e-08
1.663544e-01|1.147466e-07|1.908860e-08
1.663544e-01|9.885043e-08|1.644420e-08
1.663544e-01|9.013360e-08|1.499412e-08
1.663544e-01|1.057008e-07|1.758380e-08
1.663544e-01|1.131763e-07|1.882737e-08
1.663543e-01|1.109836e-07|1.846260e-08
1.663543e-01|9.082014e-08|1.510832e-08
1.663543e-01|9.574359e-08|1.592736e-08
1.663543e-01|8.459704e-08|1.407308e-08
1.663543e-01|8.863245e-08|1.474439e-08
1.663543e-01|1.011727e-07|1.683052e-08
It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.663542e-01|9.515093e-08|1.582876e-08

(<Axes: >, <Axes: >, <Axes: >)

Solve EMD with entropic regularization

def f(G):
    return np.sum(G * np.log(G))

def df(G):
    return np.log(G) + 1.0

reg = 1e-3

Ge = ot.optim.cg(a, b, M, reg, f, df, verbose=True)

pl.figure(3, figsize=(5, 5))
ot.plot.plot1D_mat(a, b, Ge, "OT matrix Entrop. reg")

It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.692289e-01|0.000000e+00|0.000000e+00
1.617643e-01|4.614437e-02|7.464513e-03
1.612567e-01|3.147977e-03|5.076323e-04
1.611010e-01|9.664089e-04|1.556895e-04
1.610273e-01|4.578871e-04|7.373232e-05
1.610002e-01|1.683365e-04|2.710220e-05
1.609934e-01|4.185219e-05|6.737929e-06
1.609917e-01|1.064782e-05|1.714210e-06
1.609564e-01|2.195914e-04|3.534463e-05
1.609528e-01|2.227366e-05|3.585008e-06
1.609484e-01|2.731066e-05|4.395607e-06
1.609439e-01|2.812403e-05|4.526390e-06
1.609384e-01|3.411076e-05|5.489731e-06
1.609360e-01|1.478550e-05|2.379520e-06
1.609188e-01|1.071264e-04|1.723866e-05
1.609155e-01|2.010090e-05|3.234547e-06
1.609016e-01|8.640753e-05|1.390311e-05
1.609001e-01|9.787362e-06|1.574787e-06
1.608956e-01|2.802633e-05|4.509312e-06
1.608847e-01|6.752917e-05|1.086441e-05
It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.608791e-01|3.492262e-05|5.618319e-06
1.608754e-01|2.294895e-05|3.691921e-06
1.608696e-01|3.589257e-05|5.774024e-06
1.608676e-01|1.259259e-05|2.025740e-06
1.608633e-01|2.634022e-05|4.237176e-06
1.608621e-01|7.499837e-06|1.206440e-06
1.608570e-01|3.219846e-05|5.179346e-06
1.608565e-01|2.584003e-06|4.156537e-07
1.608519e-01|2.877456e-05|4.628442e-06
1.608495e-01|1.527916e-05|2.457645e-06
1.608491e-01|1.907754e-06|3.068605e-07
1.608458e-01|2.080458e-05|3.346330e-06
1.608400e-01|3.632834e-05|5.843049e-06
1.608332e-01|4.224549e-05|6.794476e-06
1.608273e-01|3.668126e-05|5.899346e-06
1.608229e-01|2.693062e-05|4.331061e-06
1.608189e-01|2.476870e-05|3.983276e-06
1.608189e-01|4.398866e-07|7.074207e-08
1.608168e-01|1.263226e-05|2.031481e-06
1.608136e-01|2.039754e-05|3.280201e-06
It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
1.608108e-01|1.737333e-05|2.793818e-06
1.608081e-01|1.680411e-05|2.702236e-06
1.608064e-01|1.060323e-05|1.705067e-06
1.608040e-01|1.448272e-05|2.328879e-06
1.608016e-01|1.517089e-05|2.439503e-06
1.608016e-01|1.182131e-07|1.900885e-08
1.607998e-01|1.083225e-05|1.741824e-06
1.607960e-01|2.408012e-05|3.871985e-06
1.607935e-01|1.530176e-05|2.460423e-06
1.607934e-01|8.743494e-07|1.405896e-07
1.607930e-01|2.306865e-06|3.709278e-07
1.607916e-01|8.753873e-06|1.407549e-06
1.607916e-01|2.126766e-08|3.419660e-09
1.607913e-01|1.834191e-06|2.949219e-07
1.607895e-01|1.111217e-05|1.786721e-06
1.607895e-01|7.818800e-08|1.257181e-08
1.607866e-01|1.811853e-05|2.913216e-06
1.607864e-01|1.322890e-06|2.127027e-07
1.607864e-01|5.215344e-09|8.385563e-10

(<Axes: >, <Axes: >, <Axes: >)

Solve EMD with Frobenius norm + entropic regularization

def f(G):
    return 0.5 * np.sum(G**2)


def df(G):
    return G


reg1 = 1e-3
reg2 = 1e-1

Gel2 = ot.optim.gcg(a, b, M, reg1, reg2, f, df, verbose=True)

pl.figure(4, figsize=(5, 5))
ot.plot.plot1D_mat(a, b, Gel2, "OT entropic + matrix Frob. reg")
pl.show()

It.  |Loss        |Relative loss|Absolute loss
------------------------------------------------
    0|1.693084e-01|0.000000e+00|0.000000e+00
/home/circleci/project/ot/bregman/_sinkhorn.py:666: UserWarning: Sinkhorn did not converge. You might want to increase the number of iterations `numItermax` or the regularization parameter `reg`.
  warnings.warn(
    1|1.610202e-01|5.147342e-02|8.288260e-03
    2|1.610179e-01|1.406304e-05|2.264402e-06
    3|1.610174e-01|3.352083e-06|5.397436e-07
    4|1.610174e-01|0.000000e+00|0.000000e+00

Comparison of the OT matrices

nvisu = 40

pl.figure(5, figsize=(10, 4))

pl.subplot(2, 2, 1)
pl.imshow(G0[:nvisu, :], cmap="gray_r")
pl.axis("off")
pl.title("Exact OT")

pl.subplot(2, 2, 2)
pl.imshow(Gl2[:nvisu, :], cmap="gray_r")
pl.axis("off")
pl.title("Frobenius reg.")

pl.subplot(2, 2, 3)
pl.imshow(Ge[:nvisu, :], cmap="gray_r")
pl.axis("off")
pl.title("Entropic reg.")

pl.subplot(2, 2, 4)
pl.imshow(Gel2[:nvisu, :], cmap="gray_r")
pl.axis("off")
pl.title("Entropic + Frobenius reg.")

Exact OT, Frobenius reg., Entropic reg., Entropic + Frobenius reg.

Text(0.5, 1.0, 'Entropic + Frobenius reg.')

Total running time of the script: (0 minutes 0.774 seconds)

Gallery generated by Sphinx-Gallery