[SyneRBI/SIRF] Implement adding Hessian diagonal in SIRF (PR #1370)
Evgueni Ovtchinnikov
Changes in this pull request
Testing performed
Related issues
Fixes #1369
Checklist before requesting a review
- I have performed a self-review of my code
- I have added docstrings/doxygen in line with the guidance in the developer guide
- I have implemented unit tests that cover any new or modified functionality
- The code builds and runs on my machine
- CHANGES.md has been updated with any functionality change
Contribution Notes
Please read and adhere to the contribution guidelines.
Please tick the following:
- The content of this Pull Request (the Contribution) is intentionally submitted for inclusion in SIRF (the Work) under the terms and conditions of the Apache-2.0 License.
You can view, comment on, or merge this pull request online at:
https://github.com/SyneRBI/SIRF/pull/1370
Commit Summary
- d95c7bb updated files for release 3.9.0 [ci skip]
- 54505e7 Merge branch 'master' of https://github.com/SyneRBI/SIRF
- dedf83c Merged branch 'master' of https://github.com/SyneRBI/SIRF
- 17867e4 Merge branch 'master' of https://github.com/SyneRBI/SIRF
- 5a1e82f Merge branch 'master' of https://github.com/SyneRBI/SIRF
- 0980bdd started working on adding Hessian diagonal in SIRF
File Changes
(1 file)
Patch Links:
—
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 pushed 1 commit.
- ae578b8 added constructor for xSTIR RDP from Gibbs RDP
—
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#1370)I think the best strategy is to just use GibbsRelativeDifferencePenalty everywhere. It should be backwards compatible with RelativeDifferencePrior. @ColomboMatte0 is that correct? (we don't need the surrogate stuff here AFAIK)
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Matteo Colombo
ColomboMatte0 left a comment (SyneRBI/SIRF#1370)It is backwards compatible with all the methods implemented in RelativedifferencePrior , we do not have SPS stuff for the RDP.
—
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 pushed 1 commit.
- c9ee3cc trying Kris' suggestion to derive SIRF RDP from STIR Gibbs penalty
—
View it on GitHub or unsubscribe.
You are receiving this because you are subscribed to this thread.
Evgueni Ovtchinnikov
@evgueni-ovtchinnikov pushed 1 commit.
- b67e1b0 replaced STIR RDP with SIRF RDP derived from Gibbs penalty
—
View it on GitHub or unsubscribe.
You are receiving this because you are subscribed to this thread.
Evgueni Ovtchinnikov
@evgueni-ovtchinnikov pushed 1 commit.
- 77754b4 added method compute_Hessian_gradient to SIRF RelativeDifferencePrior
—
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#1370)@ColomboMatte0 can you suggest a test for this functionality?
—
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 pushed 1 commit.
- bcb5a6e removed some commented-out lines [ci skip]
—
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#1370)@ColomboMatte0 could you please suggest a Python test for this functionality. Thank you, E.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Matteo Colombo
ColomboMatte0 left a comment (SyneRBI/SIRF#1370)sorry i lost the previous notification, you mean a test for the hessian diagonal? not really sure what you mean with hessian gradient
—
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#1370)Yes. For the diagonal (not the gradient). Not sure if we have one in STIR that could be replicated?
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Matteo Colombo
ColomboMatte0 left a comment (SyneRBI/SIRF#1370)i remember i tested it numerically using finite differences from the gradient but probably those tests were not added in my PR, i can add it tomorrow, shouldn't be a big deal.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Matteo Colombo
ColomboMatte0 left a comment (SyneRBI/SIRF#1370)$$ H_{ii} \approx \frac{\nabla_i f(x + \epsilon e_i) - \nabla_i f(x - \epsilon e_i)}{2\epsilon} $$
you can test the correctness by computing this, let me know if you want me to implement it somewhere.
—
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#1370)Good idea. I suggest to use
$$ H_{ii} \approx \frac{\nabla_i f(x + \epsilon e_i) - \nabla_i f(x)}{\epsilon} $$
a bit easier to implement, and also avoids potential issues with the RDP which is only defined on non-negative images.
—
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#1370)I doubt @ColomboMatte0 will have the time for this, unfortunately. Shouldn't be too hard to build on our current "gradient of a function" tests. Just one more loop...
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Matteo Colombo
ColomboMatte0 left a comment (SyneRBI/SIRF#1370)i can do it. but shall i open a new quick PR for STIR? or you want me to add the test where exactly?
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Matteo Colombo
ColomboMatte0 left a comment (SyneRBI/SIRF#1370)here src/xSTIR/pSTIR/tests/tests_qp_lc_rdp.py?
—
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#1370)@ColomboMatte0 yes, that's a good place, please go ahead with a quick PR. Thanks a lot!
—
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 pushed 1 commit.
- eb8edbb Add numerical test for Hessian diagonal (#1385)
—
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#1370)Test failures:
src/xSTIR/pSTIR/tests/tests_qp_lc_rdp.py:97: in Hessian_diagonal_test
test.check_if_less(abs(numerical - analytical) / max_H, 0.01)
...
E ValueError: +++ test 29 failed: np.float32(0.011938471) >= 0.01
Bounds are always very hard to set in this stuff. @ColomboMatte0 what do you think? (Did you run it locally?)
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Matteo Colombo
ColomboMatte0 left a comment (SyneRBI/SIRF#1370)yes I tested it locally on a dataset and it worked on my end. I suspect the CI is running the test on a different dataset. A quick workaround could be to slightly increase the tolerance for this specific test, given that we are using a one-sided derivative.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you are subscribed to this thread.
Matteo Colombo
ColomboMatte0 left a comment (SyneRBI/SIRF#1370)@KrisThielemans I ran the test locally with the CI dataset and reproduced the failure — the one-sided first-order finite difference has O(ε) truncation error which exceeds 1% on some voxels.
Two options:
Option A — relax tolerance to 2% (already in #1386)
Option B — use a second-order one-sided finite difference, keeping the 1% tolerance:
$$ H_{ii} \approx \frac{-3 , \nabla_i f(x) + 4 , \nabla_i f(x + \varepsilon , e_i) - \nabla_i f(x + 2\varepsilon , e_i)}{2\varepsilon} $$
This is O(ε²) instead of O(ε), still one-sided (safe for RDP with non-negative images), but requires one extra gradient evaluation per test voxel.
Let me know which you prefer and I'll update the PR.
—
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#1370)thanks. Although I like the second order formulation, I think it's probably not so easy to get right with fast changing priors such as the RDP around zero. I think the main aim here is to have a consistency check and indeed code check (as indeed you found some!).
I'll therefore accept your PR.
Obviously, if we don't have this in STIR yet, it would be useful to add it there as well (sorry for the duplication!). Feel free to try the 2nd order one there :-)
thanks again
—
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 pushed 1 commit.
- b459171 Relax Hessian diagonal test tolerance from 1% to 2%
—
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.