add log and epsilon scaling stabilizations

Author: rflamary

README.md

@@ -8,7 +8,7 @@ This open source Python library provide several solvers for optimization problem
 It provides the following solvers:
 
 * OT solver for the linear program/ Earth Movers Distance [1].
-* Entropic regularization OT solver  with Sinkhorn Knopp Algorithm [2].
+* Entropic regularization OT solver with Sinkhorn Knopp Algorithm [2] and stabilized version [9][10].
 * Bregman projections for Wasserstein barycenter [3] and unmixing [4].
 * Optimal transport for domain adaptation with group lasso regularization [5]
 * Conditional gradient [6] and Generalized conditional gradient for regularized OT [7].
@@ -98,3 +98,7 @@ This toolbox benefit a lot from open source research and we would like to thank
 [7] Rakotomamonjy, A., Flamary, R., & Courty, N. (2015). Generalized conditional gradient: analysis of convergence and applications. arXiv preprint arXiv:1510.06567.
 
 [8] M. Perrot, N. Courty, R. Flamary, A. Habrard, "Mapping estimation for discrete optimal transport", Neural Information Processing Systems (NIPS), 2016.
+
+[9] Schmitzer, B. (2016). Stabilized Sparse Scaling Algorithms for Entropy Regularized Transport Problems. arXiv preprint arXiv:1610.06519.
+
+[10] Chizat, L., Peyré, G., Schmitzer, B., & Vialard, F. X. (2016). Scaling algorithms for unbalanced transport problems. arXiv preprint arXiv:1607.05816.

examples/demo_OT_1D.py

@@ -19,8 +19,8 @@
 x=np.arange(n,dtype=np.float64)
 
 # Gaussian distributions
-a=gauss(n,m=20,s=20) # m= mean, s= std
-b=gauss(n,m=60,s=60)
+a=gauss(n,m=20,s=5) # m= mean, s= std
+b=gauss(n,m=60,s=10)
 
 # loss matrix
 M=ot.dist(x.reshape((n,1)),x.reshape((n,1)))

examples/demo_optim_OTreg.py

@@ -17,8 +17,8 @@
 x=np.arange(n,dtype=np.float64)
 
 # Gaussian distributions
-a=ot.datasets.get_1D_gauss(n,m=20,s=20) # m= mean, s= std
-b=ot.datasets.get_1D_gauss(n,m=60,s=60)
+a=ot.datasets.get_1D_gauss(n,m=20,s=5) # m= mean, s= std
+b=ot.datasets.get_1D_gauss(n,m=60,s=10)
 
 # loss matrix
 M=ot.dist(x.reshape((n,1)),x.reshape((n,1)))
@@ -37,7 +37,7 @@ 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(3)
@@ -47,9 +47,9 @@ def df(G): return G
 
 def f(G): return np.sum(G*np.log(G))
 def df(G): return np.log(G)+1
-    
+
 reg=1e-3
-  
+
 Ge=ot.optim.cg(a,b,M,reg,f,df,verbose=True)
 
 pl.figure(4)