jeanfeydy / geomloss Goto Github PK

Very great package! My question is how to set the gpu device?

Hi there. Does GeomLoss support batching? I want to calculate the Sinkhorn distance between multiple pairs all at once, even better If I can do it in a one-to-many fashion.

I try this simple setup:

x = torch.randn(2, 50, 300)
y = torch.randn(2, 20, 300)

sinkhorn = SamplesLoss(loss='sinkhorn', p=2, blur=.05, scaling=.9)

sinkhorn(x, y)

and I get this error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-50-72739cb485d1> in <module>
      4 sinkhorn = SamplesLoss(loss='sinkhorn', p=2, blur=.05, scaling=.9, backend='tensorized', cost=lambda x, y: torch.cdist(x, y))
      5 
----> 6 sinkhorn(x, y)

~/.conda/envs/dlab/lib/python3.7/site-packages/torch/nn/modules/module.py in __call__(self, *input, **kwargs)
    530             result = self._slow_forward(*input, **kwargs)
    531         else:
--> 532             result = self.forward(*input, **kwargs)
    533         for hook in self._forward_hooks.values():
    534             hook_result = hook(self, input, result)

~/.conda/envs/dlab/lib/python3.7/site-packages/geomloss/samples_loss.py in forward(self, *args)
    235                     debias = self.debias, potentials = self.potentials,
    236                     labels_x = l_x, labels_y = l_y,
--> 237                     verbose = self.verbose )
    238 
    239 

~/.conda/envs/dlab/lib/python3.7/site-packages/geomloss/sinkhorn_samples.py in sinkhorn_tensorized(α, x, β, y, p, blur, reach, diameter, scaling, cost, debias, potentials, **kwargs)
     50     a_x, b_y, a_y, b_x = sinkhorn_loop( softmin_tensorized, 
     51                                         log_weights(α), log_weights(β),
---> 52                                         C_xx, C_yy, C_xy, C_yx, ε_s, ρ, debias=debias )
     53 
     54     return sinkhorn_cost(ε, ρ, α, β, a_x, b_y, a_y, b_x, batch=True, debias=debias, potentials=potentials)

~/.conda/envs/dlab/lib/python3.7/site-packages/geomloss/sinkhorn_divergence.py in sinkhorn_loop(softmin, α_logs, β_logs, C_xxs, C_yys, C_xys, C_yxs, ε_s, ρ, jumps, kernel_truncation, truncate, cost, extrapolate, debias, last_extrapolation)
    160         # "Coordinate ascent" on the dual problems:
    161         if debias:
--> 162             at_x = λ * softmin(ε, C_xx, α_log + a_x/ε )  # OT(α,α)
    163             bt_y = λ * softmin(ε, C_yy, β_log + b_y/ε )  # OT(β,β)
    164         at_y = λ * softmin(ε, C_yx, α_log + b_x/ε )  # OT(α,β) wrt. a

RuntimeError: The size of tensor a (50) must match the size of tensor b (100) at non-singleton dimension 1

I even tried this initialization:

sinkhorn = SamplesLoss(loss='sinkhorn', p=2, blur=.05, scaling=.9, backend='tensorized', cost=lambda x, y: torch.cdist(x, y))

with no luck.

Thanks in advance.

First of all, geomloss is a great llibrary. It's quite ez to use in my project. Thx for ur effort for such a great library.

Recently, I have tried SampleLoss for my data. Everything is going well at the begining. After few epochs training, it raised a ValueError:

`
File "/home/lly/prj/ASTM/20news/Gamma_T50-parameter1/model.py", line 65, in Sinkhorn_divergence

return self.OTD(alpha, beta)

File "/home/lly/anaconda3/lib/python3.8/site-packages/torch/nn/modules/module.py", line 727, in _call_impl

result = self.forward(*input, **kwargs)

File "/home/lly/anaconda3/lib/python3.8/site-packages/geomloss/samples_loss.py", line 231, in forward

values = routines[self.loss][backend]( α, x, β, y, 

File "/home/lly/anaconda3/lib/python3.8/site-packages/geomloss/sinkhorn_samples.py", line 48, in sinkhorn_tensorized

diameter, ε, ε_s, ρ = scaling_parameters( x, y, p, blur, reach, diameter, scaling )

File "/home/lly/anaconda3/lib/python3.8/site-packages/geomloss/sinkhorn_divergence.py", line 72, in scaling_parameters

ε_s = epsilon_schedule( p, diameter, blur, scaling )

File "/home/lly/anaconda3/lib/python3.8/site-packages/geomloss/sinkhorn_divergence.py", line 60, in epsilon_schedule

+ [ np.exp(e) for e in np.arange(p*np.log(diameter), p*np.log(blur), p*np.log(scaling)) ] \

ValueError: arange: cannot compute length
`

Then I switch to some vanilla distribution samples: LogNormal distribution for alpha and LogNormal with 0 mean and 1.0 std for beta. After few epochs training, this error came again. So I wonder what should I do to keep it working properly. Thank you!

Hi! In the website it is said: "full support of PyTorch’s autograd engine".

However, in the function sinkhorn_loop there is a line in the beginnig:

torch.autograd.set_grad_enabled(False)

And the gradients are turned on again only in the end.

So, as I see it, I cannot differentiate through the solution of the Optimal Transport problem (i.e. through the "flows between suppliers and demanders"). Is it correct?

The reason why I am asking is the presence of the library qpth which enables differentiation through the solution. See this example. And my goal is to pass (alpha, x, beta, y) as input where all four tensors require gradients (they are "Variables" in terms of older versions of PyTorch). So I am trying to understand if should just copy your implementation and comment the line above, or may be I am missing something...

Hello Jean,
the library looks really nice, thank you for your work!

Unfortunately, I don't manage to run the code sample in the index page. I get the following error:

line 73, in keops_lse
    log_conv = generic_logsumexp("( B - (P * " + cost + " ) )",
NameError: name 'generic_logsumexp' is not defined

It happens every time I try to use the multiscale Sinkhorn. However, tensorized versions of Sinkhorn seem to work well, both on CPU and GPU. Has this issue been referenced before?

Thanks

Hello! Thanks for your work again!

I want to see the computation graph on tensorboradx but it seems impossible because there are many "if" statements. So I was wondering that had you visualized the computation graph when computing Wasserstein loss?

Dear @jeanfeydy ,

Thanks for your great. When will be the next update? e.g., out-of-the-box ImageLoss and VolumeLoss layers

Looking forward to updating.

Kindest regards,
Jun

I am trying to computer Sinkhorn distance between two tensor of size 2x256. My data are quite sparse.
I got the folowing error:

File "/usr/local/lib/python3.6/dist-packages/geomloss/sinkhorn_samples.py", line 48, in sinkhorn_tensorized diameter, ε, ε_s, ρ = scaling_parameters( x, y, p, blur, reach, diameter, scaling ) File "/usr/local/lib/python3.6/dist-packages/geomloss/sinkhorn_divergence.py", line 72, in scaling_parameters ε_s = epsilon_schedule( p, diameter, blur, scaling ) File "/usr/local/lib/python3.6/dist-packages/geomloss/sinkhorn_divergence.py", line 61, in epsilon_schedule + [ blur**p ] ValueError: Maximum allowed size exceeded
Why this is happening and how can i solve this?
Any help is highly appreciated. Thanks

Hi,

I'm very excited about this work! I was wondering if ImagesLoss with Sinkhorn divergence would be usable (when it's ready) as a reconstruction loss for an image autoencoder for representation learning?

Pseudocode to code to illustrate what I was thinking about:

encoder = ConvNet()
decoder = DeconvNet()
optimizer = Optimizer(chain(encoder.parameters(), decoder.parameters()))
loss_fn = ImagesLoss(loss="sinkhorn", p=2, blur=.05)

