[SyneRBI/SIRF] Poisson Log Likelihood ObjectiveFunction `gradient` with and without `out` parameter has larger difference than expected (Issue #1349)
Edoardo Pasca
paskino created an issue (SyneRBI/SIRF#1349)With a Poisson Log Likelihood ObjectiveFunction one may calculate the gradient at a certain image x by calling gradient. If a pre-initialised buffer is passed in the out parameter, this is filled with the result of the gradient calculation. If no out parameter is passed the function allocates and returns the result in a newly allocated buffer. However, comparing these 2 lead to larger difference than expected.
obj_fun = pet.make_Poisson_loglikelihood(acquired_data)
g1 = obj_fun.gradient(x)
g2 = g1 * 0
obj_fun.gradient(x, out=g2)
numpy.testing.assert_allclose(g2.asarray(), g1.asarray(), atol=1e-5)
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Kris Thielemans
KrisThielemans left a comment (SyneRBI/SIRF#1349)There could be a difference in float<->double possibly? STIR works with floats. I don't think we're converting though, but @evgueni-ovtchinnikov can confirm.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Evgueni Ovtchinnikov
evgueni-ovtchinnikov left a comment (SyneRBI/SIRF#1349)@paskino have you checked dtype of g1.asarray() and g2.asarray()?
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
MaxTousss
MaxTousss left a comment (SyneRBI/SIRF#1349)Greetings,
While doing some tests for the PETRIC2 challenge, I might have stumbled upon a similar problem. After some investigation, it seems that the backward operator is slightly unstable. One can observe this by using the BSREM2 algorithm provided in SIRF contrib and computing (Note, self is the class BSREM2 for a given dataset of PETRIC2)
np.sum((self.subset_gradient(self.x, 0) - self.subset_gradient(self.x, 0)).asarray().astype('float')**2) = 0.0013965924207424846
and
np.sum((self.acq_models[0].backward(self.data[0]) - self.acq_models[0].backward(self.data[0])).asarray().astype('float')**2) = 2.0470940317509713e-06
(Value were obtained from one of the dataset provided in PETRIC2 and mostly serve to show the "amplitude", i.e., slight difference, of the problem).
Note that my attempt to use the out feature of the backward method resulted in erratic testing. For example,
b1 = self.acq_models[0].backward((self.data[0]))
b2 = self.acq_models[0].backward((self.data[0]))
self.acq_models[0].backward((self.data[0]), out=b2)
self.acq_models[0].backward((self.data[0]), out=b1)
np.sum((b2-b1).asarray().astype('float')**2)
> 1.0141001191048331e-05
self.acq_models[0].backward((self.data[0]), out=b1)
np.sum((b2-b1).asarray().astype('float')**2)
> 2187614105.957749
This was tested in the docker provided for the PETRIC2 challenges.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Evgueni Ovtchinnikov
evgueni-ovtchinnikov left a comment (SyneRBI/SIRF#1349)@paskino @casperdcl Just as I suspected, this is a typical trouble with parallel computation.
I added these lines to osem_reconstruction.py after line 112:
import numpy
x = image
obj_fun.set_up(x)
g1 = obj_fun.gradient(x)
g2 = g1 * 0
obj_fun.gradient(x, out=g2)
d = g2 - g1
print(numpy.linalg.norm(d.asarray()))
exit()
Each time I ran it I got slightly different norm of d.asarray(), e.g.:
0.017808886
0.015998868
After exporting OMP_NUM_THREADS=1, I got 0.0 each time.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Kris Thielemans
KrisThielemans left a comment (SyneRBI/SIRF#1349)I agree we'll get variations due to multithreading (it's hard to know how large these are from the numbers above, as we don't know what the norm is on f1, is course).
However, wasn't the main problem observed by @MargaretDuff that we have to to first fill g2 with 0? I have no idea what the usual numpy convention is, but suspect numpy will overwrite.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Casper da Costa-Luis
casperdcl left a comment (SyneRBI/SIRF#1349)Wait what do you mean by multithreading? Are we not using atomic operations where appropriate to ensure consistent results?
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Kris Thielemans
KrisThielemans left a comment (SyneRBI/SIRF#1349)Calling parallelproj backprojection via STIR will likely (i.e., I didn't check) do some adding of backprojected images in CPU via OpenMP. Of course, it will do stuff in parallel on the GPU. Both avoid race conditions by atomic operations. That doesn't mean that results are consistent though, once you do floats. This is because a+b+c != c+b+a != a + c +b ....
Example, if you need to add a list of numbers with a varying magnitude, you get the best numerical precision by first adding the smallest ones.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Edoardo Pasca
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Kris Thielemans
KrisThielemans left a comment (SyneRBI/SIRF#1349)Are you sure this can be closed? what about
However, wasn't the main problem observed by @MargaretDuff that we have to first fill g2 with 0? I have no idea what the usual numpy convention is, but suspect numpy will overwrite.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Kris Thielemans
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Kris Thielemans
KrisThielemans left a comment (SyneRBI/SIRF#1349)see https://discord.com/channels/1242028164105109574/1252651306779152535/1458465128117833728
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Evgueni Ovtchinnikov
evgueni-ovtchinnikov left a comment (SyneRBI/SIRF#1349)@KrisThielemans see my understanding of Margaret's problem suggested 08/01/2026 at 17:44.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Kris Thielemans
KrisThielemans left a comment (SyneRBI/SIRF#1349)Do you mean #1349 (comment)? That is about numerical precision. Margaret's problem doesn't look like that
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Evgueni Ovtchinnikov
evgueni-ovtchinnikov left a comment (SyneRBI/SIRF#1349)@KrisThielemans no, I refer specifically to my answer to @MargaretDuff of 08/01/2026 at 17:44. She never replied, so I take it she found my suggestions plausible.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Kris Thielemans
KrisThielemans left a comment (SyneRBI/SIRF#1349)Please quote it here.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Evgueni Ovtchinnikov
evgueni-ovtchinnikov left a comment (SyneRBI/SIRF#1349)For each i, x_out is overwritten showing us only what subset i of detector pairs "sees".
My understanding is that images of data subsets should be summed up, accumulating in x_tmp.
We could get better picture if you printed norms of x_out and x_tmp for each i.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.