From 11c260dcb6c8fcf742f51c0a1faa2e21f757abb7 Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Tue, 15 Nov 2022 15:24:33 +0100 Subject: [PATCH 01/13] add implementation for admm and fast admm. Based on Goldstein2014 --- modopt/opt/algorithms/admm.py | 291 ++++++++++++++++++++++++++++++++++ 1 file changed, 291 insertions(+) create mode 100644 modopt/opt/algorithms/admm.py diff --git a/modopt/opt/algorithms/admm.py b/modopt/opt/algorithms/admm.py new file mode 100644 index 00000000..d45c1d38 --- /dev/null +++ b/modopt/opt/algorithms/admm.py @@ -0,0 +1,291 @@ +"""ADMM Algorithms.""" +import numpy as np + +from modopt.opt.algorithms.base import SetUp +from modopt.opt.cost import costObj + + +class ADMM(SetUp): + r"""Fast ADMM Optimisation Algorihm. + + This class implement the ADMM algorithm (Algorithm 1 from :cite:`Goldstein2014`) + + Parameters + ---------- + A : OperatorBase + Linear operator for u + B : OperatorBase + Linear operator for v + b : array_like + Constraint vector + optimizers: 2-tuple of functions + Solvers for the u and v update, takes init_value and obs_value as argument. + and returns an estimate for: + .. math:: u_{k+1} = \argmin H(u) + \frac{\tau}{2}\|A u - y\|^2 + .. math:: v_{k+1} = \argmin G(v) + \frac{\tau}{2}\|Bv - y \|^2 + cost_funcs = 2-tuple of function + Compute the values of H and G + rho : float , optional + regularisation coupling variable default is ``1.0`` + eta : float, optional + Restart threshold, default is ``0.999`` + + Notes + ----- + The algorithm solve the problem: + + .. math:: u, v = \arg\min H(u) + G(v) + \frac{\tau}{2} \|Au + Bv - b \|_2^2 + + with the following augmented lagrangian: + + .. math :: \mathcal{L}_{\tau}(u,v, \lambda) = H(u) + G(v) + +\langle\lambda |Au + Bv -b \rangle + \frac\tau2 \| Au + Bv -b \|^{2} + + To allow easy iterative solving, the change of variable :math:`\mu=\lambda/\tau` + is used. Hence, the lagrangian of interest is: + + .. math :: \tilde{\mathcal{L}}_{\tau}(u,v, \mu) = H(u) + G(v) + + \frac\tau2 \left(\|\mu + Au +Bv - b\|^2 - \|\mu\|^2\right) + + See Also + -------- + SetUp: parent class + """ + + def __init__( + self, + u, + v, + mu, + A, + B, + b, + optimizers, + tau=1, + max_iter2=5, + cost_func=None, + **kwargs + ): + super().__init__(**kwargs) + self.A = A + self.B = B + self.b = b + self._opti_H = optimizers[0] + self._opti_G = optimizers[1] + self._tau = tau + + self._cost_func = costObj() + # patching to get the full cost + self._cost_func._calc_cost = lambda u, v: ( + cost_func[0](u) + + cost_func[1](v) + + self.xp.linalg.norm(A.op(u) + B.op(v) - b) + ) + + # init iteration variables. + self._u_old = self.xp.copy(u) + self._u_new = self.xp.copy(u) + self._v_old = self.xp.copy(v) + self._v_new = self.xp.copy(v) + self._mu_new = self.xp.copy(mu) + self._mu_old = self.xp.copy(mu) + + def _update(self): + self._u_new = self._opti_H( + init=self._u_old, + obs=self.B.op(self._v_old) + self._u_old - self.b, + ) + uA_new = self.A.op(self._u_new) + self._v_new = self.solver2( + init=self._v_old, + obs=uA_new + self._u_old - self.c, + ) + + self._mu_new = self._mu_old + (uA_new + self.B.op(self._v_new) - self.b) + + # update cycle + self._u_old = self.xp.copy(self._u_new) + self._v_old = self.xp.copy(self._v_new) + self._mu_old = self.xp.copy(self._mu_new) + + # Test cost function for convergence. + if self._cost_func: + self.converge = self.any_convergence_flag() or self._cost_func.get_cost( + self._u_new, self.v_new + ) + + def iterate(self, max_iter=150): + """Iterate. + + This method calls update until either convergence criteria is met or + the maximum number of iterations is reached. + + Parameters + ---------- + max_iter : int, optional + Maximum number of iterations (default is ``150``) + """ + self._run_alg(max_iter) + + # retrieve metrics results + self.retrieve_outputs() + # rename outputs as attributes + self.u_final = self._u_new + self.v_final = self._v_new + + def get_notify_observers_kwargs(self): + """Notify observers. + + Return the mapping between the metrics call and the iterated + variables. + + Returns + ------- + dict + The mapping between the iterated variables + """ + return { + "u_new": self._u_new, + "v_new": self._v_new, + "idx": self.idx, + } + + def retrieve_outputs(self): + """Retrieve outputs. + + Declare the outputs of the algorithms as attributes: x_final, + y_final, metrics. + """ + metrics = {} + for obs in self._observers["cv_metrics"]: + metrics[obs.name] = obs.retrieve_metrics() + self.metrics = metrics + + +class FastADMM(ADMM): + r"""Fast ADMM Optimisation Algorihm. + + This class implement the fast ADMM algorithm (Algorithm 8 from :cite:`Goldstein2014`) + + Parameters + ---------- + A : OperatorBase + Linear operator for u + B : OperatorBase + Linear operator for v + b : array_like + Constraint vector + solver1 : function + Solver for the x update, takes init_value and obs_value as argument. + ie, return an estimate for: + .. math:: u_{k+1} = \argmin H(u) + \frac{\tau}{2}\|A u - y\|^2 + solver2 : function + Solver for the z update, takes init_value and obs_value as argument. + ie return an estimate for: + .. math:: v_{k+1} = \argmin G(v) + \frac{\tau}{2}\|Bv - y \|^2 + rho : float , optional + regularisation coupling variable default is ``1.0`` + eta : float, optional + Restart threshold, default is ``0.999`` + + Notes + ----- + The algorithm solve the problem: + + .. math:: u, v = \arg\min H(u) + G(v) + \frac{\tau}{2} \|Au + Bv - b \|_2^2 + + with the following augmented lagrangian: + + .. math :: \mathcal{L}_{\tau}(u,v, \lambda) = H(u) + G(v) + +\langle\lambda |Au + Bv -b \rangle + \frac\tau2 \| Au + Bv -b \|^{2} + + To allow easy iterative solving, the change of variable :math:`\mu=\lambda/\tau` + is used. Hence, the lagrangian of interest is: + + .. math :: \tilde{\mathcal{L}}_{\tau}(u,v, \mu) = H(u) + G(v) + + \frac\tau2 \left(\|\mu + Au +Bv - b\|^2 - \|\mu\|^2\right) + + See Also + -------- + SetUp: parent class + """ + + def __init__( + self, + u, + v, + mu, + A, + B, + b, + opti_H, + opti_G, + alpha=1, + eta=0.999, + tau=1, + opti_H_kwargs=None, + opti_G_kwargs=None, + cost=None, + **kwargs + ): + super().__init__( + u=u, + v=b, + mu=mu, + A=A, + B=B, + b=b, + opti_H=opti_H, + opti_G=opti_G, + opti_H_kwargs=opti_H_kwargs, + opti_G_kwargs=opti_G, + cost=None, + **kwargs + ) + self._c_old = np.inf + self._c_new = 0.0 + self._eta = eta + self._alpha_old = alpha + self._alpha_new = alpha + self._v_hat = self.xp.copy(self._v_new) + self._mu_hat = self.xp.copy(self._mu_new) + + def _update(self): + # Classical ADMM steps + self._u_new = self._opti_H( + init=self._u_old, + obs=self.B.op(self._v_hat) + self._u_old - self.b, + ) + uA_new = self.A.op(self._u_new) + self._v_new = self.solver2( + init=self._v_hat, + obs=uA_new + self._u_hat - self.c, + ) + + self._mu_new = self._mu_hat + (uA_new + self.B.op(self._v_new) - self.b) + + # restarting condition + self._c_new = self.xp.linalg.norm(self._mu_new - self._mu_hat) + self._c_new += self._tau * self.xp.linalg.norm( + self.B.op(self._v_new - self._v_hat) + ) + if self._c_new < self._eta * self._c_old: + self._alpha_new = 1 + np.sqrt(1 + 4 * self._alpha_old**2) + update_factor = (self._alpha_new - 1) / self._alpha_old + self._v_hat = self._v_new + (self._v_new - self._v_old) * update_factor + self._mu_hat = self._mu_new + (self._mu_new - self._mu_old) * update_factor + else: + # reboot to old iteration + self._alpha_new = 1 + self._v_hat = self._v_old + self._mu_hat = self._mu_old + self._c_new = self._c_old / self._eta + + self.xp.copyto(self._u_old, self._u_new) + self.xp.copyto(self._v_old, self._v_new) + self.xp.copyto(self._mu_old, self._mu_new) + # Test cost function for convergence. + if self._cost_func: + self.converge = self.any_convergence_flag() or self._cost_func.get_cost( + self._u_new + ) From 3e0bef3217646bab2a486f3444e42ee73e8af496 Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Tue, 15 Nov 2022 15:31:21 +0100 Subject: [PATCH 02/13] add Goldstein ref. --- docs/source/refs.bib | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/docs/source/refs.bib b/docs/source/refs.bib index d8365e71..7782ca52 100644 --- a/docs/source/refs.bib +++ b/docs/source/refs.bib @@ -207,3 +207,15 @@ @article{zou2005 journal = {Journal of the Royal Statistical Society Series B}, doi = {10.1111/j.1467-9868.2005.00527.x} } + +@article{Goldstein2014, + author={Goldstein, Tom and O’Donoghue, Brendan and Setzer, Simon and Baraniuk, Richard}, + year={2014}, + month={Jan}, + pages={1588–1623}, + title={Fast Alternating Direction Optimization Methods}, + journal={SIAM Journal on Imaging Sciences}, + volume={7}, + ISSN={1936-4954}, + doi={10/gdwr49}, +} From 83459ada2082e9f2f4ba7e3a9943aeb349b88a66 Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Tue, 15 Nov 2022 16:08:45 +0100 Subject: [PATCH 03/13] WPS compliance. --- modopt/opt/algorithms/admm.py | 56 +++++++++++++++++------------------ 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/modopt/opt/algorithms/admm.py b/modopt/opt/algorithms/admm.py index d45c1d38..1830a846 100644 --- a/modopt/opt/algorithms/admm.py +++ b/modopt/opt/algorithms/admm.py @@ -8,7 +8,8 @@ class ADMM(SetUp): r"""Fast ADMM Optimisation Algorihm. - This class implement the ADMM algorithm (Algorithm 1 from :cite:`Goldstein2014`) + This class implement the ADMM algorithm + (Algorithm 1 from :cite:`Goldstein2014`) Parameters ---------- @@ -19,8 +20,8 @@ class ADMM(SetUp): b : array_like Constraint vector optimizers: 2-tuple of functions - Solvers for the u and v update, takes init_value and obs_value as argument. - and returns an estimate for: + Solvers for the u and v update, takes init_value and obs_value as + argument. each element returns an estimate for: .. math:: u_{k+1} = \argmin H(u) + \frac{\tau}{2}\|A u - y\|^2 .. math:: v_{k+1} = \argmin G(v) + \frac{\tau}{2}\|Bv - y \|^2 cost_funcs = 2-tuple of function @@ -34,15 +35,15 @@ class ADMM(SetUp): ----- The algorithm solve the problem: - .. math:: u, v = \arg\min H(u) + G(v) + \frac{\tau}{2} \|Au + Bv - b \|_2^2 + .. math:: u, v = \arg\min H(u) + G(v) + \frac\tau2 \|Au + Bv - b \|_2^2 with the following augmented lagrangian: .. math :: \mathcal{L}_{\tau}(u,v, \lambda) = H(u) + G(v) - +\langle\lambda |Au + Bv -b \rangle + \frac\tau2 \| Au + Bv -b \|^{2} + +\langle\lambda |Au + Bv -b \rangle + \frac\tau2 \| Au + Bv -b \|^2 - To allow easy iterative solving, the change of variable :math:`\mu=\lambda/\tau` - is used. Hence, the lagrangian of interest is: + To allow easy iterative solving, the change of variable + :math:`\mu=\lambda/\tau` is used. Hence, the lagrangian of interest is: .. math :: \tilde{\mathcal{L}}_{\tau}(u,v, \mu) = H(u) + G(v) + \frac\tau2 \left(\|\mu + Au +Bv - b\|^2 - \|\mu\|^2\right) @@ -95,13 +96,13 @@ def _update(self): init=self._u_old, obs=self.B.op(self._v_old) + self._u_old - self.b, ) - uA_new = self.A.op(self._u_new) + tmp = self.A.op(self._u_new) self._v_new = self.solver2( init=self._v_old, - obs=uA_new + self._u_old - self.c, + obs=tmp + self._u_old - self.c, ) - self._mu_new = self._mu_old + (uA_new + self.B.op(self._v_new) - self.b) + self._mu_new = self._mu_old + (tmp + self.B.op(self._v_new) - self.b) # update cycle self._u_old = self.xp.copy(self._u_new) @@ -110,9 +111,8 @@ def _update(self): # Test cost function for convergence. if self._cost_func: - self.converge = self.any_convergence_flag() or self._cost_func.get_cost( - self._u_new, self.v_new - ) + self.converge = self.any_convergence_flag() + self.converge |= self._cost_func.get_cost(self._u_new, self.v_new) def iterate(self, max_iter=150): """Iterate. @@ -165,7 +165,8 @@ def retrieve_outputs(self): class FastADMM(ADMM): r"""Fast ADMM Optimisation Algorihm. - This class implement the fast ADMM algorithm (Algorithm 8 from :cite:`Goldstein2014`) + This class implement the fast ADMM algorithm + (Algorithm 8 from :cite:`Goldstein2014`) Parameters ---------- @@ -192,15 +193,15 @@ class FastADMM(ADMM): ----- The algorithm solve the problem: - .. math:: u, v = \arg\min H(u) + G(v) + \frac{\tau}{2} \|Au + Bv - b \|_2^2 + .. math:: u, v = \arg\min H(u) + G(v) + \frac\tau2 \|Au + Bv - b \|_2^2 with the following augmented lagrangian: - .. math :: \mathcal{L}_{\tau}(u,v, \lambda) = H(u) + G(v) - +\langle\lambda |Au + Bv -b \rangle + \frac\tau2 \| Au + Bv -b \|^{2} + .. math:: \mathcal{L}_{\tau}(u,v, \lambda) = H(u) + G(v) + +\langle\lambda |Au + Bv -b \rangle + \frac\tau2 \| Au + Bv -b \|^2 - To allow easy iterative solving, the change of variable :math:`\mu=\lambda/\tau` - is used. Hence, the lagrangian of interest is: + To allow easy iterative solving, the change of variable + :math:`\mu=\lambda/\tau` is used. Hence, the lagrangian of interest is: .. math :: \tilde{\mathcal{L}}_{\tau}(u,v, \mu) = H(u) + G(v) + \frac\tau2 \left(\|\mu + Au +Bv - b\|^2 - \|\mu\|^2\right) @@ -256,13 +257,13 @@ def _update(self): init=self._u_old, obs=self.B.op(self._v_hat) + self._u_old - self.b, ) - uA_new = self.A.op(self._u_new) + tmp = self.A.op(self._u_new) self._v_new = self.solver2( init=self._v_hat, - obs=uA_new + self._u_hat - self.c, + obs=tmp + self._u_hat - self.c, ) - self._mu_new = self._mu_hat + (uA_new + self.B.op(self._v_new) - self.b) + self._mu_new = self._mu_hat + (tmp + self.B.op(self._v_new) - self.b) # restarting condition self._c_new = self.xp.linalg.norm(self._mu_new - self._mu_hat) @@ -271,9 +272,9 @@ def _update(self): ) if self._c_new < self._eta * self._c_old: self._alpha_new = 1 + np.sqrt(1 + 4 * self._alpha_old**2) - update_factor = (self._alpha_new - 1) / self._alpha_old - self._v_hat = self._v_new + (self._v_new - self._v_old) * update_factor - self._mu_hat = self._mu_new + (self._mu_new - self._mu_old) * update_factor + beta = (self._alpha_new - 1) / self._alpha_old + self._v_hat = self._v_new + (self._v_new - self._v_old) * beta + self._mu_hat = self._mu_new + (self._mu_new - self._mu_old) * beta else: # reboot to old iteration self._alpha_new = 1 @@ -286,6 +287,5 @@ def _update(self): self.xp.copyto(self._mu_old, self._mu_new) # Test cost function for convergence. if self._cost_func: - self.converge = self.any_convergence_flag() or self._cost_func.get_cost( - self._u_new - ) + self.converge = self.any_convergence_flag() + self.convergd |= self._cost_func.get_cost(self._u_new, self._v_new) From 23776d921027a4dd397197a8a98e163719529189 Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Tue, 15 Nov 2022 17:12:10 +0100 Subject: [PATCH 04/13] Abstract class for cost function. --- modopt/opt/cost.py | 152 +++++++++++++++++++++++++++++++++------------ 1 file changed, 114 insertions(+), 38 deletions(-) diff --git a/modopt/opt/cost.py b/modopt/opt/cost.py index 3cdfcc50..a90c87ae 100644 --- a/modopt/opt/cost.py +++ b/modopt/opt/cost.py @@ -6,6 +6,8 @@ """ +import abc + import numpy as np from modopt.base.backend import get_array_module @@ -13,8 +15,8 @@ from modopt.plot.cost_plot import plotCost -class costObj(object): - """Generic cost function object. +class AbstractcostObj(abc.ABC): + """Abstract cost function object. This class updates the cost according to the input operator classes and tests for convergence. @@ -40,7 +42,8 @@ class costObj(object): Notes ----- - The costFunc class must contain a method called ``cost``. + All child classes should implement a ``_calc_cost`` method (returning + a float) or a ``get_cost`` for more complex behavior on convergence test. Examples -------- @@ -71,7 +74,6 @@ class costObj(object): def __init__( self, - operators, initial_cost=1e6, tolerance=1e-4, cost_interval=1, @@ -80,9 +82,6 @@ def __init__( plot_output=None, ): - self._operators = operators - if not isinstance(operators, type(None)): - self._check_operators() self.cost = initial_cost self._cost_list = [] self._cost_interval = cost_interval @@ -93,30 +92,6 @@ def __init__( self._plot_output = plot_output self._verbose = verbose - def _check_operators(self): - """Check operators. - - This method checks if the input operators have a ``cost`` method. - - Raises - ------ - TypeError - For invalid operators type - ValueError - For operators without ``cost`` method - - """ - if not isinstance(self._operators, (list, tuple, np.ndarray)): - message = ( - 'Input operators must be provided as a list, not {0}' - ) - raise TypeError(message.format(type(self._operators))) - - for op in self._operators: - if not hasattr(op, 'cost'): - raise ValueError('Operators must contain "cost" method.') - op.cost = check_callable(op.cost) - def _check_cost(self): """Check cost function. @@ -167,6 +142,7 @@ def _check_cost(self): return False + @abc.abstractmethod def _calc_cost(self, *args, **kwargs): """Calculate the cost. @@ -178,14 +154,7 @@ def _calc_cost(self, *args, **kwargs): Positional arguments **kwargs : dict Keyword arguments - - Returns - ------- - float - Cost value - """ - return np.sum([op.cost(*args, **kwargs) for op in self._operators]) def get_cost(self, *args, **kwargs): """Get cost function. @@ -241,3 +210,110 @@ def plot_cost(self): # pragma: no cover """ plotCost(self._cost_list, self._plot_output) + + +class costObj(AbstractcostObj): + """Abstract cost function object. + + This class updates the cost according to the input operator classes and + tests for convergence. + + Parameters + ---------- + opertors : list, tuple or numpy.ndarray + List of operators classes containing ``cost`` method + initial_cost : float, optional + Initial value of the cost (default is ``1e6``) + tolerance : float, optional + Tolerance threshold for convergence (default is ``1e-4``) + cost_interval : int, optional + Iteration interval to calculate cost (default is ``1``). + If ``cost_interval`` is ``None`` the cost is never calculated, + thereby saving on computation time. + test_range : int, optional + Number of cost values to be used in test (default is ``4``) + verbose : bool, optional + Option for verbose output (default is ``True``) + plot_output : str, optional + Output file name for cost function plot + + Examples + -------- + >>> from modopt.opt.cost import * + >>> class dummy(object): + ... def cost(self, x): + ... return x ** 2 + ... + ... + >>> inst = costObj([dummy(), dummy()]) + >>> inst.get_cost(2) + - ITERATION: 1 + - COST: 8 + + False + >>> inst.get_cost(2) + - ITERATION: 2 + - COST: 8 + + False + >>> inst.get_cost(2) + - ITERATION: 3 + - COST: 8 + + False + """ + + def __init__( + self, + operators, + **kwargs, + ): + super().__init__(**kwargs) + + self._operators = operators + if not isinstance(operators, type(None)): + self._check_operators() + + def _check_operators(self): + """Check operators. + + This method checks if the input operators have a ``cost`` method. + + Raises + ------ + TypeError + For invalid operators type + ValueError + For operators without ``cost`` method + + """ + if not isinstance(self._operators, (list, tuple, np.ndarray)): + message = ( + 'Input operators must be provided as a list, not {0}' + ) + raise TypeError(message.format(type(self._operators))) + + for op in self._operators: + if not hasattr(op, 'cost'): + raise ValueError('Operators must contain "cost" method.') + op.cost = check_callable(op.cost) + + def _calc_cost(self, *args, **kwargs): + """Calculate the cost. + + This method calculates the cost from each of the input operators. + + Parameters + ---------- + *args : tuple + Positional arguments + **kwargs : dict + Keyword arguments + + Returns + ------- + float + Cost value + + """ + return np.sum([op.cost(*args, **kwargs) for op in self._operators]) From 38fefba4c2aed7f7dc4ba2d1ada2eb2f8b51f4f9 Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Tue, 15 Nov 2022 17:25:26 +0100 Subject: [PATCH 05/13] add custom cost operator for admm. --- modopt/opt/algorithms/admm.py | 72 ++++++++++++++++++++++++++++++----- 1 file changed, 62 insertions(+), 10 deletions(-) diff --git a/modopt/opt/algorithms/admm.py b/modopt/opt/algorithms/admm.py index 1830a846..7fef879b 100644 --- a/modopt/opt/algorithms/admm.py +++ b/modopt/opt/algorithms/admm.py @@ -1,8 +1,66 @@ """ADMM Algorithms.""" import numpy as np +from modopt.base.backend import get_array_module from modopt.opt.algorithms.base import SetUp -from modopt.opt.cost import costObj +from modopt.opt.cost import AbstractcostObj + + +class ADMMcostObj(AbstractcostObj): + r"""Cost Object for the ADMM problem class. + + Compute :math:`f(u)+g(v) + \tau \| Au +Bv - b\|^2` + + Parameters + ---------- + cost_funcs: 2-tuples of callable + f and g function. + A : OperatorBase + First Operator + B : OperatorBase + Second Operator + b : array_like + Observed data + **kwargs : dict + Extra parameters for cost operator configuration + + + See Also + -------- + AbstractcostObj: parent class + """ + + def __init__(self, cost_funcs, A, B, b, tau, **kwargs): + super().__init__(*kwargs) + self.cost_funcs = cost_funcs + self.A = A + self.B = B + self.b = b + self.tau = tau + + def _calc_cost(self, u, v): + """Calculate the cost. + + This method calculates the cost from each of the input operators. + + Parameters + ---------- + u: array_like + First primal variable of ADMM + v: array_like + Second primal variable of ADMM + + Returns + ------- + float + Cost value + + """ + xp = get_array_module(u) + cost = self.cost_func[0](u) + cost += self.cost_func[1](v) + cost += self.tau * xp.linalg.norm(self.A.op(u) + self.B.op(v) - self.b) + return cost class ADMM(SetUp): @@ -64,8 +122,8 @@ def __init__( optimizers, tau=1, max_iter2=5, - cost_func=None, - **kwargs + cost_funcs=None, + **kwargs, ): super().__init__(**kwargs) self.A = A @@ -75,13 +133,7 @@ def __init__( self._opti_G = optimizers[1] self._tau = tau - self._cost_func = costObj() - # patching to get the full cost - self._cost_func._calc_cost = lambda u, v: ( - cost_func[0](u) - + cost_func[1](v) - + self.xp.linalg.norm(A.op(u) + B.op(v) - b) - ) + self._cost_func = ADMMcostObj(cost_funcs, A, B, b, tau) # init iteration variables. self._u_old = self.xp.copy(u) From 3d710444e189a703dcfebf4bfc219bf08472f865 Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Tue, 15 Nov 2022 17:25:46 +0100 Subject: [PATCH 06/13] fix WPS compliance. --- modopt/opt/algorithms/admm.py | 16 ++++++++-------- setup.cfg | 2 ++ 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/modopt/opt/algorithms/admm.py b/modopt/opt/algorithms/admm.py index 7fef879b..5925c9ae 100644 --- a/modopt/opt/algorithms/admm.py +++ b/modopt/opt/algorithms/admm.py @@ -197,9 +197,9 @@ def get_notify_observers_kwargs(self): The mapping between the iterated variables """ return { - "u_new": self._u_new, - "v_new": self._v_new, - "idx": self.idx, + 'u_new': self._u_new, + 'v_new': self._v_new, + 'idx': self.idx, } def retrieve_outputs(self): @@ -209,7 +209,7 @@ def retrieve_outputs(self): y_final, metrics. """ metrics = {} - for obs in self._observers["cv_metrics"]: + for obs in self._observers['cv_metrics']: metrics[obs.name] = obs.retrieve_metrics() self.metrics = metrics @@ -279,7 +279,7 @@ def __init__( opti_H_kwargs=None, opti_G_kwargs=None, cost=None, - **kwargs + **kwargs, ): super().__init__( u=u, @@ -293,10 +293,10 @@ def __init__( opti_H_kwargs=opti_H_kwargs, opti_G_kwargs=opti_G, cost=None, - **kwargs + **kwargs, ) self._c_old = np.inf - self._c_new = 0.0 + self._c_new = 0 self._eta = eta self._alpha_old = alpha self._alpha_new = alpha @@ -320,7 +320,7 @@ def _update(self): # restarting condition self._c_new = self.xp.linalg.norm(self._mu_new - self._mu_hat) self._c_new += self._tau * self.xp.linalg.norm( - self.B.op(self._v_new - self._v_hat) + self.B.op(self._v_new - self._v_hat), ) if self._c_new < self._eta * self._c_old: self._alpha_new = 1 + np.sqrt(1 + 4 * self._alpha_old**2) diff --git a/setup.cfg b/setup.cfg index 87496ced..e5065be8 100644 --- a/setup.cfg +++ b/setup.cfg @@ -42,6 +42,8 @@ per-file-ignores = modopt/opt/algorithms/__init__.py: F401,F403,WPS318, WPS319, WPS412, WPS410 #Todo: x is a too short name. modopt/opt/algorithms/forward_backward.py: WPS111 + #Todo: u,v , A is a too short name. + modopt/opt/algorithms/admm.py: WPS111, N803 #Todo: Check need for del statement modopt/opt/algorithms/primal_dual.py: WPS111, WPS420 #multiline parameters bug with tuples From 53a6bdde917deb495bcf9eee61b1fbfb2428657d Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Thu, 11 May 2023 11:15:04 +0200 Subject: [PATCH 07/13] typos. --- modopt/opt/algorithms/__init__.py | 2 +- modopt/opt/algorithms/admm.py | 5 ++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/modopt/opt/algorithms/__init__.py b/modopt/opt/algorithms/__init__.py index e0ac2572..a4dc4146 100644 --- a/modopt/opt/algorithms/__init__.py +++ b/modopt/opt/algorithms/__init__.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -r"""OPTIMISATION ALGOTITHMS. +r"""OPTIMISATION ALGORITHMS. This module contains class implementations of various optimisation algoritms. diff --git a/modopt/opt/algorithms/admm.py b/modopt/opt/algorithms/admm.py index 5925c9ae..b1163e1b 100644 --- a/modopt/opt/algorithms/admm.py +++ b/modopt/opt/algorithms/admm.py @@ -66,8 +66,7 @@ def _calc_cost(self, u, v): class ADMM(SetUp): r"""Fast ADMM Optimisation Algorihm. - This class implement the ADMM algorithm - (Algorithm 1 from :cite:`Goldstein2014`) + This class implement the ADMM algorithm described in :cite:`Goldstein2014` (Algorithm 1). Parameters ---------- @@ -149,7 +148,7 @@ def _update(self): obs=self.B.op(self._v_old) + self._u_old - self.b, ) tmp = self.A.op(self._u_new) - self._v_new = self.solver2( + self._v_new = self._opti_G( init=self._v_old, obs=tmp + self._u_old - self.c, ) From 8dfa8a35387bbad209e60defa9e1b4cf4570f7de Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Thu, 11 May 2023 11:51:23 +0200 Subject: [PATCH 08/13] feat(test): add test for admm. --- modopt/opt/algorithms/__init__.py | 1 + modopt/opt/algorithms/admm.py | 37 ++++++++++++++++++------------- modopt/tests/test_algorithms.py | 34 ++++++++++++++++++++++++++-- 3 files changed, 55 insertions(+), 17 deletions(-) diff --git a/modopt/opt/algorithms/__init__.py b/modopt/opt/algorithms/__init__.py index a4dc4146..d4e7082b 100644 --- a/modopt/opt/algorithms/__init__.py +++ b/modopt/opt/algorithms/__init__.py @@ -57,3 +57,4 @@ SAGAOptGradOpt, VanillaGenericGradOpt) from modopt.opt.algorithms.primal_dual import Condat +from modopt.opt.algorithms.admm import ADMM, FastADMM diff --git a/modopt/opt/algorithms/admm.py b/modopt/opt/algorithms/admm.py index b1163e1b..6d7a59aa 100644 --- a/modopt/opt/algorithms/admm.py +++ b/modopt/opt/algorithms/admm.py @@ -38,7 +38,7 @@ def __init__(self, cost_funcs, A, B, b, tau, **kwargs): self.b = b self.tau = tau - def _calc_cost(self, u, v): + def _calc_cost(self, u, v, **kwargs): """Calculate the cost. This method calculates the cost from each of the input operators. @@ -57,8 +57,8 @@ def _calc_cost(self, u, v): """ xp = get_array_module(u) - cost = self.cost_func[0](u) - cost += self.cost_func[1](v) + cost = self.cost_funcs[0](u) + cost += self.cost_funcs[1](v) cost += self.tau * xp.linalg.norm(self.A.op(u) + self.B.op(v) - self.b) return cost @@ -70,6 +70,12 @@ class ADMM(SetUp): Parameters ---------- + u: array_like + First primal variable of ADMM + v: array_like + Second primal variable of ADMM + mu: array_like + Lagrangian multiplier. A : OperatorBase Linear operator for u B : OperatorBase @@ -83,10 +89,9 @@ class ADMM(SetUp): .. math:: v_{k+1} = \argmin G(v) + \frac{\tau}{2}\|Bv - y \|^2 cost_funcs = 2-tuple of function Compute the values of H and G - rho : float , optional - regularisation coupling variable default is ``1.0`` - eta : float, optional - Restart threshold, default is ``0.999`` + tau: float, default=1 + Coupling parameter for ADMM. + max_iter2: int Notes ----- @@ -120,7 +125,6 @@ def __init__( b, optimizers, tau=1, - max_iter2=5, cost_funcs=None, **kwargs, ): @@ -131,8 +135,10 @@ def __init__( self._opti_H = optimizers[0] self._opti_G = optimizers[1] self._tau = tau - - self._cost_func = ADMMcostObj(cost_funcs, A, B, b, tau) + if cost_funcs is not None: + self._cost_func = ADMMcostObj(cost_funcs, A, B, b, tau) + else: + self._cost_func = None # init iteration variables. self._u_old = self.xp.copy(u) @@ -150,7 +156,7 @@ def _update(self): tmp = self.A.op(self._u_new) self._v_new = self._opti_G( init=self._v_old, - obs=tmp + self._u_old - self.c, + obs=tmp + self._u_old - self.b, ) self._mu_new = self._mu_old + (tmp + self.B.op(self._v_new) - self.b) @@ -163,7 +169,7 @@ def _update(self): # Test cost function for convergence. if self._cost_func: self.converge = self.any_convergence_flag() - self.converge |= self._cost_func.get_cost(self._u_new, self.v_new) + self.converge |= self._cost_func.get_cost(self._u_new, self._v_new) def iterate(self, max_iter=150): """Iterate. @@ -182,6 +188,7 @@ def iterate(self, max_iter=150): self.retrieve_outputs() # rename outputs as attributes self.u_final = self._u_new + self.x_final = self.u_final # for backward compatibility self.v_final = self._v_new def get_notify_observers_kwargs(self): @@ -196,7 +203,7 @@ def get_notify_observers_kwargs(self): The mapping between the iterated variables """ return { - 'u_new': self._u_new, + 'x_new': self._u_new, 'v_new': self._v_new, 'idx': self.idx, } @@ -309,9 +316,9 @@ def _update(self): obs=self.B.op(self._v_hat) + self._u_old - self.b, ) tmp = self.A.op(self._u_new) - self._v_new = self.solver2( + self._v_new = self._opti_G( init=self._v_hat, - obs=tmp + self._u_hat - self.c, + obs=tmp + self._u_hat - self.b, ) self._mu_new = self._mu_hat + (tmp + self.B.op(self._v_new) - self.b) diff --git a/modopt/tests/test_algorithms.py b/modopt/tests/test_algorithms.py index 73091acd..0cb1670e 100644 --- a/modopt/tests/test_algorithms.py +++ b/modopt/tests/test_algorithms.py @@ -80,7 +80,15 @@ def build_kwargs(kwargs, use_metrics): @parametrize(use_metrics=[True, False]) class AlgoCases: - """Cases for algorithms.""" + """Cases for algorithms. + + Most of the test solves the trivial problem + + .. math:: + \\min_x \\frac{1}{2} \\| y - x \\|_2^2 \\quad\\text{s.t.} x \\geq 0 + + More complex and concrete usecases are shown in examples. + """ data1 = np.arange(9).reshape(3, 3).astype(float) data2 = data1 + np.random.randn(*data1.shape) * 1e-6 @@ -103,7 +111,8 @@ class AlgoCases: ] ) def case_forward_backward(self, kwargs, idty, use_metrics): - """Forward Backward case.""" + """Forward Backward case. + """ update_kwargs = build_kwargs(kwargs, use_metrics) algo = algorithms.ForwardBackward( self.data1, @@ -234,6 +243,27 @@ def case_grad(self, GradDescent, use_metrics, idty): algo.iterate() return algo, update_kwargs + def case_admm(self, use_metrics, idty): + """ADMM setup.""" + def optim1(init, obs): + return obs + + def optim2(init, obs): + return obs + + update_kwargs = build_kwargs({}, use_metrics) + algo = algorithms.ADMM( + u=self.data1, + v=self.data1, + mu=np.zeros_like(self.data1), + A=linear.Identity(), + B=linear.Identity(), + b=self.data1, + optimizers=(optim1, optim2), + **update_kwargs, + ) + algo.iterate() + return algo, update_kwargs @parametrize_with_cases("algo, kwargs", cases=AlgoCases) def test_algo(algo, kwargs): From 6d4e7804b4ecba7fa74008a7dfa900d866af236b Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Thu, 11 May 2023 12:08:23 +0200 Subject: [PATCH 09/13] feat(admm): improve doc. --- modopt/opt/algorithms/admm.py | 75 +++++++++++++++-------------------- 1 file changed, 31 insertions(+), 44 deletions(-) diff --git a/modopt/opt/algorithms/admm.py b/modopt/opt/algorithms/admm.py index 6d7a59aa..f8281093 100644 --- a/modopt/opt/algorithms/admm.py +++ b/modopt/opt/algorithms/admm.py @@ -3,10 +3,10 @@ from modopt.base.backend import get_array_module from modopt.opt.algorithms.base import SetUp -from modopt.opt.cost import AbstractcostObj +from modopt.opt.cost import CostParent -class ADMMcostObj(AbstractcostObj): +class ADMMcostObj(CostParent): r"""Cost Object for the ADMM problem class. Compute :math:`f(u)+g(v) + \tau \| Au +Bv - b\|^2` @@ -27,7 +27,7 @@ class ADMMcostObj(AbstractcostObj): See Also -------- - AbstractcostObj: parent class + CostParent: parent class """ def __init__(self, cost_funcs, A, B, b, tau, **kwargs): @@ -71,18 +71,18 @@ class ADMM(SetUp): Parameters ---------- u: array_like - First primal variable of ADMM + Initial value for first primal variable of ADMM v: array_like - Second primal variable of ADMM + Initial value for second primal variable of ADMM mu: array_like - Lagrangian multiplier. + Initial value for lagrangian multiplier. A : OperatorBase Linear operator for u B : OperatorBase Linear operator for v b : array_like Constraint vector - optimizers: 2-tuple of functions + optimizers: 2-tuple of callable. Solvers for the u and v update, takes init_value and obs_value as argument. each element returns an estimate for: .. math:: u_{k+1} = \argmin H(u) + \frac{\tau}{2}\|A u - y\|^2 @@ -91,7 +91,6 @@ class ADMM(SetUp): Compute the values of H and G tau: float, default=1 Coupling parameter for ADMM. - max_iter2: int Notes ----- @@ -228,45 +227,39 @@ class FastADMM(ADMM): Parameters ---------- + u: array_like + Initial value for first primal variable of ADMM + v: array_like + Initial value for second primal variable of ADMM + mu: array_like + Initial value for lagrangian multiplier. A : OperatorBase Linear operator for u B : OperatorBase Linear operator for v b : array_like Constraint vector - solver1 : function - Solver for the x update, takes init_value and obs_value as argument. - ie, return an estimate for: + optimizers: 2-tuple of callable. + Solvers for the u and v update, takes init_value and obs_value as + argument. each element returns an estimate for: .. math:: u_{k+1} = \argmin H(u) + \frac{\tau}{2}\|A u - y\|^2 - solver2 : function - Solver for the z update, takes init_value and obs_value as argument. - ie return an estimate for: .. math:: v_{k+1} = \argmin G(v) + \frac{\tau}{2}\|Bv - y \|^2 - rho : float , optional - regularisation coupling variable default is ``1.0`` - eta : float, optional - Restart threshold, default is ``0.999`` + cost_funcs = 2-tuple of function + Compute the values of H and G + tau: float, default=1 + Coupling parameter for ADMM. + eta: float, default=0.999 + Convergence parameter for ADMM. + alpha: float, default=1. + Initial value for the FISTA-like acceleration parameter. Notes ----- - The algorithm solve the problem: - - .. math:: u, v = \arg\min H(u) + G(v) + \frac\tau2 \|Au + Bv - b \|_2^2 - - with the following augmented lagrangian: - - .. math:: \mathcal{L}_{\tau}(u,v, \lambda) = H(u) + G(v) - +\langle\lambda |Au + Bv -b \rangle + \frac\tau2 \| Au + Bv -b \|^2 - - To allow easy iterative solving, the change of variable - :math:`\mu=\lambda/\tau` is used. Hence, the lagrangian of interest is: - - .. math :: \tilde{\mathcal{L}}_{\tau}(u,v, \mu) = H(u) + G(v) - + \frac\tau2 \left(\|\mu + Au +Bv - b\|^2 - \|\mu\|^2\right) + This is an accelerated version of the ADMM algorithm. The convergence hypothesis are stronger than for the ADMM algorithm. See Also -------- - SetUp: parent class + ADMM: parent class """ def __init__( @@ -277,14 +270,11 @@ def __init__( A, B, b, - opti_H, - opti_G, + optimizers, + cost_funcs=None, alpha=1, eta=0.999, tau=1, - opti_H_kwargs=None, - opti_G_kwargs=None, - cost=None, **kwargs, ): super().__init__( @@ -294,11 +284,8 @@ def __init__( A=A, B=B, b=b, - opti_H=opti_H, - opti_G=opti_G, - opti_H_kwargs=opti_H_kwargs, - opti_G_kwargs=opti_G, - cost=None, + optimizers=optimizers, + cost_funcs=cost_funcs, **kwargs, ) self._c_old = np.inf @@ -318,7 +305,7 @@ def _update(self): tmp = self.A.op(self._u_new) self._v_new = self._opti_G( init=self._v_hat, - obs=tmp + self._u_hat - self.b, + obs=tmp + self._u_old - self.b, ) self._mu_new = self._mu_hat + (tmp + self.B.op(self._v_new) - self.b) From 843e65b31abb76e889aefbcd32a563936040faa1 Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Thu, 11 May 2023 12:08:39 +0200 Subject: [PATCH 10/13] refactor: rename abstract cost to CostParent. --- modopt/opt/cost.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modopt/opt/cost.py b/modopt/opt/cost.py index a90c87ae..688a3959 100644 --- a/modopt/opt/cost.py +++ b/modopt/opt/cost.py @@ -15,7 +15,7 @@ from modopt.plot.cost_plot import plotCost -class AbstractcostObj(abc.ABC): +class CostParent(abc.ABC): """Abstract cost function object. This class updates the cost according to the input operator classes and @@ -212,7 +212,7 @@ def plot_cost(self): # pragma: no cover plotCost(self._cost_list, self._plot_output) -class costObj(AbstractcostObj): +class costObj(CostParent): """Abstract cost function object. This class updates the cost according to the input operator classes and From 676ecc0c8d6706c136f7ba92442aa06a252df785 Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Thu, 11 May 2023 12:09:02 +0200 Subject: [PATCH 11/13] feat: add test for fast admm. --- modopt/tests/test_algorithms.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/modopt/tests/test_algorithms.py b/modopt/tests/test_algorithms.py index 0cb1670e..5671b8e3 100644 --- a/modopt/tests/test_algorithms.py +++ b/modopt/tests/test_algorithms.py @@ -242,8 +242,8 @@ def case_grad(self, GradDescent, use_metrics, idty): ) algo.iterate() return algo, update_kwargs - - def case_admm(self, use_metrics, idty): + @parametrize(admm=[algorithms.ADMM,algorithms.FastADMM]) + def case_admm(self, admm, use_metrics, idty): """ADMM setup.""" def optim1(init, obs): return obs @@ -252,7 +252,7 @@ def optim2(init, obs): return obs update_kwargs = build_kwargs({}, use_metrics) - algo = algorithms.ADMM( + algo = admm( u=self.data1, v=self.data1, mu=np.zeros_like(self.data1), From 9af7fc40d6394461f70cff747fbe8555a73d5d6d Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Thu, 11 May 2023 12:20:22 +0200 Subject: [PATCH 12/13] feat(admm): improve docstrings. --- modopt/opt/algorithms/admm.py | 54 ++++++++++++++++++----------------- 1 file changed, 28 insertions(+), 26 deletions(-) diff --git a/modopt/opt/algorithms/admm.py b/modopt/opt/algorithms/admm.py index f8281093..0524ae47 100644 --- a/modopt/opt/algorithms/admm.py +++ b/modopt/opt/algorithms/admm.py @@ -9,7 +9,6 @@ class ADMMcostObj(CostParent): r"""Cost Object for the ADMM problem class. - Compute :math:`f(u)+g(v) + \tau \| Au +Bv - b\|^2` Parameters ---------- @@ -19,11 +18,14 @@ class ADMMcostObj(CostParent): First Operator B : OperatorBase Second Operator - b : array_like + b : numpy.ndarray Observed data **kwargs : dict Extra parameters for cost operator configuration + Notes + ----- + Compute :math:`f(u)+g(v) + \tau \| Au +Bv - b\|^2` See Also -------- @@ -45,9 +47,9 @@ def _calc_cost(self, u, v, **kwargs): Parameters ---------- - u: array_like + u: numpy.ndarray First primal variable of ADMM - v: array_like + v: numpy.ndarray Second primal variable of ADMM Returns @@ -70,25 +72,25 @@ class ADMM(SetUp): Parameters ---------- - u: array_like + u: numpy.ndarray Initial value for first primal variable of ADMM - v: array_like + v: numpy.ndarray Initial value for second primal variable of ADMM - mu: array_like + mu: numpy.ndarray Initial value for lagrangian multiplier. - A : OperatorBase + A : modopt.opt.linear.LinearOperator Linear operator for u - B : OperatorBase + B: modopt.opt.linear.LinearOperator Linear operator for v - b : array_like + b : numpy.ndarray Constraint vector - optimizers: 2-tuple of callable. - Solvers for the u and v update, takes init_value and obs_value as - argument. each element returns an estimate for: + optimizers: tuple + 2-tuple of callable, that are the optimizers for the u and v. + Each callable should access an init and obs argument and returns an estimate for: .. math:: u_{k+1} = \argmin H(u) + \frac{\tau}{2}\|A u - y\|^2 .. math:: v_{k+1} = \argmin G(v) + \frac{\tau}{2}\|Bv - y \|^2 - cost_funcs = 2-tuple of function - Compute the values of H and G + cost_funcs: tuple + 2-tuple of callable, that compute values of H and G. tau: float, default=1 Coupling parameter for ADMM. @@ -227,25 +229,25 @@ class FastADMM(ADMM): Parameters ---------- - u: array_like + u: numpy.ndarray Initial value for first primal variable of ADMM - v: array_like + v: numpy.ndarray Initial value for second primal variable of ADMM - mu: array_like + mu: numpy.ndarray Initial value for lagrangian multiplier. - A : OperatorBase + A : modopt.opt.linear.LinearOperator Linear operator for u - B : OperatorBase + B: modopt.opt.linear.LinearOperator Linear operator for v - b : array_like + b : numpy.ndarray Constraint vector - optimizers: 2-tuple of callable. - Solvers for the u and v update, takes init_value and obs_value as - argument. each element returns an estimate for: + optimizers: tuple + 2-tuple of callable, that are the optimizers for the u and v. + Each callable should access an init and obs argument and returns an estimate for: .. math:: u_{k+1} = \argmin H(u) + \frac{\tau}{2}\|A u - y\|^2 .. math:: v_{k+1} = \argmin G(v) + \frac{\tau}{2}\|Bv - y \|^2 - cost_funcs = 2-tuple of function - Compute the values of H and G + cost_funcs: tuple + 2-tuple of callable, that compute values of H and G. tau: float, default=1 Coupling parameter for ADMM. eta: float, default=0.999 From fd0a4f9a2a1e4c851d3bcdebda622c4d486bdfa7 Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Thu, 11 May 2023 13:32:47 +0200 Subject: [PATCH 13/13] style: remove extra line.c --- modopt/opt/algorithms/admm.py | 1 - 1 file changed, 1 deletion(-) diff --git a/modopt/opt/algorithms/admm.py b/modopt/opt/algorithms/admm.py index 0524ae47..b881b770 100644 --- a/modopt/opt/algorithms/admm.py +++ b/modopt/opt/algorithms/admm.py @@ -9,7 +9,6 @@ class ADMMcostObj(CostParent): r"""Cost Object for the ADMM problem class. - Parameters ---------- cost_funcs: 2-tuples of callable