for images in dataset:
  optimizer.zero_grad()
  encoding = encoder(images)
  image_preds = decoder(encoding)
  loss = loss_fn(images, image_preds)
  loss.backward()
  optimizer.step()

Thanks,
Agost

Hello,

Thank you for this wonderful project!

I am wondering on something, hope someone can help.

Consider a toy example with a small size.

a = torch.randn(10,8,5)

I am looking for the most efficient way to create a distance matrix D of size (10,5) such that :

D[i,j] = OTLoss(a[i,:,0], a[i,:,j])

Evidently D[i,0] = 0 for all the is
Here i runs from 0 to 9 and j runs from 0 to 4, thus D has the shape of (10,5).

Now assuming that a is a very large 3D-tensor (let's say (100k,100,100)), what is the most efficient way to get this D tensor?
I am hoping that there exist some solution without using loops.

Hello,

And thank you for the awesome library!

I was wondering if there is a way to evaluate the potential on a point that it wasn't trained on. For instance, adapting your example from here:

import torch
from geomloss import SamplesLoss

x = torch.randn(100, 3, requires_grad=True).cuda()
y = torch.randn(100, 3).cuda()
OT_solver = SamplesLoss(loss = "sinkhorn", p = 2, blur = 0.05,
                        debias = False, potentials = True)
F, G = OT_solver(x, y)  # Dual potentials

I would like e.g. to find the value of F in torch.tensor([0., 0., 0.]). I can certainly interpolate between points of x, but is there a more direct way to evaluate the network? I could not find an example of it.

Thank you so much!

Hi Jean,

Thank you very much for you amazing work. Good luck with your PhD dissertation!

When I tried to run your sample "plot_optimal_transport_2D.py" with GPU, there is an error that I could not figure out:

Error message:

Traceback (most recent call last):
File "plot_optimal_transport_2D.py", line 149, in <module> gradient_descent( SamplesLoss("sinkhorn", p=2, blur=.1) )
File "plot_optimal_transport_2D.py", line 107, in gradient_descent L_αβ = loss(x_i, y_j)
File "/home/me/anaconda3/lib/python3.7/site-packages/torch/nn/modules/module.py", line 541, in __call__
result = self.forward(*input, **kwargs)
File "/home/me/anaconda3/lib/python3.7/site-packages/geomloss/samples_loss.py", line 237, in forward verbose = self.verbose )
File "/home/me/anaconda3/lib/python3.7/site-packages/geomloss/sinkhorn_samples.py", line 102, in sinkhorn_online C_xx, C_yy, C_xy, C_yx, ε_s, ρ, debias=debias )
File "/home/me/anaconda3/lib/python3.7/site-packages/geomloss/sinkhorn_divergence.py", line 151, in sinkhorn_loop a_x = λ * softmin(ε, C_xx, α_log ) # OT(α,α)
File "/home/me/anaconda3/lib/python3.7/site-packages/geomloss/sinkhorn_samples.py", line 69, in softmin_online
return - ε * log_conv( x, y, f_y.view(-1,1), torch.Tensor([1/ε]).type_as(x) ).view(-1)
File "/home/me/anaconda3/lib/python3.7/site-packages/pykeops/torch/generic/generic_red.py", line 351, in __call__ out = GenredAutograd.apply(self.formula, self.aliases, backend, self.dtype, device_id, ranges, *args)
File "/home/me/anaconda3/lib/python3.7/site-packages/pykeops/torch/generic/generic_red.py", line 43, in forward *args)
RuntimeError: [KeOps] This KeOps shared object has been compiled without cuda support: try to set tagHostDevice to 0 or recompile the formula with a working version of cuda.

To test whether I have installed pykeops correctly, I tested the code:

import torch
import pykeops.torch as pktorch

x = torch.arange(1, 10, dtype=torch.float32).view(-1, 3)
y = torch.arange(3, 9, dtype=torch.float32).view(-1, 3)

my_conv = pktorch.Genred('SqNorm2(x-y)', ['x = Vi(3)', 'y = Vj(3)'])
print(my_conv(x, y))

and I received:

Error message

/usr/include/crt/host_config.h:121:2: error: #error -- unsupported GNU version! gcc versions later than 6 are not supported!
#error -- unsupported GNU version! gcc versions later than 6 are not supported!
^~~~~
CMake Error at keopslibKeOpstorch91c92bd508_generated_link_autodiff.cu.o.Release.cmake:219 (message):
Error generating
/home/velysianp/.cache/pykeops-1.2-cpython-37/build-
libKeOpstorch91c92bd508/CMakeFiles/keopslibKeOpstorch91c92bd508.dir/keops/core/./keopslibKeOpstorch91c92bd508_generated_link_autodiff.cu.o

make[3]: *** [CMakeFiles/keopslibKeOpstorch91c92bd508.dir/keops/core/keopslibKeOpstorch91c92bd508_generated_link_autodiff.cu.o] Error 1
make[2]: *** [CMakeFiles/keopslibKeOpstorch91c92bd508.dir/all] Error 2
make[1]: *** [CMakeFiles/libKeOpstorch91c92bd508.dir/rule] Error 2
make: *** [libKeOpstorch91c92bd508] Error 2

Currently I am using Ubuntu 18.04.3, with Python3.7 downloaded from Anaconda. Also, I my GPU is Geforce 1080 and I downloaded CUDA 10.2 with its driver and toolkit. The only thing I am worrying about is the gcc version, since now I use "gcc --version" and it returned gcc (Ubuntu 7.4.0-1ubuntu1~18.04.1) 7.4.0. I noticed that CUDA documentation mentioned that the gcc version of Ubuntu 18.04.3 with x86_64 architecture should be 7.3.0. I am wondering if I am having problem with incompatible gcc version with my system that resulted in the problem? Or do you have any other ideas or suggestions?

Thank you very much.
Elyson

Hello,

I'm trying to run a simple toy problem, where I have a target Multivariate Normal and a tunable mean, std Multivariate to try to match the target. The loss would be the sinkhorn loss between the samples from the target and the samples from the tunable model and the mean and std should be optimized by something like the ADAM optimizer. But it does not seem like the gradients are being back-propagated to the mean/std parameters to modify their values (when I check values using torch.autograd the gradients come out as None). I'm not sure what I'm doing wrong. Any help would be appreciated! Thanks!

import torch
import torch.nn as nn
import torch.optim as optim

from torch.distributions.multivariate_normal import MultivariateNormal

from geomloss import SamplesLoss

DIM_SIZE = 3
NUM_SAMPLES = 64
LR = 0.01
NUM_STEPS = 2000

target_mean = torch.ones(DIM_SIZE)
target_std = torch.ones(DIM_SIZE)*0.25
target_m = MultivariateNormal(target_mean, torch.diag(target_std))

mean = nn.Parameter(torch.zeros(DIM_SIZE, requires_grad = True))
std = nn.Parameter(torch.ones(DIM_SIZE, requires_grad = True))
optimizer = optim.Adam([mean, std], lr=LR)

loss = SamplesLoss(loss="sinkhorn", p=2, blur=.05)

for i in range(NUM_STEPS):
	m = MultivariateNormal(mean, torch.diag(std))

	target_samples = target_m.sample(sample_shape=torch.Size([NUM_SAMPLES]))
	param_samples = m.sample(sample_shape=torch.Size([NUM_SAMPLES])).requires_grad_()

	sinkhorn_loss = loss(param_samples, target_samples)
	# sinkhorn_loss = sinkhorn_loss.sum()
	optimizer.zero_grad()
	sinkhorn_loss.backward()
	optimizer.step()

	# g_mean, g_std = torch.autograd.grad(sinkhorn_loss, [mean, std])
	# print (g_mean, g_std)
	
	print('Error on mean: {} std: {} Sinkhorn Loss: {}'.format(torch.mean(torch.abs(target_mean-mean.cpu())), torch.mean(torch.abs(target_std-std.cpu())), sinkhorn_loss))

