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

Fix kwargs in PDHG #2010

Merged
merged 7 commits into from
Dec 13, 2024
Merged

Fix kwargs in PDHG #2010

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
d913589
Remove kwargs from PDHG
MargaretDuff Dec 11, 2024
01792d1
Updated based on discussions with Gemma
MargaretDuff Dec 12, 2024
52ecd89
Reference to theta
MargaretDuff Dec 12, 2024
45e0640
Updates from Laura's review
MargaretDuff Dec 13, 2024
e53731b
Merge branch 'master' into check-convergence-in-pdhg
MargaretDuff Dec 13, 2024
88d085b
Changelog and Gemma's comments
MargaretDuff Dec 13, 2024
4215402
Small error
MargaretDuff Dec 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 34 additions & 21 deletions Wrappers/Python/cil/optimisation/algorithms/PDHG.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,26 +38,25 @@ class PDHG(Algorithm):
A convex function with a "simple" proximal.
operator : LinearOperator
A Linear Operator.
sigma : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default=None
sigma : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default is 1.0/norm(K) or 1.0/ (tau*norm(K)**2) if tau is provided
Step size for the dual problem.
tau : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default=None
tau : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default is 1.0/norm(K) or 1.0/ (sigma*norm(K)**2) if sigma is provided
Step size for the primal problem.
initial : DataContainer, optional, default=None
initial : DataContainer, optional, default is a DataContainer of zeros
Initial point for the PDHG algorithm.
gamma_g : positive :obj:`float`, optional, default=None
Strongly convex constant if the function g is strongly convex. Allows primal acceleration of the PDHG algorithm.
gamma_fconj : positive :obj:`float`, optional, default=None
Strongly convex constant if the convex conjugate of f is strongly convex. Allows dual acceleration of the PDHG algorithm.

**kwargs:
Keyward arguments used from the base class :class:`Algorithm`.

max_iteration : :obj:`int`, optional, default=0
Maximum number of iterations of the PDHG.
update_objective_interval : :obj:`int`, optional, default=1
MargaretDuff marked this conversation as resolved.
Show resolved Hide resolved
Evaluates objectives, e.g., primal/dual/primal-dual gap every ``update_objective_interval``.
check_convergence : :obj:`boolean`, default=True
Checks scalar sigma and tau values satisfy convergence criterion
Checks scalar sigma and tau values satisfy convergence criterion and warns if not satisfied. Can be computationally expensive for custom sigma or tau values.
theta : Float between 0 and 1, default 1.0
Relaxation parameter for the over-relaxation of the primal variable.


Example
-------
MargaretDuff marked this conversation as resolved.
Show resolved Hide resolved
MargaretDuff marked this conversation as resolved.
Show resolved Hide resolved
Expand Down Expand Up @@ -222,7 +221,16 @@ class PDHG(Algorithm):
"""

def __init__(self, f, g, operator, tau=None, sigma=None, initial=None, gamma_g=None, gamma_fconj=None, **kwargs):
"""Initialisation of the PDHG algorithm"""

self._theta = kwargs.pop('theta', 1.0)
MargaretDuff marked this conversation as resolved.
Show resolved Hide resolved
if self.theta>1 or self.theta<0:
raise ValueError("The relaxation parameter theta must be in the range [0,1], passed theta = {}".format(self.theta))

self._check_convergence = kwargs.pop('check_convergence', True)

super().__init__(**kwargs)

self._tau = None
self._sigma = None

Expand All @@ -232,22 +240,31 @@ def __init__(self, f, g, operator, tau=None, sigma=None, initial=None, gamma_g=N
self.set_gamma_g(gamma_g)
self.set_gamma_fconj(gamma_fconj)

self.set_up(f=f, g=g, operator=operator, tau=tau, sigma=sigma, initial=initial, **kwargs)
self.set_up(f=f, g=g, operator=operator, tau=tau, sigma=sigma, initial=initial)

@property
def tau(self):
"""The primal step-size """
return self._tau

@property
def sigma(self):
"""The dual step-size """
return self._sigma

@property
def theta(self):
"""The relaxation parameter for the over-relaxation of the primal variable """
return self._theta

@property
def gamma_g(self):
"""The strongly convex constant for the function g """
return self._gamma_g

@property
def gamma_fconj(self):
"""The strongly convex constant for the convex conjugate of the function f """
return self._gamma_fconj

def set_gamma_g(self, value):
Expand Down Expand Up @@ -290,7 +307,7 @@ def set_gamma_fconj(self, value):
else:
raise ValueError("Positive float is expected for the strongly convex constant of the convex conjugate of function f, {} is passed".format(value))

def set_up(self, f, g, operator, tau=None, sigma=None, initial=None, **kwargs):
def set_up(self, f, g, operator, tau=None, sigma=None, initial=None):
"""Initialisation of the algorithm

