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.

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

Priors/Hessian method corrections #901

Merged
merged 28 commits into from
Jul 8, 2022

Conversation

robbietuk
Copy link
Collaborator

@robbietuk robbietuk commented Jun 7, 2021

This PR was going to simply add the accumulate_Hessian_times_input method to the RDP but while implementing it it was discovered that the implementations in the Quadratic and LogCosh priors were also incorrect. Additionally improves priors testing (#153).

This PR will attempt to do the following for each of the aforementioned priors:

  • Correct implementation issues with accumulate_Hessian_times_input to compute the correct Hessian vector product
  • Clearly document what is being done
  • Add a variable is_convex = false to generalised prior. Set is_convex = true for all convex priors.
  • Add test for convexity, much like Add a Hessian test to PoissonLogLikelihoodWithLinearModelForMeanAndProjData (and fix Hessian) #902 by testing x^T Hx >0
  • Add compute_Hessian() to all convex methods
  • Add test for convex priors by checking grad(x) - grad(x+eps) ~= H(x), where eps is a small perturbation.
  • Add diagonal_second_derivative and off_diagonal_second_derivative methods for use in all Hessian methods

Relative Difference Prior accumulate_Hessian_times_input

  • RDP accumulate_Hessian_times_input
  • Added RDP Hessian documentation
  • Add div by 0 safety
  • Add epsilon to Hessian function

@robbietuk
Copy link
Collaborator Author

These latest commits do not fix a possible mistake with the RDP's accumulate_Hessian_times_input but it does indicate convexity.

@robbietuk robbietuk changed the title Rdp/hessian vector product Priors/Hessian method corrections Jun 15, 2021
@robbietuk robbietuk marked this pull request as draft June 15, 2021 16:43
robbietuk added a commit to robbietuk/STIR that referenced this pull request Jun 18, 2021
…sian_concavity

This is a change to keep consistency with changes in UCL#901
@robbietuk
Copy link
Collaborator Author

robbietuk commented Jun 19, 2021

Implemented compute_Hessian for QP, logcosh, and RDP all the same style. This is complete code duplication for each of these methods in preperation for #694 changes.

I have rewritten compute_Hessian for the QP and Logcosh. Their implementations were not incorrect but required the assumption that d^2/dx^2 = - d^2/dxdy. The RDP does not abide by this constraint. I have also added the diagonal_second_derivative and off-diagonal_second_derivative methods (name may be changed) to each of the classes.

Added accumulate_Hessian_times_input to RDP and modernised the QP and logcosh priors to use the same methodology (code duplication again in prep for changes #694).

Added _is_convex and get _is_convex() to Generalised Prior. Default is false and each derived prior sets to true (if is the case)

Regarding testing:

  • Added Verbosity::set(0) calls to suppress compute_gradient info calls. When make test is run, this stuff is already suppressed but when executing test_priors the terminal output is flooded with gradient info calls without these verbosity changes.
  • Added test_gradient , test_Hessian_convexity,test_Hessian_convexity_configuration and test_Hessian_against_numerical methods. Hessian tests are only conducted if _is_convex=true

@robbietuk
Copy link
Collaborator Author

robbietuk commented Jun 19, 2021

Currently test_Hessian_against_numerical will fail for both RDP checks.

This is because of an inconsistancy between off-diagonal elements in the numerical-Hessian and the Hessian. The RDP is very sensitivie if eps (the perturbation) is too large relative to the comparison voxel values. Comparing the numerical-Hessian and the Hessian, the error is usually small (~10%), but can be large (>~60%), and occurs when both voxels values (x_j,x_k) are small w.r.t. eps, ergo effectively a large change in the numerical gradient.

I am fairly convinced that the implementation of the RDP's Hessian is correct and it is the test that is the issue.

Edit: Interesting that one of the travis jobs actually passed on test_priors. I guess that means that the numerical Hessian test for the RDP is only failing in certain input/target configurations. However, this is no stable. A constant eps in the test is not suitable for this prior.

robbietuk pushed a commit to robbietuk/STIR that referenced this pull request Jun 22, 2021
* Test PoissonLogLikelihoodWithLinearModelForMeanAndProjData Hessian

* Do subtraction in accumulate_sub_Hessian_times_input for Poisson

* Update Hessian_vector_product documentation and test error message

* Rename test_objective_function_Hessian to test_objective_function_Hessian_concavity

This is a change to keep consistency with changes in UCL#901
@robbietuk
Copy link
Collaborator Author

Create bools to toggle the use of various tests for different penalties. RDP with RDP_eps = 0 will likely allways fail the numerical Hessian test.

@robbietuk robbietuk marked this pull request as ready for review June 24, 2021 22:25
@robbietuk
Copy link
Collaborator Author

robbietuk commented Jul 15, 2021

Ready for review. There is overlap with #902.

Copy link
Collaborator

@KrisThielemans KrisThielemans left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

some comments. I didn't check the tests yet

src/include/stir/recon_buildblock/GeneralisedPrior.h Outdated Show resolved Hide resolved
src/include/stir/recon_buildblock/GeneralisedPrior.h Outdated Show resolved Hide resolved
src/include/stir/recon_buildblock/GeneralisedPrior.h Outdated Show resolved Hide resolved
src/include/stir/recon_buildblock/GeneralisedPrior.h Outdated Show resolved Hide resolved
src/include/stir/recon_buildblock/LogcoshPrior.h Outdated Show resolved Hide resolved
src/include/stir/recon_buildblock/LogcoshPrior.h Outdated Show resolved Hide resolved
src/recon_buildblock/LogcoshPrior.cxx Outdated Show resolved Hide resolved
src/include/stir/recon_buildblock/LogcoshPrior.h Outdated Show resolved Hide resolved
src/include/stir/recon_buildblock/LogcoshPrior.h Outdated Show resolved Hide resolved
@robbietuk
Copy link
Collaborator Author

some comments. I didn't check the tests yet

I have addressed all pertinent suggestions

@robbietuk
Copy link
Collaborator Author

Future progress on this PR should include an intermediate class for the impacted priors. This class should handle the for loops for each of the Hessian methods and create the virtual functions for the derivative_20 and derivative_11.

robbietuk and others added 13 commits June 24, 2022 16:00
* RDP Documententation

* QP use (off_)diagonal_second_derivative methods and documentation

* Restructure RDP and QP Hessian logic

* Improve RDP and QP documentation

* Implement log cosh accumulate_Hessian_times_input in terms of derivatives

* Documentation
* Add method get_is_convex(), which accesses prior is_convex parameter and only test convexity of priors if convex function

* compute_Hessian method added to generalised prior and errors by default but with different messages depending on is_convex 

* Test rename and create test_Hessian_convexity and test_Hessian_methods

* Add test for compute_Hessian

* RelativeDifferencePrior initialisation call set_defaults

* Modernise `compute_Hessian` for QP and LogCosh. Add `compute_Hessian` for RDP

* Correct RDP second derivative functions

* Major changes to test_Hessian_against_numerical, which loops over each voxel for perturbation response

* Add verbosity suppression to suppress gradient info calls
Reenable PLS, but disable all numerical tests - just run setup.
Make is_convex() a pure virtual method in GeneralisedPrior to be implemented in each prior.

FilterRoot is not convex as it does not have 0th or 2nd order behaviour
Co-authored-by: Kris Thielemans <KrisThielemans@users.noreply.github.com>
Co-authored-by: Kris Thielemans <KrisThielemans@users.noreply.github.com>
This removes two checks on Succeeded::no in GeneralisedObjectiveFunction
@robbietuk robbietuk force-pushed the RDP/HessianVectorProduct branch from 15aebe7 to 6596d29 Compare June 27, 2022 21:20
Copy link
Collaborator Author

@robbietuk robbietuk left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is ready to go, assuming the tests pass.

Copy link
Collaborator

@KrisThielemans KrisThielemans left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

looks good. One tiny comment. Also, can you add something to release_5.1.htm? You probably will have to merge master into here first to get that file.

src/recon_buildblock/GeneralisedPrior.cxx Outdated Show resolved Hide resolved
@KrisThielemans KrisThielemans merged commit 4d087f3 into UCL:master Jul 8, 2022
markus-jehl pushed a commit to markus-jehl/STIR that referenced this pull request Jul 25, 2022
…oduct"

This reverts commit 4d087f3, reversing
changes made to f98246f.
markus-jehl pushed a commit to markus-jehl/STIR that referenced this pull request Aug 18, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants