Hi, and thanks for maintaining this great library!

I'm currently using POT (with PyTorch backend) to compute OT-based losses.
I noticed that ot.emd2 seems not differentiable in the usual PyTorch sense,
but when I inspect the .grad field of the cost matrix M, I still get non-zero gradients after calling loss.backward().

So I would like to confirm:
Is ot.emd2 actually differentiable and supports autograd backward propagation through M or input distributions a, b?
Or are the observed gradients just residual tensors from detached operations (i.e., not truly backpropagated through the OT plan solver)?

Minimal reproducible example

import torch
import ot

device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')

n = 10
a = torch.rand(n, device=device)
a = a / a.sum()
b = torch.rand(n, device=device)
b = b / b.sum()

M = torch.randn(n, n, device=device, requires_grad=True)

loss_emd2 = ot.emd2(a, b, M)
loss_emd2.backward()
grad_emd2 = M.grad.clone()
print("EMD2 loss:", loss_emd2.item())
print("EMD2 grad:\n", grad_emd2)

M.grad.zero_()
reg = 0.1
loss_sinkhorn = ot.sinkhorn2(a, b, M, reg)
loss_sinkhorn.backward()
grad_sinkhorn = M.grad.clone()
print("\nSinkhorn loss:", loss_sinkhorn.item())
print("Sinkhorn grad:\n", grad_sinkhorn)

OUTPUT:
EMD2 loss: -1.5605539083480835
EMD2 grad:
tensor([[0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
0.0427],
[0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0095, 0.1873, 0.0095,
0.0000],
[0.0000, 0.0000, 0.0000, 0.0000, 0.0252, 0.0000, 0.0233, 0.0000, 0.0000,
0.0000],
[0.0000, 0.0000, 0.0000, 0.0925, 0.0000, 0.0000, 0.0000, 0.0000, 0.0061,
0.0987],
[0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0481, 0.0000, 0.0000,
0.0000],
[0.0000, 0.0000, 0.0552, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
0.0000],
[0.0000, 0.0126, 0.0000, 0.0000, 0.0000, 0.0000, 0.0532, 0.0000, 0.0000,
0.0000],
[0.0000, 0.0000, 0.0765, 0.0000, 0.0000, 0.0548, 0.0000, 0.0000, 0.0000,
0.0000],
[0.1010, 0.0000, 0.0175, 0.0000, 0.0020, 0.0000, 0.0000, 0.0000, 0.0000,
0.0000],
[0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0844, 0.0000, 0.0000, 0.0000,
0.0000]], device='cuda:0')

Sinkhorn loss: -1.5583062171936035
Sinkhorn grad:
tensor([[-1.0759e-15, -4.5089e-05, -2.3012e-06, -4.7454e-05, -3.1261e-15,
-8.1661e-06, -8.5109e-06, -2.6828e-12, -5.4352e-05, 4.2903e-02],
[-1.3629e-10, -5.7191e-13, -5.8734e-08, -1.4337e-03, -2.9875e-06,
-5.2210e-04, 1.2639e-02, 1.8735e-01, 8.2132e-03, -2.8584e-06],
[-1.3665e-12, -5.0650e-08, -4.7080e-04, -9.4670e-06, 2.7954e-02,
-4.7037e-09, 2.1273e-02, -9.7631e-12, -2.6613e-04, -1.6820e-12],
[-7.8764e-13, -1.0603e-12, -1.3319e-03, 9.4009e-02, -1.1775e-03,
-4.5318e-04, -6.8939e-07, -1.8201e-14, 7.6663e-03, 9.8547e-02],
[-6.4536e-11, -5.2054e-08, -3.9438e-09, -5.9460e-09, -1.3147e-07,
-1.8057e-07, 4.8088e-02, -3.1984e-11, -1.9308e-09, -2.1640e-10],
[-4.6134e-17, -7.1095e-07, 5.5260e-02, -3.3341e-17, -4.1714e-20,
-5.8153e-07, -2.3840e-05, -2.2943e-13, -7.5624e-06, -8.0583e-10],
[-1.0100e-13, 1.3838e-02, -1.4756e-08, -1.8981e-10, -2.9864e-12,
-8.6449e-05, 5.2120e-02, -1.4227e-16, -3.6979e-08, -2.9439e-11],
[-2.2352e-07, -1.1892e-05, 7.6544e-02, -2.5317e-08, -1.2087e-06,
5.4697e-02, -9.0246e-08, -8.5752e-18, -3.6054e-07, -1.7776e-11],
[ 1.0100e-01, -1.0979e-05, 1.9160e-02, -6.2655e-06, 4.4209e-04,
-1.1139e-07, -2.8346e-07, -3.8414e-05, -3.5571e-07, -3.9045e-05],
[-1.1940e-24, -1.1416e-03, -6.1396e-10, -7.2483e-15, -9.0970e-12,
8.5502e-02, -9.0134e-09, -6.0799e-18, -1.8524e-11, -1.1368e-18]],
device='cuda:0')

Please clarify whether:

ot.emd2 is intentionally non-differentiable (since it solves a linear program);
or if there is a plan to add a differentiable variant (e.g. differentiable EMD via implicit function theorem or entropic relaxation);
or if the gradients observed are accidental numerical artifacts.
Thanks a lot for your time and for maintaining this library!