Pdhg strongly convex

[epapoutsellis]

This PR adds primal OR dual acceleration for the PDHG algorithm in the case where the function f^{*} or g is strongly convex. These are **kwargs** arguments in the PDHG algorithm,i.e., gamma_g for the strongly convex constant of the function g (primal acceleration) and gamma_fconj for the strongly convex constant of the convex conjugate of f (dual acceleration). In particular:

• Checks if the strongly convex constants are positive Numbers.
• Always raises a warning for the users to be sure that the number is correct and related to the strongly convex property of the function.
• Creates a new method update_step_sizes which is called at the end of the update method. Therefore, the step-sizes change in every iteration, according to the specific case, f^{*} or g being strongly convex . If gamma_g or gamma_fconj are None, the step-sizes remain unchanged.
• Error is raised if the user passes strongly convex constants for both functions. This adaptive rule is not implemented at the moment.
• Creates a new method set_step_sizes to check the correct type of the primal-dual step-sizes, sigma and tau and set the default scalar values.
• Check if sigma, tau are positive Numbers, if they are array-like objects with the correct shape determined by the operator. If the user passes sigma and tau with the correct shape, there is no guarantee that the algorithm will converge unless it follows Lemma 2.
  - If the shape is correct, the PDHG algorithm will not run, unless use_axpby=False. Fixed by axpby blackdatacontainer fix #1080 .
• Raises a warning if the convergence condition for the case of scalar step-sizes is not satisfied.

Note The adaptive rule for primal or dual acceleration is described in Algorithm 2.

1. A denoising example can be found here where the strong convexity property comes from the function L2NormSquared.

2. A tomography example can be found here where the strong convexity property comes from the convex conjugate of the BlockFunction f. Not that the ProjectionOperator is with device=cpu to avoid any mismatch on the adjoint.

[paskino]

I guess we need to:

1. document the additional variables gamma_g and gamma_fconj in the docstring
2. add unittests (basic, not just one full denoising example). Something like set up a PDHG without the gammas, do 1 iteration and check that the sigma and tau remain the same. Set up a PDHG with the gammas, do 1 iteration and check that the sigma and tau have been changed to what you expect.

[paskino]

import numpy as np

[epapoutsellis]

done

[jakobsj]

Great to get a strong convexity version in but I'm not sure about the parameters being given as additional parameters to PDHG. Wouldn't they belong more naturally as attributes of the Functions, similar to Lipschitz? With default being None and automatically set on initialising the Function whenever possible? Then PDHG would check the f and g Functions for these parameters and if different from None the strongly convex version of the algorithm would be used.

Then the user wouldn't need to worry about understanding what these parameters are, computing them (how would they do that?) and providing them as input.

Not sure if this approach has been discussed?

[paskino]

@jakobsj 's got a great point here. If we put the gammas as properties of Functions we will be able to:

1. immediately provide a working implementation by manually setting such variables in the function
2. such implementation can be further developed when we figure out how to compute the gammas and add this logic to the relevant function classes.

[epapoutsellis]

That was my initial plan, but we decided to not follow this route. Basically, we cannot do it atm. As the convex_conjugate returns a number and not a function.

In this example, we use the strongly convexity case of PDHG, when the conjugate of the BlockFunction , i.e., separable function is strongly convex. That is the sum of the convex_conjugate of L2NormSquared and the convex_conjugate of the MixedL21Norm.

If you run the example above, you will see the following figure (label is wrong, **should be Strongly Convex (fconj) **)

[jakobsj]

Yes I remember the discussion now. But I don't see how it's not possible. I don't think we need the convex conjugate to be returned as a function. We just need the function itself to have as an attribute the strong convexity parameter of its convex conjugate.

…

On Fri, 29 Oct 2021, 09:08 Vaggelis Papoutsellis, ***@***.***> wrote: That was my initial plan, but we decided to not follow this route. Basically, we cannot do it atm. As the convex_conjugate returns a number and not a function. In this example <https://gist.github.com/epapoutsellis/f127dcd497b540e06ca3cfb3f6c32242#file-pdhg_strongly_convex_tomography-py-L84>, we use the strongly convexity case of PDHG, when the conjugate of the BlockFunction , i.e., separable function is strongly convex. That is the sum of the convex_conjugate of L2NormSquared and the convex_conjugate of the MixedL21Norm. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1030 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACMDSCDARP5EACG36N6WX33UJJJA3ANCNFSM5GYD4HFQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[epapoutsellis]

@gfardell I prefer to have the methods update, set_step_sizes, update_step_sizes, update_objective in the docs and not use _ because I believe they are important and if the user does not want to read the docs and see what is going on in these methods, can click on source and have a look at the code.

[gfardell]

I think the example code is my biggest issue. It's unmaintainable as it is, and worse than having the code snippets in __main__ as we used to. Can the examples be restricted to how to set up, modify and use the class. We have CIL demos for full CIL usage.

[gfardell]

@gfardell I prefer to have the methods update, set_step_sizes, update_step_sizes, update_objective in the docs and not use _ because I believe they are important and if the user does not want to read the docs and see what is going on in these methods, can click on source and have a look at the code.

I don't think they'll know the difference between set_step_sizes and update_step_sizes. if one is marked as _ then they can see it (you can even include it in the docs if you want), but they won't modify it. But we have not used this elsewhere in algorithms so we can leave it for now.

[epapoutsellis]

I removed the Example section and refer to the CIL-Demos repo and the papers. The default sigma and tau step-sizes are described explicitly in the beginning of the docs.

CIL/Wrappers/Python/cil/optimisation/algorithms/PDHG.py

Lines 128 to 146 in deadd5b

	- By default, the step sizes :math:`\sigma` and :math:`\tau` are positive scalars and defined as below:
	* If ``sigma`` is ``None`` and ``tau`` is ``None``:
	.. math::
	\sigma = \frac{1}{\|K\|}, \tau = \frac{1}{\|K\|}
	* If ``tau`` is ``None``:
	.. math::
	\tau = \frac{1}{\sigma\|K\|^{2}}
	* If ``sigma`` is ``None``:
	.. math::
	\sigma = \frac{1}{\tau\|K\|^{2}}

[paskino]

use assertAlmostEqual

[paskino]

move this check into the setter

[paskino]

A few minor changes