Describe the bug

In optim.gcg, the inner solver is ot.sinkhorn. Thus, if I am not mistaken, I would expect the output of

ot.optim.gcg(a, b, M, reg, 0., lambda x: 0., lambda x: np.zeros(x.shape), numInnerItermax=1000)

and

ot.sinkhorn(a, b, M, reg, numItermax=1000)

to be the same. However, it is not always the case, even when using the same tolerances.

Do I miss something?

To Reproduce

n = 100
X = ot.datasets.make_2D_samples_gauss(n, m=[0, 0], sigma=np.eye(2))
Y = ot.datasets.make_2D_samples_gauss(n, m=[1, 1], sigma=np.eye(2))
M = ot.dist(X, Y)
M /= M.max()
a = ot.unif(n)
b = ot.unif(n)
reg = 1e-2
pi1 = ot.optim.gcg(a, b, M, reg, 0., lambda x: 0., lambda x: np.zeros(x.shape), numInnerItermax=1000)
pi2 = ot.sinkhorn(a, b, M, reg, numItermax=1000)
assert np.allclose(pi1, pi2, atol=1e-5)  # Fail