Hello,
Thanks a lot for your great librairy !
We are trying to compute the transportation cost with Geomloss and we obtain inconsitent results when compared with ot. We want to use non uniforme weight with a custom cost function. Maybe you could help us to spot the error ?

OT

import numpy as np 
import ot.gpu as otgpu
size = 10000
weight = np.random.rand(size)
D_1 = weight / np.linalg.norm(weight, 1)
M = np.random.rand(size,size)
X =otgpu.sinkhorn(D_1,D_2,M_cupy, reg=0.05)
transportation_cost = np.sum(np.multiply(X,M))

We have transportation_cost=0.05004567861659549

Geomloss

def cost_function_custom(X,Y): 
    return torch.tensor(M).unsqueeze(0).cuda()

w_x = torch.tensor(D_1).cuda().clone() 
w_y = torch.tensor(D_2).cuda().clone() 
loss = SamplesLoss(loss="sinkhorn", p=2, blur=.05,backend="tensorized",cost=cost_function_custom,debias =False,scaling = 0.9999)
loss = loss(w_x,w_x.unsqueeze(1),w_y,w_y.unsqueeze(1))

We obtain loss=tensor(0.0150, device='cuda:0', dtype=torch.float64)

Are we missing something ?

Thanks a lot for your time.
Cheers

when I try to compute the loss of a batch data, I met this bug.

import torch
from geomloss import SamplesLoss  # See also ImagesLoss, VolumesLoss
cuda_device = torch.device("cuda:%d" % 0 if torch.cuda.is_available() else "cpu")
x=torch.randn(100,90,400, requires_grad=True).to(cuda_device)
y=torch.randn(100,90,400).to(cuda_device)
# Define a Sinkhorn (~Wasserstein) loss between sampled measures
loss = SamplesLoss(loss="sinkhorn", p=2, blur=.05)
L = loss(x, y)  # By default, use constant weights = 1/number of samples
print(L.item())
g_x, = torch.autograd.grad(L, [x])  # GeomLoss fully supports autograd!
print(g_x)

output

Traceback (most recent call last):
  File "/home/lowen/program/pycharm-community-2018.2.3/helpers/pydev/pydevd.py", line 1664, in <module>
    main()
  File "/home/lowen/program/pycharm-community-2018.2.3/helpers/pydev/pydevd.py", line 1658, in main
    globals = debugger.run(setup['file'], None, None, is_module)
  File "/home/lowen/program/pycharm-community-2018.2.3/helpers/pydev/pydevd.py", line 1068, in run
    pydev_imports.execfile(file, globals, locals)  # execute the script
  File "/home/lowen/program/pycharm-community-2018.2.3/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile
    exec(compile(contents+"\n", file, 'exec'), glob, loc)
  File "/devdata/new_Relation_Extraction/test_geomloss.py", line 40, in <module>
    L = loss(x, y)  # By default, use constant weights = 1/number of samples
  File "/home/lowen/anaconda3/envs/pytorch/lib/python3.7/site-packages/torch/nn/modules/module.py", line 547, in __call__
    result = self.forward(*input, **kwargs)
  File "/home/lowen/anaconda3/envs/pytorch/lib/python3.7/site-packages/geomloss/samples_loss.py", line 239, in forward
    verbose=self.verbose)
  File "/home/lowen/anaconda3/envs/pytorch/lib/python3.7/site-packages/geomloss/sinkhorn_samples.py", line 52, in sinkhorn_tensorized
    C_xx, C_yy, C_xy, C_yx, ε_s, ρ, debias=debias )
  File "/home/lowen/anaconda3/envs/pytorch/lib/python3.7/site-packages/geomloss/sinkhorn_divergence.py", line 162, in sinkhorn_loop
    at_x = λ * softmin(ε, C_xx, α_log + a_x/ε )  # OT(α,α)
RuntimeError: The size of tensor a (90) must match the size of tensor b (9000) at non-singleton dimension 1

So I debug the code and found that a_x([9000] tensors) dismatched α_log ([100,90] tensors), then I try to apply a view() operation to α_log just like this

at_x = λ * softmin(ε, C_xx, α_log + (a_x / ε).view(α_log.size()))

After I fixed all this kind of bugs with a view() operation, it successfully returns a batch loss(in my case, that's [100] tensors).

So I was wondering is it the right way to fix the bug? or there is a batter way to compute wasserstein loss for a batch data?

Thx!!!

Hello Jean !

Thanks for your amazing work ! I wonder how I should cite your geomloss framework. Do you have a paper ? Or will you intend to have one soon ?

Thanks again,

Best,

Kilian

Hi Jean,

Thanks for the library! I hope your PhD thesis writing's going fine.

I seem to have found an insidious bug. When we initialize measures as 2-D tensors with shapes like (N, 1) and (M, 1), check_shapes() tries to reshape them to be 1-D like (N,) and (M,) but these changes are not propagated to the original variables in the function calling check_shapes(). This leads to the following RunTimeError for me:

-> loss_p2 = Loss_p2(a_i, x_i, b_j, y_j)
(Pdb) c
Traceback (most recent call last):
  File "/home/rajanr/anaconda3/envs/py36/lib/python3.6/pdb.py", line 1667, in main
    pdb._runscript(mainpyfile)
  File "/home/rajanr/anaconda3/envs/py36/lib/python3.6/pdb.py", line 1548, in _runscript
    self.run(statement)
  File "/home/rajanr/anaconda3/envs/py36/lib/python3.6/bdb.py", line 434, in run
    exec(cmd, globals, locals)
  File "<string>", line 1, in <module>
  File "/home/rajanr/custom-gym-env/Wasserstein_sinkhorn_stuff.py", line 86, in <module>
    scaling=scaling, backend="multiscale")
  File "/home/rajanr/anaconda3/envs/py36/lib/python3.6/site-packages/torch/nn/modules/module.py", line 541, in __call__
    result = self.forward(*input, **kwargs)
  File "/home/rajanr/anaconda3/envs/py36/lib/python3.6/site-packages/geomloss/samples_loss.py", line 237, in forward
    verbose = self.verbose )
  File "/home/rajanr/anaconda3/envs/py36/lib/python3.6/site-packages/geomloss/sinkhorn_samples.py", line 262, in sinkhorn_multiscale
    debias = debias)
  File "/home/rajanr/anaconda3/envs/py36/lib/python3.6/site-packages/geomloss/sinkhorn_divergence.py", line 162, in sinkhorn_loop
    at_x = λ * softmin(ε, C_xx, α_log + a_x/ε )  # OT(α,α)
  File "/home/rajanr/anaconda3/envs/py36/lib/python3.6/site-packages/geomloss/sinkhorn_samples.py", line 117, in softmin_multiscale
    return - ε * log_conv( x, y, f_y.view(-1,1), torch.Tensor([1/ε]).type_as(x), ranges=ranges_xy ).view(-1)
  File "/home/rajanr/anaconda3/envs/py36/lib/python3.6/site-packages/pykeops/torch/generic/generic_red.py", line 351, in __call__
    out = GenredAutograd.apply(self.formula, self.aliases, backend, self.dtype, device_id, ranges, *args)
  File "/home/rajanr/anaconda3/envs/py36/lib/python3.6/site-packages/pykeops/torch/generic/generic_red.py", line 43, in forward
    *args)
RuntimeError: [KeOps] Wrong value of the 'j' dimension 0for arg number 2 : is 9 but was 3 in previous 'j' arguments.
Uncaught exception. Entering post mortem debugging
Running 'cont' or 'step' will restart the program

The code is here:

Would you like me to fix it and create a pull request?

Greetings,
Raghu.

do you have such example?
Thanks.

name 'generic_logsumexp' is not defined

a = torch.rand(28996, 2)
b = torch.rand(28996, 2)

samplesloss = SamplesLoss('sinkhorn')
samplesloss(a, b)

