broadinstitute / wot

A software package for analyzing snapshots of developmental processes
https://broadinstitute.github.io/wot/
BSD 3-Clause "New" or "Revised" License
139 stars 34 forks source link

duality_gap overflow #73

Closed MikeDMorgan closed 4 years ago

MikeDMorgan commented 4 years ago

Hi, I'm attempting to debug an error in the optimal transport construction step, run on the command line.

The traceback is:

Computing transport map from 0.5 to 0.75
/homes/mdmorgan/.local/lib/python3.7/site-packages/wot/ot/optimal_transport.py:145: RuntimeWarning: divide by zero encountered in true_divide
  a = (p / (K.dot(np.multiply(b, dy)))) ** alpha1 * np.exp(-u / (lambda1 + epsilon_i))
Traceback (most recent call last):
  File "/homes/mdmorgan/.local/bin/wot", line 10, in <module>
    sys.exit(main())
  File "/homes/mdmorgan/.local/lib/python3.7/site-packages/wot/__main__.py", line 23, in main
    cmd.main(args)
  File "/homes/mdmorgan/.local/lib/python3.7/site-packages/wot/commands/optimal_transport.py", line 30, in main
    tmap_out=args.out)
  File "/homes/mdmorgan/.local/lib/python3.7/site-packages/wot/ot/ot_model.py", line 188, in compute_all_transport_maps
    tmap = self.compute_transport_map(*day_pair)
  File "/homes/mdmorgan/.local/lib/python3.7/site-packages/wot/ot/ot_model.py", line 231, in compute_transport_map
    return self.compute_single_transport_map(config)
  File "/homes/mdmorgan/.local/lib/python3.7/site-packages/wot/ot/ot_model.py", line 304, in compute_single_transport_map
    tmap, learned_growth = wot.ot.compute_transport_matrix(solver=self.solver, **config)
  File "/homes/mdmorgan/.local/lib/python3.7/site-packages/wot/ot/optimal_transport.py", line 30, in compute_transport_matrix
    tmap = solver(**params)
  File "/homes/mdmorgan/.local/lib/python3.7/site-packages/wot/ot/optimal_transport.py", line 186, in optimal_transport_duality_gap
    raise RuntimeError("Overflow encountered in duality gap computation, please report this incident")
RuntimeError: Overflow encountered in duality gap computation, please report this incident

This error arises when computing the transport map between multiple different pairs of time points, so it is not specific to these two.

As far as I can tell the issue arises during the update of a and b variables, because the matrix K contains many small values very close to 0, which subsequently become Inf upon log transformation.

Any thoughts on a possible solution?

MikeDMorgan commented 4 years ago

As an update - this sparsity issue seems to be driven by include a filter gene list - using all genes circumvents the problem. However, I would expect that this makes the reduced dimensional representation much noisier.

NB: Selecting a subset of 100 random genes also works - I assume this is due to the selection of HVGs that are lowly expressed, which exacerbates the sparsity problem. I'll consider this issue resolved for now.

joshua-gould commented 4 years ago

You can also try changing the solver using the command line arguments: --solver fixed_iters

MikeDMorgan commented 4 years ago

Thanks Joshua - I have tried this, but there is still a division by zero, which suggested a sparsity issue; I wasn't confident that even running with very many iterations would get close to convergence.

geoffschieb commented 4 years ago

Hi Mike,

I think this might be due to a small value for the entropy parameter epsilon. Can you try with a bigger value for epsilon?

Best, Geoff

On Wed, Nov 13, 2019 at 5:57 AM MikeDMorgan notifications@github.com wrote:

Thanks Joshua - I have tried this, but there is still a division by zero, which suggested a sparsity issue; I wasn't confident that even running with very many iterations would get close to convergence.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/broadinstitute/wot/issues/73?email_source=notifications&email_token=ACJCQYU7OM5UUQBBDCNBR6DQTQBUBA5CNFSM4JMYJQYKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOED6HCRA#issuecomment-553414980, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACJCQYSOKSASLPBWWRWGXY3QTQBUBANCNFSM4JMYJQYA .

MikeDMorgan commented 4 years ago

Hi Geoff, I tried a range of epsilon values from the default to ~1 with no appreciable effect other than it took several more iterations for the cost matrix to fill with NAs. I suspect a very large epsilon would be required, on the order >> 10 to completely fix this, but then I'm not sure what the consequences would be for interpreting the transport map and downstream analyses. If you have any advice in that respect it would be greatly appreciated.

geoffschieb commented 4 years ago

What values are you using for lambda 1 and 2?

As epsilon goes to infinity, you would get the independent coupling. This would be uninformative — any cell could go to any cell.

Just for the purpose of debugging, it might be nice to see if the code runs for epsilon 1,2,3,...10.

How many cells do you have in both time points?

On Fri, Nov 15, 2019 at 12:26 AM MikeDMorgan notifications@github.com wrote:

Hi Geoff, I tried a range of epsilon values from the default to ~1 with no appreciable effect other than it took several more iterations for the cost matrix to fill with NAs. I suspect a very large epsilon would be required, on the order >> 10 to completely fix this, but then I'm not sure what the consequences would be for interpreting the transport map and downstream analyses. If you have any advice in that respect it would be greatly appreciated.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/broadinstitute/wot/issues/73?email_source=notifications&email_token=ACJCQYQO7HWJJWBFK3LDL5TQTZMMJA5CNFSM4JMYJQYKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEEEV6OQ#issuecomment-554262330, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACJCQYVU2PKGSGDKEMEZ56TQTZMMJANCNFSM4JMYJQYA .