Module regpy.solvers.linear.primal_dual
Classes
class PDHG (setting,
init_domain=None,
init_codomain_star=None,
tau=0,
sigma=0,
theta=1,
proximal_pars_data_fidelity_conjugate=None,
proximal_pars_penalty=None,
compute_y=True,
compute_gap=True,
logging_level=20)-
Expand source code
class PDHG(RegSolver): r"""The Primal-Dual Hybrid Gradient (PDHG) or Chambolle-Pock Algorithm For \(\theta=0)\ this is the Arrow-Hurwicz-Uzawa algorithm. Solves the minimization problem: \(\frac{1}{\alpha}\mathcal{S}_{g^{\delta}}(Tf)+\mathcal{R}(f))\ by solving the saddle-point problem: \[ \inf_f \sup_p [ - \langle Tf,p\rangle + \mathcal{R}(f)- \frac{1}{\alpha}\mathcal{S}_{g^{\delta}}^\ast(-\alpha p) ]. \] Here \(\mathcal{S}_{g^{\delta}}^\ast)\ denotes the Fenchel conjugate functional. Note: Due to a different sign convention for the dual variables, some signs in the iteration formula differ from the originial paper and most of the literature. Parameters ---------- setting : regpy.solvers.TikhonovRegularizationSetting The setting of the forward problem. The operator needs to be linear and "penalty.proximal" and "data_fid.conj.proximal" need to be implemented. init_domain : setting.op.domain [default: None] The initial guess "f". If None, it will be either initialized by 0 or using the optimality conditions in case init_codomain_star is given. init_codomain_star : setting.op.codomain [default: None] The initial guess "p". Initialization is analogous to init_domain. tau : float [default: 0] The parameter to compute the proximal operator of the penalty term. Stepsize of the primal step. Must be non-negative. If 0, a positive value is selected automatically based on the operator norm and the value of sigma. sigma : float [default: 0] The parameter to compute the proximal operator of the data-fidelity term. Stepsize of the dual step. Must be non-negative. If 0, a positive value is selected automatically based on the operator norm and the value of tau theta : float [default: 1] Relaxation parameter. For theta==0 PDHG is the Arrow-Hurwicz-Uzawa algorithm. proximal_pars_data_fidelity_conjugate : dict, optional Parameter dictionary passed to the computation of the prox-operator of the data fidelity functional. proximal_pars_penalty : dict, optional Parameter dictionary passed to the computation of the prox-operator of the penalty functional. compute_y : boolean [True] If True, the images y_k=T(x_k) are computed in each iteration. As they are not needed in the algorithm, so this may considerably increase computational costs. If False, None is returned for y_k. compute_gap : boolean [default: True] If True the duality gap is computed each Term. For this 'penalty.conj' and 'data_fid.conj' need to be implemented, if not a warning is frown and this is set to False. it is not nessesary to compute the gap for the algorithm, but it might be used in stopping criteria. """ def __init__(self, setting, init_domain=None, init_codomain_star=None, tau = 0, sigma = 0, theta= 1, proximal_pars_data_fidelity_conjugate = None, proximal_pars_penalty = None, compute_y = True,compute_gap =True, logging_level = logging.INFO ): assert isinstance(setting, TikhonovRegularizationSetting) super().__init__(setting) assert self.op.linear assert init_domain is None or init_domain in self.op.domain assert init_codomain_star is None or init_codomain_star in self.op.codomain self.log.setLevel(logging_level) if init_domain is None: if init_codomain_star is None: self.x = setting.op.domain.zeros() self.pstar = setting.op.codomain.zeros() else: self.pstar = init_codomain_star self.x = setting.dualToPrimal(self.pstar) else: self.x = init_domain if init_codomain_star is None: self.pstar = setting.primalToDual(self.x) else: self.pstar = init_codomain_star self.x_old = self.x self.compute_y = compute_y self.compute_gap = compute_gap self.y = self.op(self.x) if self.compute_y else None assert tau>=0 and sigma>=0 L = self.setting.op.norm(self.setting.h_domain,self.setting.h_codomain) if tau==0 and sigma==0: self.tau = 1/L self.sigma = 1/L elif tau==0 and sigma>0: self.tau = 1./(L**2*self.sigma) self.sigma = sigma elif sigma==0 and tau>0: self.sigma = 1./(L**2*self.tau) self.tau = tau else: self.sigma = sigma self.tau = tau self.muR = setting.penalty.convexity_param self.muSstar = self.regpar/setting.data_fid.Lipschitz if self.muR>0: if self.muSstar>0: self.mu = 2*np.sqrt(self.muR * self.muSstar)/L self.tau = self.mu/(2.*self.muR) self.sigma = self.mu/(2.*self.muSstar) self.theta = 1./(1.+self.mu) self.log.info('Using accelerated version 2 with convexity parameters mu_R={:.3e}, mu_S*={:.3e} and ||T||={:.3e}.\n Expected linear convergence rate: {:.3e}'.format(self.muR,self.muSstar,L,(1.+self.theta)/(2.+self.mu))) else: self.theta = 0 self.log.info('Using accelerated version 1 with convexity parameter mu_R={:.3e} and ||T|={:.3e}. Expected convergence rate O(1/n^2).'.format(self.muR,L)) else: self.theta = theta self.log.info('Using unaccelerated version') self.proximal_pars_data_fidelity_conjugate = proximal_pars_data_fidelity_conjugate self.proximal_pars_penalty = proximal_pars_penalty if self.compute_gap: try: self.gap = self.setting.dualityGap(primal=self.x) except: self.log.warning('missing implemtation in fuctionals to compute duality gap. The gap is not computed') self.compute_gap = False self.gap = None def _next(self): primal_step = self.x + self.tau * self.h_domain.gram_inv(self.op.adjoint(self.pstar)) self.x = self.penalty.proximal(primal_step, self.tau, self.proximal_pars_penalty) self.y = self.op(self.x) if self.compute_y else None dual_step = -self.pstar + self.sigma * self.h_codomain.gram(self.op( self.x+self.theta*(self.x-self.x_old) )) self.pstar = (-1./self.regpar)*self.data_fid.conj.proximal(self.regpar*dual_step, self.regpar*self.sigma, self.proximal_pars_data_fidelity_conjugate) self.x_old = self.x if self.muR>0 and self.muSstar==0: self.theta = 1./np.sqrt(1+self.muR*self.tau) self.tau *= self.theta self.sigma /= self.theta if self.compute_gap: self.gap = self.setting.dualityGap(primal=self.x,dual= self.pstar)
The Primal-Dual Hybrid Gradient (PDHG) or Chambolle-Pock Algorithm For (\theta=0)\ this is the Arrow-Hurwicz-Uzawa algorithm.
Solves the minimization problem: (\frac{1}{\alpha}\mathcal{S}{g^{\delta}}(Tf)+\mathcal{R}(f))\ by solving the saddle-point problem: \inf_f \sup_p [ - \langle Tf,p\rangle + \mathcal{R}(f)- \frac{1}{\alpha}\mathcal{S}_{g^{\delta}}^\ast(-\alpha p) ]. Here (\mathcal{S}^\ast)\ denotes the Fenchel conjugate functional.}
Note: Due to a different sign convention for the dual variables, some signs in the iteration formula differ from the originial paper and most of the literature.
Parameters
setting
:TikhonovRegularizationSetting
- The setting of the forward problem. The operator needs to be linear and "penalty.proximal" and "data_fid.conj.proximal" need to be implemented.
init_domain
:setting.op.domain [default: None]
- The initial guess "f". If None, it will be either initialized by 0 or using the optimality conditions in case init_codomain_star is given.
init_codomain_star
:setting.op.codomain [default: None]
- The initial guess "p". Initialization is analogous to init_domain.
tau
:float [default: 0]
- The parameter to compute the proximal operator of the penalty term. Stepsize of the primal step. Must be non-negative. If 0, a positive value is selected automatically based on the operator norm and the value of sigma.
sigma
:float [default: 0]
- The parameter to compute the proximal operator of the data-fidelity term. Stepsize of the dual step. Must be non-negative. If 0, a positive value is selected automatically based on the operator norm and the value of tau
theta
:float [default: 1]
- Relaxation parameter. For theta==0 PDHG is the Arrow-Hurwicz-Uzawa algorithm.
proximal_pars_data_fidelity_conjugate
:dict
, optional- Parameter dictionary passed to the computation of the prox-operator of the data fidelity functional.
proximal_pars_penalty
:dict
, optional- Parameter dictionary passed to the computation of the prox-operator of the penalty functional.
compute_y
:boolean [True]
- If True, the images y_k=T(x_k) are computed in each iteration. As they are not needed in the algorithm, so this may considerably increase computational costs. If False, None is returned for y_k.
compute_gap
:boolean [default: True]
- If True the duality gap is computed each Term. For this 'penalty.conj' and 'data_fid.conj' need to be implemented, if not a warning is frown and this is set to False. it is not nessesary to compute the gap for the algorithm, but it might be used in stopping criteria.
Ancestors
Inherited members
class DouglasRachford (setting,
init_h,
tau=1,
regpar=1,
proximal_pars_data_fidelity=None,
proximal_pars_penalty=None)-
Expand source code
class DouglasRachford(RegSolver): r"""The Douglas-Rashford Splitting Algorithm Minimizes \(\mathcal{S}(Tf)+\alpha*\mathcal{R}(f)\) Parameters ---------- setting : regpy.solvers.RegularizationSetting The setting of the forward problem, both penalty and data fidelity need prox-operators. The operator needs to be linear. And the data_fid term contains the the operator for example `data_fid = HilbertNorm(h_space=L2) * (op - data)`, i.e. it is mapping from the domain of the operator. init_h : array_like The initial guess "f". Must be in setting.op.domain. tau : float , optional The parameter to compute the proximal operator of the penalty term. Must be positive. (Default: 1) regpar : float, optional The regularization parameter. Must be positive. (Default: 1) proximal_pars_data_fidelity : dict, optional Parameter dictionary passed to the computation of the prox-operator of the data fidelity functional. (Default: None) proximal_pars_penalty : dict, optional Parameter dictionary passed to the computation of the prox-operator of the penalty functional. (Default: None)) """ def __init__(self, setting, init_h, tau = 1, regpar = 1, proximal_pars_data_fidelity = None, proximal_pars_penalty = None): super().__init__(setting) assert init_h in self.op.domain self.h = init_h self.tau = tau self.regpar = regpar self.proximal_pars_data_fidelity = proximal_pars_data_fidelity self.proximal_pars_penalty = proximal_pars_penalty self.x = self.penalty.proximal(self.h, self.tau*self.regpar, self.proximal_pars_penalty) self.y = self.op(self.x) def _next(self): self.h += self.data_fid.proximal(2*self.x-self.h, self.tau, self.proximal_pars_data_fidelity) - self.x self.x = self.penalty.proximal(self.h, self.tau*self.regpar, self.proximal_pars_penalty) self.y = self.op(self.x)
The Douglas-Rashford Splitting Algorithm
Minimizes \mathcal{S}(Tf)+\alpha*\mathcal{R}(f)
Parameters
setting
:RegularizationSetting
- The setting of the forward problem, both penalty and data fidelity need prox-operators. The operator needs to be linear.
And the data_fid term contains the the operator for example
data_fid = HilbertNorm(h_space=L2) * (op - data)
, i.e. it is mapping from the domain of the operator. init_h
:array_like
- The initial guess "f". Must be in setting.op.domain.
tau
:float
, optional- The parameter to compute the proximal operator of the penalty term. Must be positive. (Default: 1)
regpar
:float
, optional- The regularization parameter. Must be positive. (Default: 1)
proximal_pars_data_fidelity
:dict
, optional- Parameter dictionary passed to the computation of the prox-operator of the data fidelity functional. (Default: None)
proximal_pars_penalty
:dict
, optional- Parameter dictionary passed to the computation of the prox-operator of the penalty functional. (Default: None))
Ancestors
Inherited members