Hi, I'm the author of PyTorch Geometry and I would like to discuss with you potential collaborations between packages. We already have a defined losses module which we would like to expand with your stuff. We have a consolidated testing framework and supporting cuda pretty soon in our contiguous integration pipeline. I couldn't find any email contact. Please, try to reach me out.

Bonjour Jean,

First of all, thanks a lot for sharing this wonderful work with the community. Great job!

I am trying to use OT for mesh registration. I basically use the code of your tutorial here: https://www.kernel-operations.io/geomloss/_auto_examples/optimal_transport/plot_interpolation_3D.html

I just start from my own template mesh and keep track of the mesh triangulation. Note that the target and template are not closed meshes.

Here is the template (black) and target meshes before registration:

The final match is given here:

And the taget here:

As you can see the result are not good between the thumb and the 4 fingers. Would you be able to tell me if I (probably) did something wrong or if this kind of failure could be expected with OT? In the first case, any idea of what is wrong and how to solve it? (I know it quite an open question! :))

I tried to add a mesh constraint in the objective function (just add a L2 norm which penalizes changes in the edge lengths), but it does not solve this "split" problem.

I also tried to use landmarks by adding a template/target L2 mismatch term in the OF, but it did not work well. Is there a better way to use landmarks with the sinkhorn loss?

Thanks a lot in advance and do not hesitate to answer "Sorry dude but it could be anything and I do not have time and enough info to figure out what is wrong" :)

Merci,
Théophile

Hi, I met error "The size of tensor a (2048) must match the size of tensor b (4096) at non-singleton dimension 1" when I run the following code:

points1 = torch.rand(2, 2048, 3, requires_grad=True).cuda()
points2 = torch.rand(2, 2048, 3).cuda()
geomloss.SamplesLoss(loss='sinkhorn', p=2, blur=0.05)(points1 , points2)

It seems that the code does not process one batch input well. I wonder how to solve it. Thanks!