Parameters
Expand All @@ -301,14 +318,12 @@ def set_up(self, f, g, operator, tau=None, sigma=None, initial=None, **kwargs):
A convex function with a "simple" proximal.
operator : LinearOperator
A Linear Operator.
sigma : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default=None
sigma : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default is 1.0/norm(K) or 1.0/ (tau*norm(K)**2) if tau is provided
Step size for the dual problem.
tau : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default=None
tau : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default is 1.0/norm(K) or 1.0/ (sigma*norm(K)**2) if sigma is provided
Step size for the primal problem.
initial : DataContainer, optional, default=None
Initial point for the PDHG algorithm.
theta : Relaxation parameter, Number, default 1.0
"""
initial : DataContainer, optional, default is a DataContainer of zeros
Initial point for the PDHG algorithm. """
log.info("%s setting up", self.__class__.__name__)
# Triplet (f, g, K)
self.f = f
Expand All @@ -317,7 +332,7 @@ def set_up(self, f, g, operator, tau=None, sigma=None, initial=None, **kwargs):

self.set_step_sizes(sigma=sigma, tau=tau)

if kwargs.get('check_convergence', True):
if self._check_convergence:
self.check_convergence()

if initial is None:
Expand All @@ -330,8 +345,6 @@ def set_up(self, f, g, operator, tau=None, sigma=None, initial=None, **kwargs):
self.y = self.operator.range_geometry().allocate(0)
self.y_tmp = self.operator.range_geometry().allocate(0)

# relaxation parameter, default value is 1.0
self.theta = kwargs.get('theta',1.0)

if self.gamma_g is not None:
warnings.warn("Primal Acceleration of PDHG: The function g is assumed to be strongly convex with positive parameter `gamma_g`. You need to be sure that gamma_g = {} is the correct strongly convex constant for g. ".format(self.gamma_g))
Expand Down Expand Up @@ -447,14 +460,14 @@ def update_step_sizes(self):
"""
# Update sigma and tau based on the strong convexity of G
if self.gamma_g is not None:
self.theta = 1.0/ np.sqrt(1 + 2 * self.gamma_g * self.tau)
self._theta = 1.0/ np.sqrt(1 + 2 * self.gamma_g * self.tau)
self._tau *= self.theta
self._sigma /= self.theta

# Update sigma and tau based on the strong convexity of F
# Following operations are reversed due to symmetry, sigma --> tau, tau -->sigma
if self.gamma_fconj is not None:
self.theta = 1.0 / np.sqrt(1 + 2 * self.gamma_fconj * self.sigma)
self._theta = 1.0 / np.sqrt(1 + 2 * self.gamma_fconj * self.sigma)
self._sigma *= self.theta
self._tau /= self.theta

Expand Down
47 changes: 45 additions & 2 deletions Wrappers/Python/test/test_algorithms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@
if has_cvxpy:
import cvxpy


import warnings
class TestGD(CCPiTestClass):
def setUp(self):

Expand Down Expand Up @@ -700,6 +700,10 @@ def test_update(self):
beta= ((self.data - self.initial-alpha*(self.data-self.initial)).norm()**2)/4
self.assertNumpyArrayEqual(self.alg.p.as_array(), ((self.data - self.initial-alpha*(self.data-self.initial))+beta*(self.data-self.initial)).as_array())
class TestPDHG(CCPiTestClass):





def test_PDHG_Denoising(self):
# adapted from demo PDHG_TV_Color_Denoising.py in CIL-Demos repository
Expand Down Expand Up @@ -876,13 +880,32 @@ def test_PDHG_step_sizes(self):
with self.assertRaises(AttributeError):
pdhg = PDHG(f=f, g=g, operator=operator,
tau="tau")



# check warning message if condition is not satisfied
sigma = 4
sigma = 4/operator.norm()
tau = 1/3
with self.assertWarnsRegex(UserWarning, "Convergence criterion"):
pdhg = PDHG(f=f, g=g, operator=operator, tau=tau,
sigma=sigma)

# check no warning message if check convergence is false
sigma = 4/operator.norm()
tau = 1/3
with warnings.catch_warnings(record=True) as warnings_log:
pdhg = PDHG(f=f, g=g, operator=operator, tau=tau,
sigma=sigma, check_convergence=False)
self.assertEqual(warnings_log, [])

# check no warning message if condition is satisfied
sigma = 1/operator.norm()
tau = 1/3
with warnings.catch_warnings(record=True) as warnings_log:
pdhg = PDHG(f=f, g=g, operator=operator, tau=tau,
sigma=sigma)
self.assertEqual(warnings_log, [])


def test_PDHG_strongly_convex_gamma_g(self):
ig = ImageGeometry(3, 3)
Expand Down Expand Up @@ -970,7 +993,27 @@ def test_PDHG_strongly_convex_both_fconj_and_g(self):
except ValueError as err:
log.info(str(err))

def test_pdhg_theta(self):
ig = ImageGeometry(3, 3)
data = ig.allocate('random')

f = L2NormSquared(b=data)
g = L2NormSquared()
operator = IdentityOperator(ig)

pdhg = PDHG(f=f, g=g, operator=operator)
self.assertEqual(pdhg.theta, 1.0)

pdhg = PDHG(f=f, g=g, operator=operator,theta=0.5)
self.assertEqual(pdhg.theta, 0.5)

with self.assertRaises(ValueError):
PDHG( f=f, g=g, operator=operator, theta=-0.5)

with self.assertRaises(ValueError):
PDHG( f=f, g=g, operator=operator, theta=5)


class TestSIRT(CCPiTestClass):

def setUp(self):
Expand Down
Loading