From 5d9328928f5a6df58fe1d1c3817dd26f192437ae Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Wed, 28 Jan 2026 11:42:09 +0000 Subject: [PATCH 1/7] add ProxSKip algo --- .../cil/optimisation/algorithms/ProxSkip.py | 104 ++++++++++++++++++ .../cil/optimisation/algorithms/__init__.py | 3 +- 2 files changed, 106 insertions(+), 1 deletion(-) create mode 100644 Wrappers/Python/cil/optimisation/algorithms/ProxSkip.py diff --git a/Wrappers/Python/cil/optimisation/algorithms/ProxSkip.py b/Wrappers/Python/cil/optimisation/algorithms/ProxSkip.py new file mode 100644 index 0000000000..27be18f011 --- /dev/null +++ b/Wrappers/Python/cil/optimisation/algorithms/ProxSkip.py @@ -0,0 +1,104 @@ +from cil.optimisation.algorithms import Algorithm +import numpy as np +import logging + + +class ProxSkip(Algorithm): + + + r"""Proximal Skip (ProxSkip) algorithm, see "ProxSkip: Yes! Local Gradient Steps Provably Lead to Communication Acceleration! Finally!†" + + Parameters + ---------- + + initial : DataContainer + Initial point for the ProxSkip algorithm. + f : Function + A smooth function with Lipschitz continuous gradient. + g : Function + A convex function with a "simple" proximal. + prob : positive :obj:`float` + Probability to skip the proximal step. If :code:`prob=1`, proximal step is used in every iteration. + step_size : positive :obj:`float` + Step size for the ProxSkip algorithm and is equal to 1./L where L is the Lipschitz constant for the gradient of f. + + """ + + + def __init__(self, initial, f, g, step_size, prob, seed=None, **kwargs): + """ Set up of the algorithm + """ + + super(ProxSkip, self).__init__(**kwargs) + + self.f = f # smooth function + self.g = g # proximable + self.step_size = step_size + self.prob = prob + self.rng = np.random.default_rng(seed) + self.set_up(initial, f, g, step_size, prob, **kwargs) + self.thetas = [] + + + def set_up(self, initial, f, g, step_size, prob, **kwargs): + + logging.info("{} setting up".format(self.__class__.__name__, )) + + self.initial = initial[0] + + self.x = initial.copy() + self.xhat_new = initial.copy() + self.x_new = initial.copy() + self.ht = initial.copy() + + self.configured = True + + logging.info("{} configured".format(self.__class__.__name__, )) + + + def update(self): + r""" Performs a single iteration of the ProxSkip algorithm + """ + + self.f.gradient(self.x, out=self.xhat_new) + self.xhat_new -= self.ht + self.x.sapyb(1., self.xhat_new, -self.step_size, out=self.xhat_new) + + theta = self.rng.random() < self.prob + # convention: use proximal in the first iteration + if self.iteration==0: + theta = 1 + self.thetas.append(theta) + + if theta==1: + # Proximal step is used + self.g.proximal(self.xhat_new - (self.step_size/self.prob)*self.ht, self.step_size/self.prob, out=self.x_new) + self.ht.sapyb(1., (self.x_new - self.xhat_new), (self.prob/self.step_size), out=self.ht) + else: + self.x_new.fill(self.xhat_new) + + def _update_previous_solution(self): + """ Swaps the references to current and previous solution based on the :func:`~Algorithm.update_previous_solution` of the base class :class:`Algorithm`. + """ + tmp = self.x_new + self.x = self.x_new + self.x = tmp + + def get_output(self): + " Returns the current solution. " + return self.x + + + def update_objective(self): + + """ Updates the objective + + .. math:: f(x) + g(x) + + """ + + fun_g = self.g(self.x) + fun_f = self.f(self.x) + p1 = fun_f + fun_g + self.loss.append( p1 ) + diff --git a/Wrappers/Python/cil/optimisation/algorithms/__init__.py b/Wrappers/Python/cil/optimisation/algorithms/__init__.py index 0ca481beaa..9237b2b91e 100644 --- a/Wrappers/Python/cil/optimisation/algorithms/__init__.py +++ b/Wrappers/Python/cil/optimisation/algorithms/__init__.py @@ -22,10 +22,11 @@ from .GD import GD from .FISTA import FISTA from .FISTA import ISTA +from .ProxSkip import ProxSkip from .FISTA import ISTA as PGD from .APGD import APGD from .PDHG import PDHG from .ADMM import LADMM from .SPDHG import SPDHG from .PD3O import PD3O -from .LSQR import LSQR \ No newline at end of file +from .LSQR import LSQR From 4dbaea7cfd98ea043093805244fe5d2fab584d7a Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Wed, 28 Jan 2026 14:12:08 +0000 Subject: [PATCH 2/7] initial tests, need more, CVXpy not working on mac --- .../cil/optimisation/algorithms/ProxSkip.py | 5 +- Wrappers/Python/test/test_algorithms.py | 49 +++++++++++++++++-- 2 files changed, 50 insertions(+), 4 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/algorithms/ProxSkip.py b/Wrappers/Python/cil/optimisation/algorithms/ProxSkip.py index 27be18f011..414b2fb3d1 100644 --- a/Wrappers/Python/cil/optimisation/algorithms/ProxSkip.py +++ b/Wrappers/Python/cil/optimisation/algorithms/ProxSkip.py @@ -38,13 +38,16 @@ def __init__(self, initial, f, g, step_size, prob, seed=None, **kwargs): self.rng = np.random.default_rng(seed) self.set_up(initial, f, g, step_size, prob, **kwargs) self.thetas = [] + + ## TODO raise warning if p=1, better to use non-skipped version + # because update of control variate is not used. def set_up(self, initial, f, g, step_size, prob, **kwargs): logging.info("{} setting up".format(self.__class__.__name__, )) - self.initial = initial[0] + self.initial = initial self.x = initial.copy() self.xhat_new = initial.copy() diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py index 7908debf66..4a9e0d4474 100644 --- a/Wrappers/Python/test/test_algorithms.py +++ b/Wrappers/Python/test/test_algorithms.py @@ -38,7 +38,7 @@ from cil.optimisation.functions import MixedL21Norm, BlockFunction, L1Norm, KullbackLeibler, IndicatorBox, LeastSquares, ZeroFunction, L2NormSquared, OperatorCompositionFunction, TotalVariation, SGFunction, SVRGFunction, SAGAFunction, SAGFunction, LSVRGFunction, ScaledFunction -from cil.optimisation.algorithms import Algorithm, GD, CGLS, SIRT, FISTA, ISTA, SPDHG, PDHG, LADMM, PD3O, PGD, APGD , LSQR +from cil.optimisation.algorithms import Algorithm, GD, CGLS, SIRT, FISTA, ISTA, SPDHG, PDHG, LADMM, PD3O, PGD, APGD , LSQR, ProxSkip from scipy.optimize import minimize, rosen @@ -340,8 +340,7 @@ def test_provable_convergence(self): with self.assertRaises(NotImplementedError): alg.is_provably_convergent() - - + class TestFISTA(CCPiTestClass): @@ -533,6 +532,50 @@ def get_step_size(self, algorithm): self.assertEqual(alg.step_size, 0.99/2) self.assertEqual(alg.step_size, 0.99/2) + +class testProxSkip(CCPiTestClass): + + def setUp(self): + + np.random.seed(10) + n = 50 + m = 500 + + A = np.random.uniform(0, 1, (m, n)).astype('float32') + b = (A.dot(np.random.randn(n)) + 0.1 * + np.random.randn(m)).astype('float32') + + self.Aop = MatrixOperator(A) + self.bop = VectorData(b) + + self.f = LeastSquares(self.Aop, b=self.bop, c=0.5) + self.g = ZeroFunction() + self.h = L1Norm() + + self.ig = self.Aop.domain + + self.initial = self.ig.allocate() + + def tearDown(self): + pass + + @unittest.skipUnless(has_cvxpy, "CVXpy not installed") + def test_with_cvxpy(self): + + ista = ProxSkip(initial=self.initial, f=self.f, + g=self.g, prob = 0.1) + ista.run(2000, verbose=0) + + u_cvxpy = cvxpy.Variable(self.ig.shape[0]) + objective = cvxpy.Minimize( + 0.5 * cvxpy.sum_squares(self.Aop.A @ u_cvxpy - self.bop.array)) + p = cvxpy.Problem(objective) + p.solve(verbose=True, solver=cvxpy.SCS, eps=1e-4) + + np.testing.assert_allclose(p.value, ista.objective[-1], atol=1e-3) + np.testing.assert_allclose( + u_cvxpy.value, ista.solution.array, atol=1e-3) + class testISTA(CCPiTestClass): def setUp(self): From 0df1ff6099f120ec86088ce8cd5b7053ce4997a5 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Wed, 28 Jan 2026 14:30:32 +0000 Subject: [PATCH 3/7] ISTA vs ProxSkip, more tests needed, test signature --- Wrappers/Python/test/test_algorithms.py | 32 ++++++++++++------------- 1 file changed, 15 insertions(+), 17 deletions(-) diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py index 4a9e0d4474..09ae84f552 100644 --- a/Wrappers/Python/test/test_algorithms.py +++ b/Wrappers/Python/test/test_algorithms.py @@ -57,8 +57,8 @@ log = logging.getLogger(__name__) -if has_cvxpy: - import cvxpy +# if has_cvxpy: + # import cvxpy from unittest.mock import MagicMock, patch class TestGD(CCPiTestClass): @@ -533,7 +533,7 @@ def get_step_size(self, algorithm): self.assertEqual(alg.step_size, 0.99/2) -class testProxSkip(CCPiTestClass): +class TestProxSkip(CCPiTestClass): def setUp(self): @@ -550,7 +550,7 @@ def setUp(self): self.f = LeastSquares(self.Aop, b=self.bop, c=0.5) self.g = ZeroFunction() - self.h = L1Norm() + self.h = 0.5 * L1Norm() self.ig = self.Aop.domain @@ -559,22 +559,20 @@ def setUp(self): def tearDown(self): pass - @unittest.skipUnless(has_cvxpy, "CVXpy not installed") - def test_with_cvxpy(self): + def test_ista_vs_proxskip(self): - ista = ProxSkip(initial=self.initial, f=self.f, - g=self.g, prob = 0.1) - ista.run(2000, verbose=0) - - u_cvxpy = cvxpy.Variable(self.ig.shape[0]) - objective = cvxpy.Minimize( - 0.5 * cvxpy.sum_squares(self.Aop.A @ u_cvxpy - self.bop.array)) - p = cvxpy.Problem(objective) - p.solve(verbose=True, solver=cvxpy.SCS, eps=1e-4) + prox = ProxSkip(initial=self.initial, f=self.f, + g=self.h, step_size = 1.99/self.f.L, prob = 0.1) + prox.run(2000, verbose=0) - np.testing.assert_allclose(p.value, ista.objective[-1], atol=1e-3) + ista = ISTA(initial=self.initial, f=self.f, + g=self.h, step_size = 1.99/self.f.L) + ista.run(1000, verbose=0) + + np.testing.assert_allclose(ista.objective[-1], prox.objective[-1], atol=1e-3) np.testing.assert_allclose( - u_cvxpy.value, ista.solution.array, atol=1e-3) + prox.solution.array, ista.solution.array, atol=1e-3) + class testISTA(CCPiTestClass): From b23f4050377e04e10d3bdfd3c07bd713749dfde5 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 29 Jan 2026 10:12:15 +0000 Subject: [PATCH 4/7] more tests, problem with random seed --- .../cil/optimisation/algorithms/ProxSkip.py | 23 +++--- Wrappers/Python/test/test_algorithms.py | 79 ++++++++++++++++++- 2 files changed, 89 insertions(+), 13 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/algorithms/ProxSkip.py b/Wrappers/Python/cil/optimisation/algorithms/ProxSkip.py index 414b2fb3d1..43da3b14b4 100644 --- a/Wrappers/Python/cil/optimisation/algorithms/ProxSkip.py +++ b/Wrappers/Python/cil/optimisation/algorithms/ProxSkip.py @@ -1,6 +1,7 @@ from cil.optimisation.algorithms import Algorithm import numpy as np import logging +from warnings import warn class ProxSkip(Algorithm): @@ -20,7 +21,7 @@ class ProxSkip(Algorithm): prob : positive :obj:`float` Probability to skip the proximal step. If :code:`prob=1`, proximal step is used in every iteration. step_size : positive :obj:`float` - Step size for the ProxSkip algorithm and is equal to 1./L where L is the Lipschitz constant for the gradient of f. + Step size for the ProxSkip algorithm. It is equal to 1./L for strongly convex f and 2./L for convex f, where L is the Lipschitz constant for the gradient of f. """ @@ -35,20 +36,23 @@ def __init__(self, initial, f, g, step_size, prob, seed=None, **kwargs): self.g = g # proximable self.step_size = step_size self.prob = prob - self.rng = np.random.default_rng(seed) - self.set_up(initial, f, g, step_size, prob, **kwargs) + self.rng = np.random.default_rng(seed=seed) self.thetas = [] - ## TODO raise warning if p=1, better to use non-skipped version - # because update of control variate is not used. - + if self.prob<=0: + raise ValueError("Need a positive probability") + if self.prob==1: + raise warn("If p=1, ProxSkip is equivalent to ISTA/PGD. Please use ISTA/PGD to avoid computing updates of the control variate that is not used.") + + self.set_up(initial, f, g, step_size, prob, **kwargs) + def set_up(self, initial, f, g, step_size, prob, **kwargs): logging.info("{} setting up".format(self.__class__.__name__, )) + ## TODO better to use different initials for x and h. self.initial = initial - self.x = initial.copy() self.xhat_new = initial.copy() self.x_new = initial.copy() @@ -70,10 +74,10 @@ def update(self): theta = self.rng.random() < self.prob # convention: use proximal in the first iteration if self.iteration==0: - theta = 1 + theta = True self.thetas.append(theta) - if theta==1: + if theta: # Proximal step is used self.g.proximal(self.xhat_new - (self.step_size/self.prob)*self.ht, self.step_size/self.prob, out=self.x_new) self.ht.sapyb(1., (self.x_new - self.xhat_new), (self.prob/self.step_size), out=self.ht) @@ -104,4 +108,5 @@ def update_objective(self): fun_f = self.f(self.x) p1 = fun_f + fun_g self.loss.append( p1 ) + diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py index 09ae84f552..9a27335dc3 100644 --- a/Wrappers/Python/test/test_algorithms.py +++ b/Wrappers/Python/test/test_algorithms.py @@ -549,8 +549,8 @@ def setUp(self): self.bop = VectorData(b) self.f = LeastSquares(self.Aop, b=self.bop, c=0.5) - self.g = ZeroFunction() - self.h = 0.5 * L1Norm() + self.g = 0.5 * L1Norm() + self.step_size = 1.99/self.f.L self.ig = self.Aop.domain @@ -559,14 +559,85 @@ def setUp(self): def tearDown(self): pass + def test_signature(self): + + # check required arguments (initial, f, g, step size, and prob) + with np.testing.assert_raises(TypeError): + proxskip = ProxSkip(initial = self.initial, f=self.f, g=self.g, step_size=self.step_size) + + # test neg prob + with np.testing.assert_raises(ValueError): + proxskip = ProxSkip(initial = self.initial, f=self.f, g=self.g, step_size=self.step_size, prob=-0.1) + + # zero prob + with np.testing.assert_raises(ValueError): + proxskip = ProxSkip(initial = self.initial, f=self.f, g=self.g, step_size=self.step_size, prob=0.) + + def test_coin_flip(self): + + seed = 10 + num_it = 100 + prob = 0.3 + proxskip1 = ProxSkip(initial = self.initial, f=self.f, g=self.g, + step_size=self.step_size, prob=prob, seed=seed) + proxskip1.run(num_it, verbose=0) + + thetas1 = [] + rng = np.random.default_rng(seed) + for i in range(num_it-1): + if i==0: + thetas1.append(True) + theta = proxskip1.rng.random() < prob + thetas1.append(theta) + + # thetas2 = [] + # rng = np.random.default_rng(seed) + # for i in range(num_it-1): + # if i==0: + # thetas2.append(True) + # theta = rng.random() < prob + # thetas2.append(theta) + + # np.testing.assert_allclose(thetas1, thetas2) + np.testing.assert_allclose(proxskip1.thetas, thetas1) + + + + + + def test_seeds(self): + + # same seeds + proxskip1 = ProxSkip(initial = self.initial, f=self.f, g=self.g, step_size=self.step_size, prob=0.1, seed=10) + proxskip1.run(100, verbose=0) + + proxskip2 = ProxSkip(initial = self.initial, f=self.f, g=self.g, step_size=self.step_size, prob=0.1, seed=10) + proxskip2.run(100, verbose=0) + + np.testing.assert_allclose(proxskip2.thetas, proxskip1.thetas) + + # different seeds + proxskip2 = ProxSkip(initial = self.initial, f=self.f, g=self.g, step_size=self.step_size, + prob=0.1, seed=20) + proxskip2.run(100, verbose=0) + + assert not np.array_equal(proxskip2.thetas, proxskip1.thetas) + + + + + + + + def test_ista_vs_proxskip(self): prox = ProxSkip(initial=self.initial, f=self.f, - g=self.h, step_size = 1.99/self.f.L, prob = 0.1) + g=self.h, step_size = self.step_size, prob = 0.1) prox.run(2000, verbose=0) ista = ISTA(initial=self.initial, f=self.f, - g=self.h, step_size = 1.99/self.f.L) + g=self.h, step_size = self.step_size) ista.run(1000, verbose=0) np.testing.assert_allclose(ista.objective[-1], prox.objective[-1], atol=1e-3) From aa270f415609c031e95ce356aad77c5cf7e295b5 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 29 Jan 2026 10:23:13 +0000 Subject: [PATCH 5/7] fix coin flip test --- Wrappers/Python/test/test_algorithms.py | 39 +++++++------------------ 1 file changed, 11 insertions(+), 28 deletions(-) diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py index 9a27335dc3..65d337e59a 100644 --- a/Wrappers/Python/test/test_algorithms.py +++ b/Wrappers/Python/test/test_algorithms.py @@ -577,34 +577,23 @@ def test_coin_flip(self): seed = 10 num_it = 100 - prob = 0.3 - proxskip1 = ProxSkip(initial = self.initial, f=self.f, g=self.g, - step_size=self.step_size, prob=prob, seed=seed) - proxskip1.run(num_it, verbose=0) + prob = 0.2 + + proxskip1 = ProxSkip(initial=self.initial, f=self.f, g=self.g, + step_size=self.step_size, prob=prob, seed=seed) + proxskip1.run(num_it, verbose=0) - thetas1 = [] rng = np.random.default_rng(seed) - for i in range(num_it-1): - if i==0: - thetas1.append(True) - theta = proxskip1.rng.random() < prob + + thetas1 = [] + for k in range(num_it): + tmp = rng.random() < prob + theta = True if k == 0 else tmp thetas1.append(theta) - # thetas2 = [] - # rng = np.random.default_rng(seed) - # for i in range(num_it-1): - # if i==0: - # thetas2.append(True) - # theta = rng.random() < prob - # thetas2.append(theta) + assert np.array_equal(proxskip1.thetas, thetas1) - # np.testing.assert_allclose(thetas1, thetas2) - np.testing.assert_allclose(proxskip1.thetas, thetas1) - - - - def test_seeds(self): # same seeds @@ -624,12 +613,6 @@ def test_seeds(self): assert not np.array_equal(proxskip2.thetas, proxskip1.thetas) - - - - - - def test_ista_vs_proxskip(self): prox = ProxSkip(initial=self.initial, f=self.f, From 259b3244540dd34291133956a0f237bfab134332 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 29 Jan 2026 10:31:35 +0000 Subject: [PATCH 6/7] fix input function --- Wrappers/Python/test/test_algorithms.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py index 65d337e59a..7f3cf7d6de 100644 --- a/Wrappers/Python/test/test_algorithms.py +++ b/Wrappers/Python/test/test_algorithms.py @@ -616,11 +616,11 @@ def test_seeds(self): def test_ista_vs_proxskip(self): prox = ProxSkip(initial=self.initial, f=self.f, - g=self.h, step_size = self.step_size, prob = 0.1) + g=self.g, step_size = self.step_size, prob = 0.1) prox.run(2000, verbose=0) ista = ISTA(initial=self.initial, f=self.f, - g=self.h, step_size = self.step_size) + g=self.g, step_size = self.step_size) ista.run(1000, verbose=0) np.testing.assert_allclose(ista.objective[-1], prox.objective[-1], atol=1e-3) From d95f7bdc28e7cd4c6058825965aee127449ed817 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 29 Jan 2026 11:03:55 +0000 Subject: [PATCH 7/7] local prob with cvxpy, uncomment --- Wrappers/Python/test/test_algorithms.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py index 7f3cf7d6de..34d24660df 100644 --- a/Wrappers/Python/test/test_algorithms.py +++ b/Wrappers/Python/test/test_algorithms.py @@ -57,8 +57,8 @@ log = logging.getLogger(__name__) -# if has_cvxpy: - # import cvxpy +if has_cvxpy: + import cvxpy from unittest.mock import MagicMock, patch class TestGD(CCPiTestClass):