I am trying to use the sklearn.metrics.pairwise.haversine_distances as the cost for SamplesLoss, in order to calculate the EMD for spherical distributions (see [https://arxiv.org/abs/2012.11116](EMLight: Lighting Estimation via Spherical Distribution)).

I am feeding in the sample locations in spherical coordinates (angles - dimensions [M,2]). The problem is when computing C_xx, C_yy, C_xy, C_yx, the cost function cannot do this for batches. What is the best way to efficiently implement this?

First of all thanks for the great library!

I just tried to run SampleLoss with batches of data and it did not work.

So, I have two tensors x,y of the same shape [batch_dim, n_points, feature_dim] and wish to compute the sinkhorn divergence between the point clouds x[0] and y[0], x[1] and y[1] in a batched way (to prevent slow loops) in order to return a tensor of shape [batch_dim].

However, when trying this out with SampleLoss() I receive a shape error.

To reproduce a minimal example I add the following collab here:
https://colab.research.google.com/drive/1NqagWVIv-a8YN258NcFEBXbRBFAVuuiR?usp=sharing

Test code

import torch
from geomloss import SamplesLoss

A = torch.randn((5, 100, 1))
B = torch.randn((5, 100, 1))

L_ab = SamplesLoss('sinkhorn', p=1, blur=0.1, reach=None, debias=True, potentials=True, backend='tensorized')

f, g = L_ab(A, B)

C = A[0]
D = B[0]

h, i = L_ab(C, D)

a1 = f[0] - torch.mean(f[0])
a2 = h - torch.mean(h)

err = torch.norm(a1 - a2) / torch.norm(a2)

print(err)

Currently an error is thrown.

RuntimeError: The size of tensor a (100) must match the size of tensor b (500) at non-singleton dimension 1

After debugging a little bit.
The problem seems to be in sinkhorn_samples.py line 33, where the batch is not treated properly.
The problme is resolved when replacing line 33 by

return - ε * ( f.view(B,1,-1) - C/ε ).logsumexp(2).view(B, -1)

And the output now are like

tensor(0.0048)

which is reasonable to me.

What do you think?

Hi, great work! I love the visualizations and comparisons of MMD, Sinkhorn, and OT.

When I run the code sample here, I'm getting a NameError. Any tips for how to resolve it?

can you suggest a way we can use cosine similarity loss for the cost function in sinkhorn divergence?
currently i can only see option for l1 and l2 norm.

Hi,

First of all, thanks for the great library!

I have a question. The web site says that the divergences are symmetric and positive definite. Can we have any claims about if/when triangular inequality may hold for the divergences implemented here?

I am looking into using the library in a problem which needs a metric measure, but I couldn't find any information if there is support for this.

Thanks!

Essentially whenever I calculate the sinkhorn loss, pytorch's internal gradient tracking turns back on.

I believe this is to do with line 280 in geomloss/sinkhorn_divergence.py:
torch.autograd.set_grad_enabled(True)

This is obviously not an issue during training, but it is causing my memory to explode during inference or validation as the gradients are not cleared as there typically is no backwards pass.

I think this could be fixed by using something like:
cur_grad_status = torch.is_grad_enabled() ## At beginning of sinkhorn_loop()
torch.autograd.set_grad_enabled(cur_grad_status) ## At end of sinkhorn_loop()

First of all, great work!

At the moment, I am running a test scenario, where precision matters. Hence, I want to use double tensors (I am aware that this is not efficient on most GPUs). I am trying to compute the energy distance (but the problem is the same with other loss types), which works fine as long as the backend is tensorized. Scaling this up to larger measures is mainly a memory issue and not really a computation speed issue.
If online backend is used instead, I get an error, see the attached error log (independent of CPU or GPU use). The problem seems to be somewhere in the interplay between pykeops and PyTorch. Is this behavior intended ? I am running PyTorch 1.3, pykeops 1.2 and geomloss 0.2.3, all installed using conda and pip (I also checked some older combinations without any succes). If necessary, I can try to build a minimal running example.

Thanks,
Sebastian
error_log_geomloss.txt

Edit: I found the issue: Inside the kernel_samples.py file, there is the datatype missing when calling generic_sum. Hence, the default value float32 is used and types are not matching.

Hi,
First of all, thank you for making this library publicly available.

As indicated in the title, I have encountered examples where the loss from SamplesLoss returns a negative value.
I have saved such an example and it consists of two sets of 50 points living in R^160. They are quite close as the maximum distance between the points of those two sets is about 0.25 (see notebook below).

The loss is the default sinkhorn divergence and I have tried both the default parameters and more accurate ones.

default_loss = geomloss.SamplesLoss()
accurate_loss = geomloss.SamplesLoss(blur=1e-3, scaling=0.999)

The results are the following:
Default Loss: -0.00108
More Accurate Loss: -0.000130

I have reproduced such behavior on my computer and on Google Colab. Here is the notebook.

It is possible that it is just me who misconfigured something. Please let me know if this is the case.
Otherwise, if you think that it shouldn't happen with those parameters, I can look at the code more thoroughly to try to fix it.

Thank you,
Bests,
Waïss

Hi, I'm using this code for density estimation and generative modelling. I think Sinkhorn divergence is perfect for such tasks because of it unbiased nature. However, in practice, I find it is slight biased.

Basically, I have several datapoints and a Gaussian distribution. I want to match them so the datapoints can be regarded as samples from the Gaussian distribution. This is a basic setting in generative modelling. Ideally, this will work well because Sinkhorn divergence gets rid of entropic bias of regularized W-distance by adding two self correlation terms. However, the results are still biased.

I present my simple test code here.

import os
import time

import numpy as np
import torch

from geomloss import SamplesLoss


eps = 0.1
sample = 2000
lr = 0.001
epoch = 400
p = 1

# Synthesis data points.
device = 'cuda'
tdtype = torch.float

rng = np.random.RandomState(0)
Train_all = rng.randn(1000, 2)
train_tensor = torch.tensor(Train_all, device=device, dtype=tdtype).requires_grad_(True)

d=Train_all.shape[1]


# Gaussian model
class Gaussian():
    def __init__(self, Sample_n, u, sigma, eps=0.5, p=1, tdtype=torch.float, device='cuda'):



        self.u = torch.tensor(u, requires_grad=True, device=device, dtype=tdtype)
        self.sigma = torch.tensor(sigma, requires_grad=True, device=device, dtype=tdtype)


        self.L_ab = SamplesLoss('energy', potentials=False, backend='tensorized', scaling=0.9)
        # self.L_ab = SamplesLoss('sinkhorn', p=p, blur=eps, debias=True, potentials=False, backend='tensorized', scaling=0.9)
        # self.L_ab = SamplesLoss('sinkhorn', p=p, blur=eps, debias=False, potentials=False, backend='tensorized', scaling=0.9)


        self.dim = self.u.shape[-1]

        self.n1 = Sample_n 
        self.eps = eps
        self.tdtype = tdtype
        self.device = device


    def logp_x(self, x):
        out = - self.dim / 2 * np.log(2 * np.pi) - torch.sum(torch.log(self.sigma), -1) - 0.5 * torch.sum((torch.unsqueeze(x, 0) - self.u)** 2 / (self.sigma ** 2), -1)
        return out       

    def Sample(self):
        # with torch.no_grad():
        Sample = torch.randn((1, self.n1, self.dim), device=self.device, dtype=self.tdtype)
        tmp = Sample * self.sigma + self.u
        out = torch.reshape(tmp, [-1, self.dim])
        return out


    def loss(self, D):

        x = self.Sample()
        out = self.L_ab(x, D)
        return out


# Init the Gaussian
alpha = np.ones(1) 
u_ = np.zeros((1, 1, d))
S_ = np.ones((1, 1, d))
model = Gaussian(sample, u_, S_, eps=eps, p=1)


a_opt = torch.optim.RMSprop([
                            {'params': model.u},
                            {'params': model.sigma}
                            ] , lr=lr, alpha=0.9)

p_opt = torch.optim.RMSprop([
                            {'params': train_tensor},
                            ] , lr=lr, alpha=0.9)


for epoch in range(epoch):

    lp = model.loss(train_tensor)

    # Uncomment to optimize datapoint
    # p_opt.zero_grad()
    # lp.backward()
    # p_opt.step()
    # print('Data mean:\t', train_tensor.detach().mean(0).cpu().numpy(), '\t Data std:\t', train_tensor.detach().std(0).cpu().numpy())
    # print('Gaussian mean:\t', model.u[:,0,:].detach().cpu().numpy(), '\t \t \t Gaussian std:\t', model.sigma.detach()[:,0,:].cpu().numpy())

    # Uncomment to optimize Gaussian dsitribution
    a_opt.zero_grad()
    lp.backward()
    a_opt.step()
    print('Data mean:\t', train_tensor.detach().mean(0).cpu().numpy(), '\t Data std:\t', train_tensor.detach().std(0).cpu().numpy())
    print('Gaussian mean:\t', model.u[:,0,:].detach().cpu().numpy(), '\t Gaussian std:\t', model.sigma.detach()[:,0,:].cpu().numpy())

I present several experimental facts of this code:

• Sinkhorn divergence suffers from bias, and this bias decrease as \epsilon increase, in the limit, MMD works perfectly.
  When I fix the datapoint and learn the Gaussian distribution. The std of data is 0.97487277. With eps=0.01, 0.1, 0.5, 1, the learned stds of Gaussian are 0.96189475, 0.9575431, 0.96259254, 0.97405905, which are all smaller than the data std. With MMD (uncomment line 38), the learned std are near data std, it could be larger or smaller than than data std in different runs.

• This is also true if I learn the data while fixed the Gaussian.
  Uncomment line 92-96.
  The std of Gaussian distribution is 1. With eps=0.01, 0.1, 0.5, 1, the data stds are 0.98536015, 0.9891638, 0.9949019, 0.9966444. With MMD, the learned data std is 1.0002198.

Excuse my poor understanding towards EMD. I test geomloss with the following code:

points1 = torch.rand(2, 1024, 2, requires_grad=True).cuda()
points2 = torch.rand(2, 1024, 2).cuda()
emd_sinkhorn = geomloss.SamplesLoss(loss='sinkhorn', p=2, blur=0.01, backend='auto')(points1, points2)

Also, I try to use another EMD library (https://github.com/Colin97/MSN-Point-Cloud-Completion/tree/master/emd) which uses auction algorithm:

emd-auction, assignment = emd(points1, points2, eps=0.05, iters=2000)

However, two experiments give different results. How can I change the usage of geomloss to match another one? Thanks for your help very much!

Hi,

First, thanks a lot for a great package that is very straightforward to use!

I am confused by some results I am getting where the Sinkhorn distance (with default parameters) is negative. From your website I understand that this shouldn't be the case.

Minimal example:

import numpy as np
import torch
from geomloss import SamplesLoss

x0 = torch.from_numpy(np.load("x0.npy"))
x1 = torch.from_numpy(np.load("x1.npy"))

sinkhorn = SamplesLoss("sinkhorn")

print(sinkhorn(x0, x1).item())

with the data files here:
mwe.zip

With geomloss v0.2.3 on the CPU this prints -0.022376153618097305.

Am I misunderstanding the Sinkhorn divergence, is there something anomalous in my data, or is this some numerical artifact?

Thanks a lot,
Johann

Hi,

I've tried following the links in the Matchs and algorithms section of the docs, and I am getting an Apache server error 404 for the author's website at the ENS.

For example, the link to the gradients flow tutorial.

Thanks,

Hi there,

Love the package. One thing I noticed is that I wanted to try using a few custom cost functions. One I used: C(x,y) = 0.5*|| x - y ||^2 only gets the same as the inbuilt p=2 for some sample points, but can occasionally disagree.

Any thoughts on this?

Many thanks,
Jake

Hi!
I am using the suggested !pip install geomloss[full] on google colab. When I try to run the color transfer example I get:

Compiling libKeOpstorch05a7ba1328 in /root/.cache/pykeops-1.4.2-cpython-37:
       formula: Max_SumShiftExp_Reduction(( B - (P * (SqDist(X,Y) / IntCst(2)) ) ),0)
       aliases: X = Vi(0,3); Y = Vj(1,3); B = Vj(2,1); P = Pm(3,1); 
       dtype  : float32
... 
--------------------- MAKE DEBUG -----------------
Command '['cmake', '--build', '.', '--target', 'libKeOpstorch05a7ba1328', '--', 'VERBOSE=1']' returned non-zero exit status 2.
/usr/local/lib/python2.7/dist-packages/cmake/data/bin/cmake -H/usr/local/lib/python3.7/dist-packages/pykeops -B/root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328 --check-build-system CMakeFiles/Makefile.cmake 0
/usr/bin/make -f CMakeFiles/Makefile2 libKeOpstorch05a7ba1328
make[1]: Entering directory '/root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328'
/usr/local/lib/python2.7/dist-packages/cmake/data/bin/cmake -H/usr/local/lib/python3.7/dist-packages/pykeops -B/root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328 --check-build-system CMakeFiles/Makefile.cmake 0
/usr/local/lib/python2.7/dist-packages/cmake/data/bin/cmake -E cmake_progress_start /root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328/CMakeFiles 5
/usr/bin/make -f CMakeFiles/Makefile2 CMakeFiles/libKeOpstorch05a7ba1328.dir/all
make[2]: Entering directory '/root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328'
/usr/bin/make -f CMakeFiles/keopslibKeOpstorch05a7ba1328.dir/build.make CMakeFiles/keopslibKeOpstorch05a7ba1328.dir/depend
make[3]: Entering directory '/root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328'
cd /root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328 && /usr/local/lib/python2.7/dist-packages/cmake/data/bin/cmake -E cmake_depends "Unix Makefiles" /usr/local/lib/python3.7/dist-packages/pykeops /usr/local/lib/python3.7/dist-packages/pykeops /root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328 /root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328 /root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328/CMakeFiles/keopslibKeOpstorch05a7ba1328.dir/DependInfo.cmake --color=
make[3]: Leaving directory '/root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328'
/usr/bin/make -f CMakeFiles/keopslibKeOpstorch05a7ba1328.dir/build.make CMakeFiles/keopslibKeOpstorch05a7ba1328.dir/build
make[3]: Entering directory '/root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328'
[ 20%] Building CUDA object CMakeFiles/keopslibKeOpstorch05a7ba1328.dir/keops/core/link_autodiff.cu.o
/usr/local/cuda/bin/nvcc -ccbin=/usr/bin/c++ -DCUDA_BLOCK_SIZE=192 -DC_CONTIGUOUS=1 -DENABLECHUNK=1 -DKERNEL_GEOM_TYPE=0 -DKERNEL_SIG_TYPE=0 -DKERNEL_SPHERE_TYPE=0 -DMAXIDGPU=0 -DMAXTHREADSPERBLOCK0=1024 -DMODULE_NAME=libKeOpstorch05a7ba1328 -DMODULE_NAME_FSHAPE_SCP=fshape_scp_gaussiangaussiangaussian_unoriented_float -DSHAREDMEMPERBLOCK0=49152 -DSUM_SCHEME=1 -DUSE_CUDA=1 -DUSE_DOUBLE=0 -DUSE_HALF=0 -D_FORCE_INLINES -D_GLIBCXX_USE_CXX11_ABI=0 -D__TYPEACC__=float -D__TYPE__=float -DkeopslibKeOpstorch05a7ba1328_EXPORTS -I/usr/local/cuda/targets/x86_64-linux/include -I/usr/local/lib/python3.7/dist-packages/pykeops -I/usr/local/lib/python3.7/dist-packages/pykeops/keops -I/root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328 -I/usr/local/lib/python3.7/dist-packages/torch/include -I/usr/local/lib/python3.7/dist-packages/torch/include/torch/csrc/api/include  -gencode arch=compute_70,code=sm_70 --use_fast_math --compiler-options=-fPIC -O3 -DNDEBUG -Xcompiler=-fPIC   -include libKeOpstorch05a7ba1328.h -std=c++14 -x cu -c /usr/local/lib/python3.7/dist-packages/pykeops/keops/core/link_autodiff.cu -o CMakeFiles/keopslibKeOpstorch05a7ba1328.dir/keops/core/link_autodiff.cu.o
CMakeFiles/keopslibKeOpstorch05a7ba1328.dir/build.make:62: recipe for target 'CMakeFiles/keopslibKeOpstorch05a7ba1328.dir/keops/core/link_autodiff.cu.o' failed
make[3]: Leaving directory '/root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328'
CMakeFiles/Makefile2:331: recipe for target 'CMakeFiles/keopslibKeOpstorch05a7ba1328.dir/all' failed
make[2]: Leaving directory '/root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328'
CMakeFiles/Makefile2:306: recipe for target 'CMakeFiles/libKeOpstorch05a7ba1328.dir/rule' failed
make[1]: Leaving directory '/root/.cache/pykeops-1.4.2-cpython-37/build-libKeOpstorch05a7ba1328'
Makefile:196: recipe for target 'libKeOpstorch05a7ba1328' failed

Would you mind taking a look?

When I tried to run the python script

import torch

from geomloss import SamplesLoss
WD = SamplesLoss(loss='sinkhorn', p=2, blur=.05)
a = torch.randn((8,4096,3))
b = torch.randn((8,4096,3))
c= WD(a,b)
print(c)

with the latest release (0.2.3), an error occurs

Traceback (most recent call last):
File "t.py", line 7, in
c= WD(a,b)
File "/home/miniconda3/envs/py37/lib/python3.7/site-packages/torch/nn/modules/module.py", line 493, in call
result = self.forward(*input, **kwargs)
File "/home/miniconda3/envs/py37/lib/python3.7/site-packages/geomloss/samples_loss.py", line 237, in forward
verbose = self.verbose )
File "/home/miniconda3/envs/py37/lib/python3.7/site-packages/geomloss/sinkhorn_samples.py", line 52, in sinkhorn_tensorized
C_xx, C_yy, C_xy, C_yx, ε_s, ρ, debias=debias )
File "/home/miniconda3/envs/py37/lib/python3.7/site-packages/geomloss/sinkhorn_divergence.py", line 162, in sinkhorn_loop
at_x = λ * softmin(ε, C_xx, α_log + a_x/ε ) # OT(α,α)
RuntimeError: The size of tensor a (4096) must match the size of tensor b (32768) at non-singleton dimension 1

but the same script works well with an older version (0.2.1) in the same environment.

Hello Jean, thank you for this incredible library. I read somewhere that you planned on releasing several things among which barycenter support but that you still had to go from prototype to the feature. Is there a way to have access to the prototype ? Or would some help be useful in releasing this feature or another one ?

Thanks a lot !

Hi Jean,

Thank you for building up this awesome library! I would like to do MMD loss for comparing images, do you know how to do it? I saw there are ImageLoss and VolumeLoss in the documentation, do you plan to implement them? If not what is the workaround?

Best,
Daiqing

Hi Jean,

Very nice work congratulations!

I am quite new to optimal transport but I was wondering why you only implemented 2 scale Sinkhorn algorithm for the multi-scale approach? It seems to me that just adding a for loop with the number of loop being computed on the fly based on input parameters or distribution characteristic would not change much and could have some benefits. Maybe it is actually more complex to propagate the approximate solution over multiple scales?

Thank you in advance for your help.

Best,

Samuel

Hi,

This line detaches the gradient for the weights associated with the locations for (I believe) no real reason as it serves the same purpose as the cost matrix in the fixed point iteration.

It's probably worth detaching a_x, a_y, b_x and b_y only.

Hello, I wrote a small test where I calculate the Sinkhorn distance between a target distribution (univariate normal with mean 0 and variance 1) and a bunch of sample distributions (univariate normals with different means and variances), and the results I get are quite surprising.

If my sampled distribution is a normal with mean 0 and variance < 1, the Sinkhorn distance from the target distribution is lower than if my sampled distribution has variance 1. This seems to indicate that the loss is biased, as I would expect Sinkhorn to be at a minimum when the two distributions have same mean and variance. I played a bit with the parameters of SamplesLoss but the result didn't change significantly.

I attached a script below that does a mean and variance sweep and saves a bunch of plots with the corresponding Sinkhorn distance values. I'd be grateful if you could shed some light on this.

import matplotlib.pyplot as plt
import numpy as np
import torch

from geomloss import SamplesLoss


np.random.seed(31)
torch.manual_seed(31)
torch.backends.cudnn.deterministic = True
torch.backends.cudnn.benchmark = False


def normal(imgshape, mu, sigma):
    """
    Sample a random tensor from a normal
    """  
    return torch.randn(imgshape) * sigma + mu 


def main():

    N = 10000

    mus = np.arange(-5, 5, 0.5)
    sigmas = np.arange(0.1, 2, 0.1)

    target = normal((1, N), mu=0, sigma=1)

    musweep = np.zeros(len(mus))
    sigmasweep = np.zeros(len(sigmas))

    sinkhorn_loss = SamplesLoss(loss="sinkhorn", p=1, blur=.05)
    
    # Cycle thru all mus
    for i, mu in enumerate(mus):
        sample = normal((1, N), mu=mu, sigma=1)

        sinkhorn_D = sinkhorn_loss(sample, target)
        musweep[i] = sinkhorn_D

        print("mu:", mu, "sinkhorn:", sinkhorn_D)

        # Log histograms
        fig = plt.figure(figsize=(5, 4), dpi=500)
        fig.tight_layout(pad=0)
        ax = fig.add_subplot(111)
        plt.subplots_adjust(wspace=0, hspace=0)
        ax.margins(0)
        plt.hist(sample.flatten().cpu().numpy(), bins=np.arange(-4, 4 + 0.05, 0.05), alpha=0.5)
        plt.hist(target.flatten().cpu().numpy(), bins=np.arange(-4, 4 + 0.05, 0.05), alpha=0.5)
        plt.title(f'sample vs target')
        plt.legend(["sample", "target"])
        plt.xlim(-4, 4)
        fig.text(.5, 0, 
            f"sinkhorn={float(sinkhorn_D):.4}", ha="center")
        plt.savefig(f"mu_{mu}.png")
        plt.close()
        
    # Cycle thru all sigmas
    for i, sigma in enumerate(sigmas):
        sample = normal((1, N), mu=0, sigma=sigma)

        sinkhorn_D = sinkhorn_loss(sample, target)
        sigmasweep[i] = sinkhorn_D

        print("sigma:", sigma, "sinkhorn:", sinkhorn_D)

        # Log histograms
        fig = plt.figure(figsize=(5, 4), dpi=500)
        fig.tight_layout(pad=0)
        ax = fig.add_subplot(111)
        plt.subplots_adjust(wspace=0, hspace=0)
        ax.margins(0)
        plt.hist(sample.flatten().cpu().numpy(), bins=np.arange(-4, 4 + 0.05, 0.05), alpha=0.5)
        plt.hist(target.flatten().cpu().numpy(), bins=np.arange(-4, 4 + 0.05, 0.05), alpha=0.5)
        plt.title(f'sample vs target')
        plt.legend(["sample", "target"])
        plt.xlim(-4, 4)
        fig.text(.5, 0, 
            f"sinkhorn={float(sinkhorn_D):.4}", ha="center")
        plt.savefig(f"sigma_{sigma:.2}.png")
        plt.close()

        im = plt.imread(f"sigma_{sigma:.2}.png")
    
    # Plot sinkhorn vs mu
    fig = plt.figure(figsize=(5, 4), dpi=500)
    fig.tight_layout(pad=0)
    ax = fig.add_subplot(111)
    plt.subplots_adjust(wspace=0, hspace=0)
    ax.margins(0)
    plt.plot(mus, musweep)
    plt.title(f'sinkhorn vs mu')
    plt.xlim(-4, 4)
    plt.savefig(f"sinkhorn_vs_mu.png")
    plt.close()

    # Plot sinkhorn vs sigma
    fig = plt.figure(figsize=(5, 4), dpi=500)
    fig.tight_layout(pad=0)
    ax = fig.add_subplot(111)
    plt.subplots_adjust(wspace=0, hspace=0)
    ax.margins(0)
    plt.plot(sigmas, sigmasweep)
    plt.title(f'sinkhorn vs sigma')
    plt.xlim(-4, 4)
    plt.savefig(f"sinkhorn_vs_sigma.png")
    plt.close()


if __name__ == "__main__":
    main()

The following code causes the crash - a degenerate case of loss between Dirac distribution and itself:

import torch
from geomloss import SamplesLoss

loss_func = SamplesLoss("sinkhorn")

t1_samples = torch.tensor([[1.0, ], ])
t2_samples = torch.tensor([[1.0, ], ])

loss_value = loss_func(t1_samples, t2_samples)

Hello,

When I was working on calculating the Sinkhorn loss, I am getting the below issue

truth is a matrix of shape (100,100)
I would appreciate the help! Thanks!

With regards,
Vishal

Hello!
How to get the transport plan for wasserstein distance? For example:

    import torch
    from geomloss import SamplesLoss  # See also ImagesLoss, VolumesLoss
    from tensorboardX import SummaryWriter
    # Create some large point clouds in 3D
    x = torch.randn(100, 90, requires_grad=True).cuda()
    y = torch.randn(100, 90).cuda()
    x_weight =torch.softmax(torch.randn(100, requires_grad=True).cuda(),dim=-1)
    y_weight = torch.softmax(torch.randn(100, requires_grad=True).cuda(),dim=-1)
    # Define a Sinkhorn (~Wasserstein) loss between sampled measures
    loss = SamplesLoss(loss="sinkhorn", p=2, blur=.05)
    L = loss(x_weight, x,y_weight, y)  # By default, use constant weights = 1/number of samples
    g_x_weight, = torch.autograd.grad(L, [x_weight])  # GeomLoss fully supports autograd!
    print(g_x_weight)

output:
    tensor([38.3446, 40.1507, 54.9848, 24.2403, 36.6134, 21.7442, 40.9704, 37.5197,
        36.8129, 30.8433, 41.8184, 37.8540, 29.5628, 32.1728, 35.4814, 34.3499,
        30.5784, 29.4510, 30.2563, 35.3756, 34.6008, 21.2113, 49.6258, 33.0905,
        26.7779, 39.9584, 27.6789, 33.0485, 35.1336, 31.8450, 33.2152, 49.8619,
        27.9988, 31.8357, 32.6381, 36.5462, 22.4289, 34.1696, 23.5495, 37.9015,
        37.1510, 30.8962, 42.6643, 41.9360, 32.8569, 52.4191, 48.5090, 32.2841,
        39.2593, 32.8378, 24.9272, 29.9081, 27.0655, 25.6891, 36.1147, 38.3208,
        43.5676, 40.3591, 32.9711, 22.7926, 33.5953, 26.1191, 36.2213, 26.2793,
        29.8464, 37.0326, 43.1056, 30.1160, 42.6618, 31.8065, 36.9675, 24.9423,
        47.4061, 27.5967, 22.2735, 25.6721, 36.1677, 33.0408, 28.5409, 32.0098,
        31.1463, 33.6331, 28.8633, 20.3415, 32.0224, 32.7113, 29.5103, 49.1259,
        38.1893, 29.8794, 41.5620, 31.8520, 30.0871, 29.5324, 17.6336, 38.2577,
        38.8874, 28.6118, 36.6356, 31.2856], device='cuda:0')

And I want to check the meanning of g_x_weight. I mean, when I backward the loss and update params, what's the change in the model (asume x and x_weight are outputs of the model)? After update, will both x and x_weight be close to y and y_weight next time?

Hi all, I have successfully installed geomloss. However, when I run the script, the following Bug results. Any idea about it? thank you

command

python transfer_labels.py

log:

Fetching the atlas... Done.
Fetching the atlas... Done.
Fetching the atlas... Done.
Data loaded.
[pyKeOps] Compiling libKeOpstorch5fbe6e35e9 in /home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0:
      formula: ArgMin_Reduction(Sum(Square((Var(0,60,0) - Var(1,60,1)))),0)
      aliases: Var(0,60,0); Var(1,60,1);
      dtype  : float32
...
/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/cmake_scripts/script_keops_formula/../../keops/core/formulas/maths/TensorDot.h:201:267: error: expansion pattern ‘tao::seq::map_t<typename tao::seq::difference<typename tao::seq::impl::memoize_sequence<long unsigned int, tao::seq::impl::sequence_size<I>::value>::type, CONTFA>::type, DIMFA>’ contains no argument packs
/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/cmake_scripts/script_keops_formula/../../keops/core/formulas/maths/TensorDot.h:201:410: error: expansion pattern ‘tao::seq::map_t<typename tao::seq::difference<typename tao::seq::impl::memoize_sequence<long unsigned int, tao::seq::impl::sequence_size<DIMFB>::value>::type, CONTFB>::type, DIMFB>’ contains no argument packs
compilation terminated due to -fmax-errors=2.
make[3]: *** [CMakeFiles/copy_KeOps_formula.dir/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/keops/core/link_autodiff.cu.o] Error 1
make[2]: *** [CMakeFiles/copy_KeOps_formula.dir/all] Error 2
make[1]: *** [CMakeFiles/KeOps_formula.dir/rule] Error 2
make: *** [KeOps_formula] Error 2

--------------------- MAKE DEBUG -----------------
Command '['cmake', '--build', '.', '--target', 'KeOps_formula', '--', 'VERBOSE=1']' returned non-zero exit status 2.
/usr/bin/cmake -S/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/cmake_scripts/script_keops_formula -B/home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82 --check-build-system CMakeFiles/Makefile.cmake 0
/usr/bin/make  -f CMakeFiles/Makefile2 KeOps_formula
make[1]: Entering directory '/home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82'
/usr/bin/cmake -S/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/cmake_scripts/script_keops_formula -B/home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82 --check-build-system CMakeFiles/Makefile.cmake 0
/usr/bin/cmake -E cmake_progress_start /home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82/CMakeFiles 1
/usr/bin/make  -f CMakeFiles/Makefile2 CMakeFiles/KeOps_formula.dir/all
make[2]: Entering directory '/home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82'
/usr/bin/make  -f CMakeFiles/copy_KeOps_formula.dir/build.make CMakeFiles/copy_KeOps_formula.dir/depend
make[3]: Entering directory '/home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82'
cd /home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82 && /usr/bin/cmake -E cmake_depends "Unix Makefiles" /home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/cmake_scripts/script_keops_formula /home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/cmake_scripts/script_keops_formula /home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82 /home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82 /home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82/CMakeFiles/copy_KeOps_formula.dir/DependInfo.cmake --color=
make[3]: Leaving directory '/home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82'
/usr/bin/make  -f CMakeFiles/copy_KeOps_formula.dir/build.make CMakeFiles/copy_KeOps_formula.dir/build
make[3]: Entering directory '/home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82'
[100%] Building CUDA object CMakeFiles/copy_KeOps_formula.dir/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/keops/core/link_autodiff.cu.o
/usr/local/cuda/bin/nvcc  -ccbin=/usr/bin/c++ -DCUDA_BLOCK_SIZE=192 -DC_CONTIGUOUS=1 -DMAXIDGPU=0 -DMAXTHREADSPERBLOCK0=1024 -DSHAREDMEMPERBLOCK0=49152 -DUSE_CUDA=1 -DUSE_DOUBLE=0 -DUSE_HALF=0 -D_FORCE_INLINES -D_GLIBCXX_USE_CXX11_ABI=0 -D__TYPE__=float -I/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/cmake_scripts/script_keops_formula/../.. -I/home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82 -I/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/cmake_scripts/script_keops_formula/../../keops -I/usr/local/cuda/targets/x86_64-linux/include -gencode arch=compute_61,code=sm_61 --use_fast_math --compiler-options=-fPIC -DUSE_OPENMP -Xcompiler -fopenmp,-Wall,-Wno-unknown-pragmas,-fmax-errors=2 -O3 -DNDEBUG -Xcompiler=-fPIC -include KeOps_formula.h -std=c++14 -x cu -c /home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/keops/core/link_autodiff.cu -o CMakeFiles/copy_KeOps_formula.dir/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/keops/core/link_autodiff.cu.o
CMakeFiles/copy_KeOps_formula.dir/build.make:76: recipe for target 'CMakeFiles/copy_KeOps_formula.dir/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/keops/core/link_autodiff.cu.o' failed
make[3]: Leaving directory '/home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82'
CMakeFiles/Makefile2:109: recipe for target 'CMakeFiles/copy_KeOps_formula.dir/all' failed
make[2]: Leaving directory '/home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82'
CMakeFiles/Makefile2:90: recipe for target 'CMakeFiles/KeOps_formula.dir/rule' failed
make[1]: Leaving directory '/home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0/build-7e6973de82'
Makefile:124: recipe for target 'KeOps_formula' failed

--------------------- ----------- -----------------
[pyKeOps] Compiling pybind11 template libKeOps_template_354d16605e in /home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0 ... Traceback (most recent call last):
 File "transfer_labels.py", line 253, in <module>
   C_i_flat, labs_i = KMeans(X_i_flat, C_i_flat)
 File "transfer_labels.py", line 223, in KMeans
   labs_i = nn_search(x_i, c_j, ranges=ranges)
 File "transfer_labels.py", line 210, in nn_search
   return D_ij.argmin(dim=1).view(-1)
 File "/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/common/lazy_tensor.py", line 1929, in argmin
   return self.reduction("ArgMin", axis=axis, **kwargs)
 File "/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/common/lazy_tensor.py", line 744, in reduction
   return res()
 File "/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/common/lazy_tensor.py", line 929, in __call__
   return self.callfun(*args, *self.variables, **self.kwargs)
 File "/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/torch/generic/generic_red.py", line 579, in __call__
   *args
 File "/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/torch/generic/generic_red.py", line 48, in forward
   formula, aliases, dtype, "torch", optional_flags, include_dirs
 File "/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/common/keops_io.py", line 48, in __init__
   self._safe_compile()
 File "/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/common/utils.py", line 75, in wrapper_filelock
   func_res = func(*args, **kwargs)
 File "/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/common/keops_io.py", line 63, in _safe_compile
   self.build_folder,
 File "/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/common/compile_routines.py", line 246, in compile_generic_routine
   dtype, lang, include_dirs, use_prebuilt_formula=True
 File "/home/gamorosino/miniconda3/lib/python3.7/site-packages/pykeops/common/compile_routines.py", line 67, in get_or_build_pybind11_template
   os.mkdir(template_build_folder)
FileExistsError: [Errno 17] File exists: '/home/gamorosino/.cache/pykeops-1.5-cpython-37-gpu0//build-pybind11_template-libKeOps_template_354d16605e'

I am aware that Geomloss supports batching, however my distributions have different number of samples, and the only way I could think of creating the batches was to pad the tensors. But, obviously the results are incorrect this way.

Any suggestions to overcome this problem?

Self contained example:

import torch
import torch.nn.functional as F
from geomloss import SamplesLoss

N = 20

# I want to compute the loss between d and every s_i.
# Instead of iterating: [loss(d, s1), loss(d, s2), loss(d, s3)],
# I want to compute it using batching: loss([d, d, d], [s1, s2, s3])
# The problem is that the number of samples are different between s_i.
# I pad the tensors, however, the results are obviously wrong.

# Original tensors
d = torch.rand(size=(80, 100))

s1 = torch.rand(size=(10, 100))
s2 = torch.rand(size=(15, 100))
s3 = torch.rand(size=(20, 100))

# Padded tensors
d_ = d.unsqueeze(0).repeat(3,1,1)

s_ = torch.stack([
    F.pad(s1, (0, 0, 0, N - s1.shape[0])),
    F.pad(s2, (0, 0, 0, N - s2.shape[0])),
    F.pad(s3, (0, 0, 0, N - s3.shape[0]))
])

# Sinkhorn
loss = SamplesLoss()

assert(loss(d_, s_)[0] == loss(d, s1))

Hi,

Wouldn't it be better to have the reach param depend on the diameter of the data? It's probably easier to think of saturation of the distance in terms of %diameter than it is in the direct unit of the distance. What do you think?

Issue

Given the program sinkhorn_loss.py:

from geomloss import SamplesLoss
import torch

loss= SamplesLoss(loss="sinkhorn", p=2, blur=.1)

x = torch.linspace(0,1)[:,None]
y = torch.linspace(0,1)[:,None]
l = loss.forward(x,y)

I would expect an exception to either be told to send my tensors to GPU (if that is required), or for the program to attempt to automatically send the tensors to the correct device.

However, what i get is:

$ python sinkhorn_loss.py
$ Illegal instruction (core dumped)

Running the same code, but sending x,y to the GPU first produces no errors and gives the correct result (loss=0).

Suggested Solution

Perform a check in SamplesLoss to see that the input tensors are on the correct device, and either send an exception or transfer them.