From acd1d77d380aec1fde34ada822bdd6bbd9be3b4f Mon Sep 17 00:00:00 2001 From: Thomas Moreau Date: Tue, 19 Sep 2023 07:31:58 +0900 Subject: [PATCH] MTN use benchopt 1.5 API (#56) --- objective.py | 17 +++++++++-------- solvers/ADMM.py | 8 ++++---- solvers/Celer.py | 2 +- solvers/ChambollePockPDsplit.py | 15 +++++++-------- solvers/CondatVu.py | 13 ++++++------- solvers/DPGD.py | 13 ++++++------- solvers/FP.py | 23 +++++++++++------------ solvers/ISTA.py | 25 ++++++++++++------------- solvers/PGD.py | 27 ++++++++++++++------------- solvers/skglm.py | 2 +- 10 files changed, 71 insertions(+), 74 deletions(-) diff --git a/objective.py b/objective.py index ea6939c..03a6cdc 100644 --- a/objective.py +++ b/objective.py @@ -9,13 +9,14 @@ class Objective(BaseObjective): - min_benchopt_version = "1.3" name = "TV1D" + min_benchopt_version = "1.5" - parameters = {'reg': [0.5], - 'delta': [0.9], - 'data_fit': ['quad', 'huber'] - } + parameters = { + 'reg': [0.5], + 'delta': [0.9], + 'data_fit': ['quad', 'huber'] + } def set_data(self, A, y, x): self.A, self.y, self.x = A, y, x @@ -23,7 +24,7 @@ def set_data(self, A, y, x): self.c = self.get_c(S, self.delta) self.reg_scaled = self.reg*self.get_reg_max(self.c) - def compute(self, u): + def evaluate_result(self, u): R = self.y - self.A @ u reg_TV = abs(np.diff(u)).sum() if self.data_fit == 'quad': @@ -34,8 +35,8 @@ def compute(self, u): return dict(value=loss + self.reg_scaled * reg_TV, norm_x=norm_x) - def get_one_solution(self): - return np.zeros(self.A.shape[1]) + def get_one_result(self): + return dict(u=np.zeros(self.A.shape[1])) def get_objective(self): return dict(A=self.A, reg=self.reg_scaled, y=self.y, c=self.c, diff --git a/solvers/ADMM.py b/solvers/ADMM.py index ab9b420..d161f54 100644 --- a/solvers/ADMM.py +++ b/solvers/ADMM.py @@ -44,7 +44,7 @@ def set_objective(self, A, reg, y, c, delta, data_fit): def run(self, callback): p = self.A.shape[1] - u = self.c * np.ones(p) + self.u = u = self.c * np.ones(p) z = np.zeros(p - 1) mu = np.zeros(p - 1) gamma = self.gamma @@ -63,7 +63,7 @@ def run(self, callback): self.A @ x) - gamma * np.diff( np.diff(x), append=0, prepend=0)) - while callback(u): + while callback(): z_old = z if self.data_fit == 'quad': u_tmp = (Aty + np.diff(mu, append=0, prepend=0) @@ -95,7 +95,7 @@ def jac(u): gamma *= 2 if s > 10 * r: gamma /= 2 - self.u = u + self.u = u def get_result(self): - return self.u + return dict(u=self.u) diff --git a/solvers/Celer.py b/solvers/Celer.py index 9dc052f..5e95111 100644 --- a/solvers/Celer.py +++ b/solvers/Celer.py @@ -71,4 +71,4 @@ def get_next(previous): return previous + 1 def get_result(self): - return self.u + return dict(u=self.u) diff --git a/solvers/ChambollePockPDsplit.py b/solvers/ChambollePockPDsplit.py index f24a34f..b911aa1 100644 --- a/solvers/ChambollePockPDsplit.py +++ b/solvers/ChambollePockPDsplit.py @@ -36,13 +36,13 @@ def run(self, callback): sigma_v = 1.0 / (self.ratio * LD) sigma_w = 1.0 / (self.ratio * LA) # Init variables - u = np.zeros(p) v = np.zeros(p - 1) # we consider non-cyclic finite difference w = np.zeros(n) - u_bar = u.copy() + self.u = np.zeros(p) + u_bar = self.u.copy() - while callback(u): - u_old = u + while callback(): + u_old = self.u v = np.clip(v + sigma_v * np.diff(u_bar), -self.reg, self.reg) w_tmp = w + sigma_w * self.A @ u_bar @@ -55,12 +55,11 @@ def run(self, callback): else: w = (w_tmp - sigma_w * self.y) / (1.0 + sigma_w) # grad.T = -div, hence + sign - u = u + tau * np.diff(v, prepend=0, append=0) - tau * self.A.T @ w - u_bar = u + self.theta * (u - u_old) - self.u = u + self.u += tau * (np.diff(v, prepend=0, append=0) - self.A.T @ w) + u_bar = self.u + self.theta * (self.u - u_old) def get_result(self): - return self.u + return dict(u=self.u) def _prox_huber(self, u, mu): return np.where( diff --git a/solvers/CondatVu.py b/solvers/CondatVu.py index a93c64d..51077f6 100644 --- a/solvers/CondatVu.py +++ b/solvers/CondatVu.py @@ -36,20 +36,19 @@ def run(self, callback): tau = 1 / (LA ** 2 / 2 + sigma * LD ** 2) eta = self.eta # initialisation - u = self.c * np.ones(p) + self.u = self.c * np.ones(p) v = np.zeros(p - 1) - while callback(u): - u_tmp = (u - tau * self.grad(self.A, u) + while callback(): + u_tmp = (self.u - tau * self.grad(self.A, self.u) - tau * (-np.diff(v, append=0, prepend=0))) - v_tmp = np.clip(v + sigma * np.diff(2 * u_tmp - u), + v_tmp = np.clip(v + sigma * np.diff(2 * u_tmp - self.u), -self.reg, self.reg) - u = eta * u_tmp + (1 - eta) * u + self.u = eta * u_tmp + (1 - eta) * self.u v = eta * v_tmp + (1 - eta) * v - self.u = u def get_result(self): - return self.u + return dict(u=self.u) def grad(self, A, u): R = self.y - A @ u diff --git a/solvers/DPGD.py b/solvers/DPGD.py index 137e896..c183225 100644 --- a/solvers/DPGD.py +++ b/solvers/DPGD.py @@ -57,14 +57,14 @@ def run(self, callback): # alpha / rho stepsize = self.alpha / (np.linalg.norm(DA_inv, ord=2)**2) # initialisation - u = self.c * np.ones(p) + self.u = self.c * np.ones(p) v = np.zeros(p - 1) v_tmp = np.zeros(p) v_old = v.copy() v_acc = v.copy() t_new = 1 - while callback(u): + while callback(): if self.use_acceleration: t_old = t_new t_new = (1 + np.sqrt(1 + 4 * t_old ** 2)) / 2 @@ -83,11 +83,10 @@ def run(self, callback): v_acc[:] = v + (t_old - 1.) / t_new * (v - v_old) if isinstance(self.A, np.ndarray): - u = AtA_inv @ (Aty + np.diff(v, append=0, prepend=0)) + self.u = AtA_inv @ (Aty + np.diff(v, append=0, prepend=0)) else: - u, _ = cg(AtA, Aty + np.diff(v, append=0, prepend=0), - x0=u, tol=tol_cg) - self.u = u + self.u, _ = cg(AtA, Aty + np.diff(v, append=0, prepend=0), + x0=self.u, tol=tol_cg) def get_result(self): - return self.u + return dict(u=self.u) diff --git a/solvers/FP.py b/solvers/FP.py index cb6f2cb..64687bf 100644 --- a/solvers/FP.py +++ b/solvers/FP.py @@ -34,30 +34,29 @@ def run(self, callback): # alpha / rho stepsize = self.alpha / (n * np.max((AL**2).sum(axis=1))) # initialisation - z = np.zeros(p) - z[0] = self.c + self.z = np.zeros(p) + self.z[0] = self.c mu = np.zeros((p, p)) nu = np.zeros(p) - z_old = z.copy() - z_acc = z.copy() + z_old = self.z.copy() + z_acc = self.z.copy() t_new = 1 - while callback(np.cumsum(z)): + while callback(): if self.use_acceleration: t_old = t_new t_new = (1 + np.sqrt(1 + 4 * t_old ** 2)) / 2 - z_old[:] = z - z[:] = z_acc - mu = z - stepsize * (n * self.grad(AL, z) * AL.T).T + z_old[:] = self.z + self.z[:] = z_acc + mu = self.z - stepsize * (n * self.grad(AL, self.z) * AL.T).T nu = np.mean(mu, axis=0) - z = prox_z(nu, stepsize * self.reg) + self.z = prox_z(nu, stepsize * self.reg) if self.use_acceleration: - z_acc[:] = z + (t_old - 1.) / t_new * (z - z_old) - self.u = np.cumsum(z) + z_acc[:] = self.z + (t_old - 1.) / t_new * (self.z - z_old) def get_result(self): - return self.u + return dict(u=np.cumsum(self.z)) def grad(self, A, u): R = self.y - A @ u diff --git a/solvers/ISTA.py b/solvers/ISTA.py index 15de3e0..3f01a36 100644 --- a/solvers/ISTA.py +++ b/solvers/ISTA.py @@ -11,7 +11,7 @@ class Solver(BaseSolver): """Proximal gradient descent for synthesis formulation.""" name = 'Primal PGD synthesis (ISTA)' - stopping_strategy = 'callback' + sampling_strategy = 'callback' # any parameter defined here is accessible as a class attribute parameters = {'alpha': [1., 1.9], @@ -36,26 +36,25 @@ def run(self, callback): # alpha / rho stepsize = self.alpha / (np.linalg.norm(AL, ord=2)**2) # initialisation - z = np.zeros(p) - z[0] = self.c - z_old = z.copy() - z_acc = z.copy() + self.z = np.zeros(p) + self.z[0] = self.c + z_old = self.z.copy() + z_acc = self.z.copy() t_new = 1 - while callback(np.cumsum(z)): + while callback(): if self.use_acceleration: t_old = t_new t_new = (1 + np.sqrt(1 + 4 * t_old ** 2)) / 2 - z_old[:] = z - z[:] = z_acc - z = prox_z(z - stepsize * self.grad(AL, z), - self.reg * stepsize) + z_old[:] = self.z + self.z[:] = z_acc + self.z = prox_z(self.z - stepsize * self.grad(AL, self.z), + self.reg * stepsize) if self.use_acceleration: - z_acc[:] = z + (t_old - 1.) / t_new * (z - z_old) - self.u = np.cumsum(z) + z_acc[:] = self.z + (t_old - 1.) / t_new * (self.z - z_old) def get_result(self): - return self.u + return dict(u=np.cumsum(self.z)) def grad(self, A, u): R = self.y - A @ u diff --git a/solvers/PGD.py b/solvers/PGD.py index 95ecc59..c05d2a3 100644 --- a/solvers/PGD.py +++ b/solvers/PGD.py @@ -23,8 +23,10 @@ class Solver(BaseSolver): ) # any parameter defined here is accessible as a class attribute - parameters = {'alpha': [1.], - 'use_acceleration': [False, True]} + parameters = { + 'alpha': [1.], + 'use_acceleration': [False, True] + } def set_objective(self, A, reg, y, c, delta, data_fit): self.reg = reg @@ -38,26 +40,25 @@ def run(self, callback): # alpha / rho stepsize = self.alpha / get_l2norm(self.A)**2 # initialisation - u = self.c * np.ones(p) - u_acc = u.copy() - u_old = u.copy() + self.u = self.c * np.ones(p) + u_acc = self.u.copy() + u_old = self.u.copy() t_new = 1 - while callback(u): + while callback(): if self.use_acceleration: t_old = t_new t_new = (1 + np.sqrt(1 + 4 * t_old ** 2)) / 2 - u_old[:] = u - u[:] = u_acc - u = ptv.tv1_1d( - u - stepsize * self.grad(self.A, u), + u_old[:] = self.u + self.u[:] = u_acc + self.u = ptv.tv1_1d( + self.u - stepsize * self.grad(self.A, self.u), self.reg * stepsize, method='condat') if self.use_acceleration: - u_acc[:] = u + (t_old - 1.) / t_new * (u - u_old) - self.u = u + u_acc[:] = self.u + (t_old - 1.) / t_new * (self.u - u_old) def get_result(self): - return self.u + return dict(u=self.u) def grad(self, A, u): R = self.y - A @ u diff --git a/solvers/skglm.py b/solvers/skglm.py index 8abe55c..ccc623f 100644 --- a/solvers/skglm.py +++ b/solvers/skglm.py @@ -77,4 +77,4 @@ def get_next(previous): return previous + 1 def get_result(self): - return self.u + return dict(u=self.u)