Non repeatable behavior of Image.backproject

[remy-abergel]

Describe the bug:

Calling Image.backproject several times on the same input signal does not lead to the same output.

To Reproduce:

#####################
# loading libraries #
#####################
import numpy as np
import aspire
from aspire.image import Image
from aspire.utils import Rotation
from aspire.volume import Volume

##################################################################
# compute yref as a set of n_img random images with size LxL and #
# back-project it                                                #
##################################################################
n_img,L = 100,32
yref = aspire.image.Image(np.array(np.random.rand(n_img,L,L),dtype='float32'))
rots = Rotation.generate_random_rotations(n=n_img, seed=31415)
vref = yref.backproject(rots.matrices)

#####################################################################
# repeat the back-projection of yref several times and compare each #
# output back-projected volume to vref --> the diffrence should be  #
# strictly zero (not satisfied here)                                #
#####################################################################
Nsimu = 1000
maxabserr = -np.inf
for id in range(Nsimu):
    v = yref.backproject(rots.matrices)
    abserr = np.max(np.abs(v.asnumpy()[0]-vref.asnumpy()[0]))
    maxabserr = np.max([abserr,maxabserr])
    print("simulation %d/%d : absolute error = %.3e, max. absolute error = %.3e" % (id+1,Nsimu,abserr,maxabserr))

Expected behavior:

The absolute error between v and vref should be exactly zero at each simulation run. On my laptop, I end up with small (but nonzero) errors roughly equal to $5\cdot 10^{-7}$.

Environment:

• OS: Ubuntu 22.04
• Python: 3.10.12
• ASPIRE Version: 0.12.1 dev
• CPU Type: Intel® Core™ i7-7920HQ CPU @ 3.10GHz × 8
• GPU Extensions: yes

Additional context:

I applied the same test to Volume.project leading to a perfectly repeatable behavior for this operator (see script below)

############################################################
# same experiment with the volume projection operator (OK) #
############################################################
_vref = aspire.volume.Volume(np.array(np.random.rand(L,L,L),dtype='float32'))
_yref = _vref.project(rots)
Nsimu = 1000
maxabserr = -np.inf
for id in range(Nsimu):
    _y = _vref.project(rots)
    abserr = np.max(np.abs(_y.asnumpy()-_yref.asnumpy()))
    maxabserr = np.max([abserr,maxabserr])
    print("simulation %d/%d : absolute error = %.3e, max. absolute error = %.3e" % (id+1,Nsimu,abserr,maxabserr))

At first sight, the non-repeatability issue with Image.backproject does not seem to be related to potential CPU parallelization (I observed the same non-repeatable behavior when working with only one CPU). Is there any source of randomness occurring during the Image.backproject operation?

[garrettwrong]

The ASPIRE code for this is basically just FFT and NUFFT. I think that non-determinism may come primarily from the finufft call itself. We have brought up determinism to them a few times, but it has generally been noted as within their expectations. IE, there is not a guarantee of run to run repro. You can run in doubles with a smaller epsilon to bring the errs down to -e15 or so. However, the default we use is already just beyond the max eps for float (you can see I dial it back slightly in the code below).

Using the debugger, I dropped out the arrays in your scripts anufft call, like so

(Pdb) np.save('im_f.npy', im_f)
(Pdb) np.save('pts_rot', pts_rot)
(Pdb) L
32

Then in the following script the non-determinism is evident in only a few runs on my machine.

import numpy as np

from aspire.nufft import anufft
import finufft


im_f = np.load("im_f.npy")
pts_rot = np.load("pts_rot.npy")
L = 32
NN = 10

ma = 0
# ASPIRE wrapper of FINUFFT
eps = 1e-8  # default
res0 = anufft(im_f, pts_rot, (L, L, L), epsilon=eps, real=True)
for i in range(NN):
    res = anufft(im_f, pts_rot, (L, L, L), real=True)
    a = np.max(np.abs(res - res0))
    ma = max(a, ma)
    print(f"iter {i:3} : {a}")

print(f"ASPIRE wrapper Max Abs Err: {ma}")

# Direct FINUFFT call
cdtype = np.complex64  # complex single prec
eps = max(eps, np.finfo(pts_rot.dtype).eps)
plan = finufft.Plan(
    nufft_type=1,
    n_modes_or_dim=(L, L, L),
    eps=eps,
    n_trans=1,
    dtype=pts_rot.dtype,
)

