KullbackLeibler gradient

[epapoutsellis]

At the moment, there are two implementations of the gradient of the KullbackLeibler

Numba: https://github.com/TomographicImaging/CIL/blob/701c4999c53e77c6ae3b587b1776f74accff2335/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py#L123-L130

and Numpy

https://github.com/TomographicImaging/CIL/blob/701c4999c53e77c6ae3b587b1776f74accff2335/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py#L344-L353

The numpy implementation above is not correct, see #1269 and will be corrected by PR #1272.

However, in the case where x and eta have zero values, numba implementation fails as

File ~/Edo/build/INSTALL/python/cil/optimisation/functions/KullbackLeibler.py:487, in KullbackLeibler_numba.gradient(self, x, out)
    485 else:
    486     out_np = out.as_array()
--> 487     kl_gradient(x.as_array(), self.b.as_array(), out_np, self.eta.as_array())
    488     out.fill(out_np)

ZeroDivisionError: division by zero

Reproduce: NEMA phantom + KullbackLeibler + FISTA.

[paskino]

So this is issue is now Handle ZeroDivisionError in KullbackLeibler gradient method? Division by zero originate when both x and eta have zeros in the same location.

https://github.com/TomographicImaging/CIL/blob/701c4999c53e77c6ae3b587b1776f74accff2335/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py#L123-L130

[KrisThielemans]

There have been changes to the code so the above isn't up-to-date anymore. Current status:

Documentation

https://github.com/TomographicImaging/CIL/blob/2f92286b00381e703f4871e14b97cd87be7ea4ad/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py#L37-L44 (Note that the last line of this comment is superfluous, it's handled in the "ifs" above).

numpy implementation:

https://github.com/TomographicImaging/CIL/blob/2f92286b00381e703f4871e14b97cd87be7ea4ad/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py#L138-L143 Therefore, this will actually return 0 when u>0, v=0, which is not what was documented at the start, and strictly speaking incorrect for KL. (Note that the gradient is consistent with the above choice)

numba implementation

https://github.com/TomographicImaging/CIL/blob/2f92286b00381e703f4871e14b97cd87be7ea4ad/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py#L350-L360 Does return infinity.

conclusion

numby and numpy implementations are still inconsistent. The numpy implementation is a practical choice that avoids problems, and is also what STIR/SIRF does. I don't really like it, but it does happen for instance when the reconstructed image size is smaller than the FOV. The forward projections in the "outer" part of the sinogram will therefore be zero, while the measured data might not be. This would cause problems tracking KL improvement. (Not any actual problems in the reconstruction, as the backprojection of those LORs will never contribute to the image). Of course, in real life, the scatter/randoms estimate should not be zero, and therefore this situation shouldn't ever occur in the "proper" case, but people do crazy things...

I'd be fine with either choice (return inf or 0), but I would expect both implementations to give the same result, and to match the documentation.

Maybe @MargaretDuff has an opinion...