mikgroup/sigpy

Sensitivity map estimates are not reproducible

ad12 opened this issue · 5 comments

ad12 commented

Describe the bug
Setting the seeds for random number generators in relevant libraries (cupy, numpy, random, etc.) should guarantee that the sensitivity map estimates are identical when using the espirit and jsense calibration apps; but this is not the case.

To Reproduce
Colab: https://colab.research.google.com/drive/113Xkjkuadl476thWy8nlKMjb5VHzBD96?usp=sharing

Expected behavior
Setting the seed for the random number generators for all relevant libraries should guarantee that the maps are the same.

Reproduce
This code is reproducible on colab using the most recent SigPy version (0.1.22) and cupy 7.4.0. Note colab does not support higher versions of cupy just yet, but when this was run locally using sigpy v0.1.22 and cupy 9.0.0, the same issues were seen.

I did briefly look at your code and notice a couple of things which might point to this error.

  • The code uses np.all( ) but you should probably be using np.allclose.
  • The error is quite small, relative to the signal

My first thought is that this may be due to a parallelization issue. With floating point addition/sumation, there can be differences in the floating point error depending on the order of sums (e.g. (A+B)+C != A + ( B + C) ).

ad12 commented

Great points - adding some thoughts below

  1. we tried with np.allclose earlier, and we get the same issue even with higher relative and absolute tolerances (1e-4). Given that complex128 can represent floating point numbers with precision of 1e-15, I would think that we wouldnt see this with a 1e-4 tolerance.
  2. On another note, the espirit calibration would have these floating point operations, in which case we would expect to see some degree of random behavior; however, we consistently get reproducible results with espirit. This is also the case for a few other datasets we have tested locally.
  3. The reproducibility issues with Jsense are also seen when testing this on the cpu (i.e. sp.Device(-1)), so the non-deterministic behavior that is sometimes characteristic of gpu ops may not be the culprit

Hi, I finally get the chance to run through the notebook. I think the errors are within reasonable tolerance. In particular, when I look at element-wise difference I got np.abs(maps - old_maps).max(): 9.840061e-05

Because we are doing a lot of FFTs and iFFTs over data with high dynamic range (k-space centers are orders of magnitude larger than outer k-space), I think ~1e-4 is reasonable for normalized sensitivity maps (max is 1). For reproducibility, you can set a higher threshold, say 1e-3.

ad12 commented

Makes sense - could be good to add this as a note in the method. I can make a small PR for that

@ad12 Please refer to this when you make a PR! I'm closing this issue for now.