# enforce pts convention
pts = np.ascontiguousarray(np.mod(pts_rot + np.pi, 2 * np.pi) - np.pi)

plan.setpts(*pts)

ma = 0
for i in range(NN):
    fires = np.real(plan.execute(im_f.astype(cdtype)))
    a = np.max(np.abs(fires - res0))
    ma = max(a, ma)
    print(f"iter {i:3} : {a}")

print(f"Direct Max Abs Err: {ma}")

 python anufft_det.py 
iter   0 : 1.9073486328125e-05
iter   1 : 1.52587890625e-05
iter   2 : 1.9073486328125e-05
iter   3 : 1.52587890625e-05
iter   4 : 1.1444091796875e-05
iter   5 : 1.52587890625e-05
iter   6 : 1.9073486328125e-05
iter   7 : 1.52587890625e-05
iter   8 : 1.52587890625e-05
iter   9 : 2.288818359375e-05
ASPIRE wrapper Max Abs Err: 2.288818359375e-05
iter   0 : 1.52587890625e-05
iter   1 : 1.52587890625e-05
iter   2 : 1.52587890625e-05
iter   3 : 2.288818359375e-05
iter   4 : 2.288818359375e-05
iter   5 : 1.52587890625e-05
iter   6 : 1.52587890625e-05
iter   7 : 1.9073486328125e-05
iter   8 : 1.52587890625e-05
iter   9 : 1.52587890625e-05
Direct Max Abs Err: 2.288818359375e-05

The FINUFFT developers may have some tricks to improve this if you are interested in pursuing it with them. Let me know what you think.

Oh and thanks again for the nice report, it has been a pleasure looking into your concerns... and I really don't get to say that often :).

[remy-abergel]

Thank you for having isolated the issue to the finufft modules.

I do believe that the error magnitude is small enough for not being an issue in my own applications. However, I do not understand the reason for such non-determinism in this nonuniform FFT library.

Being more familiar with the TU Chemnitz NFFT library, I never observed non-deterministic behavior in the evaluation of NFFT and its adjoint.

Quickly looking at the reference papers related to fiNUFFT, I did not find any mention to stochastic algorithms that could have explained the non-determinism. My first guess would be an incorrect parallelization issue, I will investigate further this point.

[remy-abergel]

I think there is indeed a parallelization issue. When launching python using

taskset -c 0 python

then running again

#####################
# loading libraries #
#####################
import numpy as np
import aspire
from aspire.image import Image
from aspire.utils import Rotation
from aspire.volume import Volume

##################################################################
# compute yref as a set of n_img random images with size LxL and #
# back-project it                                                #
##################################################################
n_img,L = 100,32
yref = aspire.image.Image(np.array(np.random.rand(n_img,L,L),dtype='float32'))
rots = Rotation.generate_random_rotations(n=n_img, seed=31415)
vref = yref.backproject(rots.matrices)

#####################################################################
# repeat the back-projection of yref several times and compare each #
# output back-projected volume to vref --> the diffrence should be  #
# strictly zero (not satisfied here)                                #
#####################################################################
Nsimu = 1000
maxabserr = -np.inf
for id in range(Nsimu):
    v = yref.backproject(rots.matrices)
    abserr = np.max(np.abs(v.asnumpy()[0]-vref.asnumpy()[0]))
    maxabserr = np.max([abserr,maxabserr])
    print("simulation %d/%d : absolute error = %.3e, max. absolute error = %.3e" % (id+1,Nsimu,abserr,maxabserr))

I end up with a perfectly repeatable behavior. I will report the issue to the finufft team.

[garrettwrong]

Hi, I doubt it is due to the algorithm, more implementation. I think this may be related to threading or vectorization (or both).

When codes are compiled with threading, or with advanced vector optimization they often drop guarantees about dispatch orders etc, and this can have small numerical effects. You may also try export OMP_NUM_THREADS=1 and/or compiling their code from source with different compiler optimization options if you want to explore this more. In other codes, going back to -msseX usually fixed most vector repro issues (by disabling non-determism caused by mixing with newer AVX instructions), but again, losing performance. Threading is a different animal :).

I'll follow that other issue, best of luck.

[remy-abergel]

Your comment was perfectly right, this issue was not a good one.
Best,
Rémy