Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

accumulate_Hessian_times_input for Poisson...ProjData fails for input containing negatives #1461

Closed
KrisThielemans opened this issue Jun 19, 2024 · 3 comments · Fixed by #1477
Closed

accumulate_Hessian_times_input for Poisson...ProjData fails for input containing negatives #1461

KrisThielemans opened this issue Jun 19, 2024 · 3 comments · Fixed by #1477
Assignees
@KrisThielemans
Labels
bug
Milestone
v6.2

Comments

@KrisThielemans
Copy link
Collaborator

KrisThielemans commented Jun 19, 2024 •
edited
Loading

If input is filled with -1, the output is 0. This will be because we use divide_and_truncate (e.g. here). This function is designed to work for the PoissonLL operations, and truncates negatives to zero. That is incorrect for the multiplication with the Hessian though, as the "input" vector can be anything.

Unfortunately, it doesn't seem to be easy to make the Hessian calculation consistent with the truncation strategy used for the gradient (and value) of the PoissonLL (or at least, do it without becoming quite slow).
@KrisThielemans KrisThielemans self-assigned this Jun 19, 2024
@KrisThielemans KrisThielemans mentioned this issue Jul 4, 2024
Open
@KrisThielemans KrisThielemans added this to the v6.2 milestone Jul 9, 2024
@KrisThielemans
Copy link
Collaborator Author

value is computed as

STIR/src/buildblock/recon_array_functions.cxx

Lines 367 to 371 in 12bfa87

const float new_estimate = max(estimated_projections[r][b], projection_data[r][b] / max_quotient);
if (projection_data[r][b] <= small_value)
sub_result += -double(new_estimate);
else
sub_result += projection_data[r][b] * log(double(new_estimate)) - double(new_estimate);
. Gradient is consistent with that.

@KrisThielemans KrisThielemans mentioned this issue Jul 12, 2024
Open
@KrisThielemans
Copy link
Collaborator Author

Gradient is consistent with that.

not really, see #1472

@KrisThielemans KrisThielemans changed the title accumulate_Hessian_times_input fails for input containing negatives accumulate_Hessian_times_input for Poisson...ProjData fails for input containing negatives Jul 22, 2024
KrisThielemans added a commit to KrisThielemans/STIR-1 that referenced this issue Jul 22, 2024
@KrisThielemans
[BUGFIX] remove thresholding in Hessian for Poisson...ProjData
0835074 
accumulate_sub_Hessian_times_input_without_penalty used divide_and_truncate,
but this is incorrect as its threshold strategies are not appropriate
for (measured * forward_input)/forward_current^2. Most obvious example
is if forward_input contains negatives. We now use 0/x=0 (first thresholding
measured data to be non-negative).

UCL#1461
KrisThielemans added a commit to KrisThielemans/STIR-1 that referenced this issue Jul 22, 2024
@KrisThielemans
Hessian: throw error if result would be incorrect
7ab7390 
The Hessian calculations use divide_and_truncate, which sets negative
numerators to zero. This is incorrect when the forward projection
of the "input" is negative.
We now check this and throw an error if it occurs. A proper fix
will have to be for later.

See UCL#1461
KrisThielemans added a commit to KrisThielemans/STIR-1 that referenced this issue Jul 22, 2024
@KrisThielemans
Hessian times input: throw error if result would be incorrect
01daddd 
The Hessian calculations for the Poisson log-lik for projdata
use divide_and_truncate, which sets negative
numerators to zero. This is incorrect when the forward projection
of the "input" contains negatives.
We now check this and throw an error if it occurs. A proper fix
will have to be for later.

See UCL#1461
KrisThielemans added a commit to KrisThielemans/STIR-1 that referenced this issue Jul 22, 2024
@KrisThielemans
Hessian times input: throw error if result would be incorrect
2992bc5 
The Hessian calculations for the Poisson log-lik for projdata
use divide_and_truncate, which sets negative
numerators to zero. This is incorrect when the forward projection
of the "input" contains negatives.
We now check this and throw an error if it occurs. A proper fix
will have to be for later.

See UCL#1461
@KrisThielemans KrisThielemans mentioned this issue Jul 22, 2024
Merged
@KrisThielemans KrisThielemans linked a pull request Jul 22, 2024 that will close this issue
Hessian times input for PoissonLL...ProjData: throw error if result would be incorrect #1477
Merged
@KrisThielemans
Copy link
Collaborator Author

#1477 raises an error in this case (it doesn't really fix it though)

KrisThielemans added a commit to KrisThielemans/STIR-1 that referenced this issue Jul 25, 2024
@KrisThielemans
[BUGFIX] remove thresholding in Hessian for Poisson...ProjData
437b1e6 
accumulate_sub_Hessian_times_input_without_penalty used divide_and_truncate,
but this is incorrect as its threshold strategies are not appropriate
for (measured * forward_input)/forward_current^2. Most obvious example
is if forward_input contains negatives. We now use 0/x=0 (first thresholding
measured data to be non-negative).

UCL#1461
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@KrisThielemans KrisThielemans

Labels
bug
Projects
None yet
Milestone
v6.2
1 participant
@KrisThielemans