Refitting solutions promoted by $\ell_{12}$ sparse analysis regularization with block penalties

TL;DR

This paper introduces a novel refitting framework that removes bias from $\\ell_{12}$ analysis regularized solutions in inverse problems, using block penalties that act on the co-support to improve estimation accuracy.

Contribution

It proposes a new refitting block penalty and an efficient algorithm for bias removal in $\\ell_{12}$ regularized inverse problems, enhancing solution accuracy.

Findings

The new block penalty effectively reduces bias in estimations.

The algorithm efficiently computes both biased and refitted solutions.

Experiments demonstrate improved estimation quality with the proposed method.

Abstract

In inverse problems, the use of an $\ell_{12}$ analysis regularizer induces a bias in the estimated solution. We propose a general refitting framework for removing this artifact while keeping information of interest contained in the biased solution. This is done through the use of refitting block penalties that only act on the co-support of the estimation. Based on an analysis of related works in the literature, we propose a new penalty that is well suited for refitting purposes. We also present an efficient algorithmic method to obtain the refitted solution along with the original (biased) solution for any convex refitting block penalty. Experiments illustrate the good behavior of the proposed block penalty for refitting.

Tables

Table 1. Table 1: Properties satisfied by the considered block penalties ϕ italic-ϕ \phi .

Properties	HO	HD	QO	SD
1	$\surd$	$\surd$	$\surd$	$\surd$
2		$\surd$		$\surd$
3			$\surd$	$\surd$

Table 2. Table 2: Convex conjugates and proximal operators of the studied block penalties ϕ italic-ϕ \phi .

$\varphi$	${\varphi }^{\ast }(z,\hat{z})$	${\mathrm{prox}}_{\sigma {\varphi }^{\ast }}(z_0,\hat{z})$
HO		$z_0-P_{\hat{z}}\left(z_0\right)$
HD
QO		$\frac{\lambda }{\lambda +\sigma \Vert \hat{z}\Vert }\left(z_0-P_{\hat{z}}\left(z_0\right)\right)$
SD		$\lambda \left(\frac{z_0+\lambda \frac{\hat{z}}{\Vert \hat{z}\Vert }}{\max (\lambda ,\Vert z_0+\lambda \frac{\hat{z}}{\Vert \hat{z}\Vert }\Vert )}-\frac{\hat{z}}{\Vert \hat{z}\Vert }\right)$

Equations

$$\hat{x}\in \underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert \Phi x-y{\Vert }^2+\lambda \Vert \Gamma x{\Vert }_{1,2}.$$

$$\mathrm{supp}(\Gamma x)=\left\{1\le i\le m:(\Gamma x{)}_i\ne 0_b\right\}.$$

$${\tilde{x}}^{\mathrm{supp}}\in \underset{x;\mathrm{supp}(\Gamma x)\subseteq \hat{\mathcal{I}}}{\mathrm{argmin}}\frac{1}{2}\Vert \Phi x-y{\Vert }^2,\text{vspace}-0.1cm$$

$$\left\{\begin{array}[]{ll}\hat{z}_{i}=\tfrac{(\Gamma\hat{x})_{i}}{\|(\Gamma\hat{x})_{i}\|}&\textrm{if }i\in\hat{\mathcal{I}}=\operatorname{supp}(\Gamma\hat{x})\enspace,\\ \|\hat{z}_{i}\|\leq 1&\textrm{otherwise}\enspace,\end{array}\right.$$

$$D_i(\Gamma x)=\Vert (\Gamma x{)}_i\Vert -⟨{\hat{z}}_i,(\Gamma x{)}_i⟩=D_{\Vert \cdot \Vert }((\Gamma x{)}_i,(\Gamma \hat{x}{)}_i),$$

$$D_{\Vert \Gamma \cdot {\Vert }_{1,2}}(x,\hat{x})=\sum _{i=1}^m\left(\Vert (\Gamma x{)}_i\Vert -⟨{\hat{z}}_i,(\Gamma x{)}_i⟩\right)=\sum _{i=1}^m{D}_i(\Gamma x).\text{vspace}-0.2cm$$

$$\forall i\in \hat{\mathcal{I}},D_i(\Gamma x)=0\text{ iff }\exists {\alpha }_i\ge 0\text{ such that }(\Gamma x{)}_i={\alpha }_i(\Gamma \hat{x}{)}_i.$$

$$\forall i\in {\hat{\mathcal{I}}}^c,D_i(\Gamma x)=0\text{ with }{\hat{z}}_i\in \partial \Vert \cdot \Vert ((\Gamma \hat{x}{)}_i)\text{ s.t. }\Vert {\hat{z}}_i\Vert <1\text{ iff }(\Gamma x{)}_i=0_b.$$

$${\tilde{x}}_{l+1}\in \underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert \Phi x-y{\Vert }^2+\lambda D_{\Vert \Gamma \cdot {\Vert }_{1,2}}(x,{\tilde{x}}_l).$$

$${\tilde{x}}^{\mathrm{IB}}={\tilde{x}}_2\in \underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert \Phi x-y{\Vert }^2+\lambda D_{\Vert \Gamma \cdot {\Vert }_{1,2}}(x,\hat{x}).$$

$${\tilde{x}}^{\mathrm{HD}}\in \underset{x;\hat{p}\in \partial \Vert \Gamma \cdot {\Vert }_{1,2}(x)}{\mathrm{argmin}}\frac{1}{2}\Vert \Phi x-y{\Vert }^2,$$

$${\tilde{x}}^{\mathrm{HO}}\in \underset{x;\pm \hat{p}\in \partial \Vert \Gamma \cdot {\Vert }_{1,2}(x)}{\mathrm{argmin}}\frac{1}{2}\Vert \Phi x-y{\Vert }^2.$$

$${\tilde{x}}_{\gamma }^{\mathrm{HD}}\in \underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert \Phi x-y{\Vert }^2+\gamma D_{\Vert \Gamma \cdot {\Vert }_{1,2}}(x,\hat{x}).$$

$${\tilde{x}}^{\mathrm{QO}}\in \underset{x;\mathrm{supp}(\Gamma x)\subseteq \hat{\mathcal{I}}}{\mathrm{argmin}}\frac{1}{2}\Vert \Phi x-y{\Vert }^2+\sum _{i\in \hat{\mathcal{I}}}\frac{\lambda }{2\Vert (\Gamma \hat{x}{)}_i\Vert }{\Vert (\Gamma x{)}_i-⟨(\Gamma x{)}_i,\frac{(\Gamma \hat{x}{)}_i}{\Vert (\Gamma \hat{x}{)}_i\Vert }⟩\frac{(\Gamma \hat{x}{)}_i}{\Vert (\Gamma \hat{x}{)}_i\Vert }\Vert }^2.$$

$${\tilde{x}}^{\varphi }\in \underset{x;\mathrm{supp}(\Gamma x)\subseteq \hat{\mathcal{I}}}{\mathrm{argmin}}\frac{1}{2}\Vert \Phi x-y{\Vert }^2+\sum _{i\in \hat{\mathcal{I}}}\varphi ((\Gamma x{)}_i,(\Gamma \hat{x}{)}_i),$$

$$\cos (z,\hat{z})=⟨\frac{z}{\Vert z\Vert },\frac{\hat{z}}{\Vert \hat{z}\Vert }⟩=\frac{1}{\Vert z\Vert \Vert \hat{z}\Vert }\sum _{j=1}^b{z}_j{\hat{z}}_j,$$

$$\text{and}{P}_{\hat{z}}(z)=⟨z,\frac{\hat{z}}{\Vert \hat{z}\Vert }⟩\frac{\hat{z}}{\Vert \hat{z}\Vert }=\frac{\Vert z\Vert }{\Vert \hat{z}\Vert }\cos (z,\hat{z})\hat{z},$$

$${\varphi }_{\mathrm{HD}}(z,\hat{z})={\iota }_{\{z\in {\mathbb{R}}^b:\cos (z,\hat{z})=1\}}$$

$${\varphi }_{\mathrm{HO}}(z,\hat{z})={\iota }_{\{z\in {\mathbb{R}}^b:|\cos (z,\hat{z})|=1\}}$$

$${\varphi }_{\mathrm{QO}}(z,\hat{z})=\frac{\lambda }{2\Vert \hat{z}\Vert }{\left|\left|z-P_{\hat{z}}(z)\right|\right|}^2=\frac{\lambda }{2}\frac{\Vert z{\Vert }^2}{\Vert \hat{z}\Vert }(1-{\cos }^2(z,\hat{z})).$$

$${\varphi }_{\mathrm{SD}}(z,\hat{z})=\lambda \Vert z\Vert (1-\cos (z,\hat{z})).$$

$${\tilde{x}}^{\mathrm{IB}}$$

$${\tilde{x}}^{\mathrm{HD}}$$

$${\tilde{x}}^{\mathrm{SD}}\in \underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert \Phi x-y{\Vert }^2+\lambda \sum _{i\in \hat{\mathcal{I}}}{D}_i(\Gamma x),$$

$$\underset{x}{\min }\underset{z}{\max }\frac{1}{2}\Vert \Phi x-y{\Vert }^2+⟨\Gamma x,z⟩-{\iota }_{B_2^{\lambda }}(z),$$

$${\omega }_{\varphi }(z,\hat{z},\hat{\mathcal{I}})=\sum _{i\in {\hat{\mathcal{I}}}^c}{\iota }_{\{0\}}(z_i)+\sum _{i\in \hat{\mathcal{I}}}\varphi (z_i,{\hat{z}}_i),\text{vspace}-0.1cm$$

$${\tilde{x}}^{\varphi }\in \underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert \Phi x-y{\Vert }^2+{\omega }_{\varphi }(\Gamma x,\Gamma \hat{x},\hat{\mathcal{I}}).\text{vspace}-0.1cm$$

$$\underset{x}{\min }\underset{z}{\max }\frac{1}{2}\Vert \Phi x-y{\Vert }^2+⟨\Gamma x,z⟩-{\omega }_{\varphi }^{\ast }(z,\Gamma \hat{x},\hat{\mathcal{I}}),\text{vspace}-0.1cm$$

$$\normalsize\begin{array}[]{@{}ll@{\hspace{.4em}}|@{\hspace{.4em}}ll@{}}\\[-13.99995pt] \hat{z}^{k+1}_{i}&=\frac{\hat{z}^{k}_{i}+\sigma(\Gamma\hat{v}^{k})_{i}}{\max(\lambda,\|\hat{z}^{k}_{i}+\sigma(\Gamma\hat{v}^{k})_{i}\|)}\hfil\hskip 3.99994pt\vrule\hskip 3.99994pt&\,\hat{\mathcal{I}}^{k}&=\left\{i\in[m]\;:\;\|\hat{z}^{k}_{i}+\sigma(\Gamma\hat{v}^{k})_{i}\|>\lambda\right\}\\ &\hfil\hskip 3.99994pt\vrule\hskip 3.99994pt&\,\tilde{z}^{k+1}&=\mathrm{prox}_{\sigma\omega_{\phi}^{*}(\cdot,\Psi(\hat{z}^{k},\hat{v}^{k}),\hat{\mathcal{I}}^{k})}(\tilde{z}^{k}+\sigma\Gamma\tilde{v}^{k})\\ \hat{x}^{k+1}&=\Phi_{\tau}^{-}\left(\hat{x}^{k}+\tau(\Phi^{\top}y-\Gamma^{\top}\hat{z}^{k+1})\right)\hfil\hskip 3.99994pt\vrule\hskip 3.99994pt&\,\tilde{x}^{k+1}&=\Phi_{\tau}^{-}\left(\tilde{x}^{k}+\tau(\Phi^{\top}y-\Gamma^{\top}\tilde{z}^{k+1})\right),\\ \hat{v}^{k+1}&=\hat{x}^{k+1}+\theta(\hat{x}^{k+1}-\hat{x}^{k})\hfil\hskip 3.99994pt\vrule\hskip 3.99994pt&\,\tilde{v}^{k+1}&=\tilde{x}^{k+1}+\theta(\tilde{x}^{k+1}-\tilde{x}^{k}),\vspace{-1em}\end{array}\vspace{1em}$$

$${\mathrm{prox}}_{\sigma {\omega }_{\varphi }^{\ast }}(z^0,\Gamma \hat{x},\hat{\mathcal{I}}{)}_i$$

$$\Psi ({\hat{z}}^k,{\hat{v}}^k{)}_i=\frac{\Vert {\hat{\nu }}_i^k\Vert -\lambda }{\sigma \Vert {\hat{\nu }}_i^k\Vert }{\hat{\nu }}_i^k\text{where}{\hat{\nu }}_i^k=({\hat{z}}^k+\sigma \Gamma {\hat{v}}^k{)}_i.$$

$${\varphi }^{\ast }({\alpha }_0,{\beta }_0)=\underset{\alpha ,\beta }{\sup }{\alpha }_0\alpha +⟨{\beta }_0,\beta ⟩-\varphi (\alpha ,\beta ).$$

Full text

1 Introduction

We consider linear inverse problems of the form $y=\Phi x+w$, where $y\in\mathbb{R}^{p}$ is an observed degraded image, $x\in\mathbb{R}^{n}$ the unknown clean image, $\Phi:\mathbb{R}^{n}\to\mathbb{R}^{p}$ a linear operator and $w\in\mathbb{R}^{p}$ a noise component, typically a zero-mean white Gaussian random vector with standard deviation $\sigma>0$. To reduce the effect of noise and the potential ill-conditioning of $\Phi$, we consider a regularized least square problem with an $\ell_{12}$ structured sparse analysis term of the form

[TABLE]

where $\lambda>0$ is a regularization parameter, $\Gamma:\mathbb{R}^{n}\to\mathbb{R}^{m\times b}$ is a linear analysis operator mapping an image over $m$ blocks of size $b$ and $\|z\|_{1,2}=\sum_{i=1}^{m}\|z_{i}\|=\sum_{i=1}^{m}(\sum_{j=1}^{b}z_{i,j}^{2})^{1/2}$, with $z_{i}=\{z_{i,j}\}_{j=1}^{b}\in\mathbb{R}^{b}$. This model is known to recover co-sparse solutions, *i.e., *such that $(\Gamma x)_{i}=0_{b}$ for most blocks $1\leq i\leq m$. A typical example is the one of isotropic total-variation (TViso) with $\Gamma=\nabla$ being the operator which extracts $m=n$ image gradient vectors of size $b=2$ (for volumes $b=3$, and so on). The anisotropic total-variation is another example corresponding to $\Gamma$ the operator which concatenates the vertical and horizontal components of the gradients into a vector of size $m=2n$, hence $b=1$.

1.1 Refitting

The co-support of an image $x$ (or support of $\Gamma x$) is the set of its non-zero blocks:

[TABLE]

While in some cases, the estimate $\hat{x}$ obtained by structured sparse analysis regularization (1) recovers correctly the co-support of the underlying signal $x$, it nevertheless suffers from a systematical bias in the estimated amplitudes $\hat{x}_{i}$. With TViso, this bias is reflected by a loss of contrast (see Fig. 1(b)). A standard strategy to reduce this effect, called refitting [8, 11, 1, 9], consists in approximating $y$ through $\Phi$ by an image sharing the same co-support as $\hat{x}$:

[TABLE]

where $\hat{\mathcal{I}}=\operatorname{supp}(\Gamma\hat{x})$. While this strategy works well for blocks of size $b=1$, it suffers from an excessive increase of variance whenever $b\geq 2$, *e.g., *for TViso. This is due to the fact that solutions do not only present sharp edges, but may involve gradual transitions. To cope with this issue, additional features of $\hat{x}$ than its co-support must be preserved by a refitting procedure. For the LASSO ($\Gamma=\mathrm{Id}$, $m=n$ and $b=1$), a pointwise preservation of the sign of $\hat{x}_{i}$ onto the support improves the numerical performances of the refitting [5].

1.2 Outline and contributions

In this paper, we introduce a new framework for refitting solutions promoted by $\ell_{12}$ structured sparse analysis regularization (1). In Section 2, we present related works, illustrated in Figure 1, that include Bregman iterations[10] or debiasing approaches [7, 2]. In Section 3, we describe our general variational refitting method for block penalties and show how the works presented in Section 2 can be described with such a framework. We discuss suitable properties a refitting block penalty should satisfy and introduce the Soft-penalized Direction model (SD), a new flexible refitting block penalty inheriting the advantages of the Bregman-based approaches. In Section 4, we propose a stable and one-step algorithm to compute our refitting strategy for any convex refitting block penalty, including the models in [2] and [7]. Experiments in Section 5 illustrate the practical benefits for the SD refitting and the potential of our framework for image processing.

2 Related re-fitting works

We first present some properties of Bregman divergences used all along the paper.

2.1 Properties of Bregman divergence of $\ell_{12}$ structured regularizers

In the literature [10, 2], Bregman divergences have proven to be well suited to measure the discrepancy between the biased solution $\hat{x}$ and its refitting $\tilde{x}$. We recall that for a convex function $\psi$, the associated (generalized) Bregman divergence between $x$ and $\hat{x}$ is, for any subgradient $\hat{p}\in\partial\psi(\hat{x})$: $D_{\psi}(x,\hat{x})=\psi(x)-\psi(\hat{x})-\langle\hat{p},x-\hat{x}\rangle\geq 0\enspace.$ If $\psi$ is an absolutely 1-homogeneous function ($\psi(\alpha x)=|\alpha|\psi(x)$, $\forall\alpha\in\mathbb{R}$) then $p\in\partial\psi(x)\Rightarrow\psi(x)=\left\langle p,\,x\right\rangle$ and the Bregman divergence simplifies into $D_{\psi}(x,\hat{x})=\psi(x)-\langle\hat{p},x\rangle\geq 0$. For regularizers of the form $\psi(x)=\|\Gamma x\|_{1,2}$, a subgradient $\hat{p}\in\partial\|\Gamma\cdot\|_{1,2}(\hat{x})$ satisfies $\hat{p}=\Gamma^{\top}\hat{z}$ with [3]:

[TABLE]

where we have $\Gamma^{\top}\hat{z}\in\partial\|\Gamma\cdot\|_{1,2}(\hat{x})\Leftrightarrow\hat{z}_{i}\in\partial\|\cdot\|((\Gamma\hat{x})_{i}),\,\forall i\in[m].$ Now denoting

[TABLE]

there exists an interesting relation between the global Bregman divergence on $x$’s and the local ones on $(\Gamma x)_{i}$’s:

[TABLE]

Combining relations (4) and (5), we see that on the co-support $i\in\hat{\mathcal{I}}$, such divergence measures the fit of directions between $(\Gamma x)_{i}$ and $(\Gamma\hat{x})_{i}$:

[TABLE]

This divergence also partially captures the co-support as we have from (5)

[TABLE]

2.2 Bregman-based Refitting

We now review some refitting methods based on Bregman divergences.

2.2.1 Flexible Iterative Bregman regularization

The Bregman process [10] reduces the bias of solutions of (1) by successively solving problems of the form

[TABLE]

We here consider a fixed $\lambda$, but different refitting strategies can be considered with decreasing parameters $\lambda_{l}$ as in [12, 13].

Setting $\tilde{x}_{0}=0_{n}$, we have $0_{n}\in\partial\|\Gamma\cdot\|_{1,2}(0_{n})$ so that $D_{\|\Gamma\cdot\|_{1,2}}(x,0_{n})=\|\Gamma x\|_{1,2}$ and we recover the biased solution of (1) as $\tilde{x}_{1}=\hat{x}$. We denote by $\tilde{x}^{\mathrm{IB}}=\tilde{x}_{2}$ the refitting obtained after $2$ steps of the Iterative Bregman (IB) procedure (9):

[TABLE]

As underlined in relation (7), by minimizing $D_{\|\Gamma\cdot\|_{1,2}}(x,\hat{x})$, one aims at preserving the direction of $\Gamma\hat{x}$ on the support $\hat{\mathcal{I}}$, without ensuring $\operatorname{supp}(\Gamma\tilde{x}^{\textrm{IB}})\subseteq\hat{\mathcal{I}}$. This issue can be observed in the background of Cameraman in Fig. 1(c), where noise reinjection is visible. For the iterative framework, the co-support of the previous solution may indeed not be preserved ($\|(\Gamma\tilde{x}_{l})_{i}\|=0\nRightarrow\|(\Gamma\tilde{x}_{l+1})_{i}\|=0$) and can hence grow. The support of $\Gamma x_{0}$ for $\tilde{x}_{0}=0_{n}$ is for instance totally empty whereas the one of $\hat{x}=\tilde{x}_{1}$ may not (and should not) be empty. For $l\to\infty$, the process actually converges to some $x$ such that $\Phi x=y$.

2.2.2 Hard-constrained refitting without explicit support identification.

In order to respect the support of the biased solution $\hat{x}$ and to keep track of the direction $\Gamma\hat{x}$ during the refitting, the authors of [2] proposed the following model:

[TABLE]

for $\hat{p}\in\partial\|\Gamma\cdot\|_{1,2}(\hat{x})$. This model enforces the Bregman divergence to be [math], since $\hat{p}\in\partial\|\Gamma\cdot\|_{1,2}(x)\Rightarrow\|\Gamma x\|_{1,2}=\langle\hat{p},x\rangle\Rightarrow D_{\|\Gamma\cdot\|_{1,2}}(x,\hat{x})=0$.

We see from (7) that for $i\in\hat{\mathcal{I}}$, the direction of $(\Gamma\hat{x})_{i}$ is preserved in the refitted solution. Following relations (4) and (8), the co-support is also preserved for any $i\in\hat{\mathcal{I}}^{c}$ such that $\|\hat{z}_{i}\|<1$. Note though that extra elements in the co-support $\Gamma\tilde{x}^{\textrm{HD}}$ may be added at coordinates $i\in\hat{\mathcal{I}}^{c}$ such that $\|\hat{z}_{i}\|=1$. We denote this model as HD, for Hard-constrained Direction. To get ride of the direction dependency, a Hard-constrained Orientation (HO) model is also proposed in [2]:

[TABLE]

The orientation model may nevertheless involve contrast inversions between biased and refitted solutions, as shown in Fig. 1(d) with the banana dark shape in the white region. In practice, relaxations are used in [2] by solving, for a large value $\gamma>0$

[TABLE]

The main advantage of this refitting strategy is that no support identification is required since everything is implicitly encoded in the subgradient $\hat{p}$. This makes the process stable even if the estimation of $\hat{x}$ is not highly accurate. The support of $\Gamma\hat{x}$ is nevertheless only approximately preserved, since the constraint $D_{\|\Gamma\cdot\|_{1,2}}(x,\hat{x})=0$ can never be ensured numerically with a finite value of $\gamma$. Finally, as shown in Fig. 1(d-e), such constrained approaches lack of flexibility since the orientation of $\Gamma\tilde{x}$ cannot deviate from the one of $\Gamma\hat{x}$ (for complex signals, such as Cameraman, amplitudes remain significantly biased and less details are recovered).

2.3 Flexible Quadratic refitting without support identification

We now describe an alternative way for performing variational refitting. When specialized to $\ell_{1,2}$ sparse analysis regularization, CLEAR, a general refitting framework [7], consists in computing

[TABLE]

This model promotes refitted solutions preserving to some extent the orientation $\Gamma\hat{x}$ of the biased solution. It also shrinks the amplitude of $\Gamma x$ all the more that the amplitude of $\Gamma\hat{x}$ are small. As this model penalizes changes of orientation, we refer to it as QO for Quadratic-penalized Orientation. This penalty does not promote any kind of direction preservation, and as for the HO model, contrast inversions may be observed between biased and refitted solutions (see Fig. 1(f)). The quadratic term also over-penalizes large changes of orientation.

3 Refitting with block penalties

As mentioned in the previous section, the methods of [2] and [7] have proposed variational refitting formulations that not only aim at preserving the co-support $\hat{\mathcal{I}}$ but also the orientation of $(\Gamma\hat{x})_{i\in\hat{\mathcal{I}}}$. In this paper, we propose to express these (two-steps) refitting procedures in the following general framework

[TABLE]

where $\phi:\mathbb{R}^{b}\times\mathbb{R}^{b}\to\mathbb{R}$ is a refitting block penalty ($b\geq 1$ is the size of the blocks) promoting $\Gamma x$ to share information with $\Gamma\hat{x}$ in some sense to be specified. To refer to some features of the vector $\Gamma\hat{x}$, let us first define properly the notions of relative orientation, direction and projection between two vectors.

Definition 1

Let $z$ and $\hat{z}$ being two vectors in $\mathbb{R}^{b}$, we define

[TABLE]

where $P_{\hat{z}}(z)$ is the orthogonal projection of $z$ onto $\operatorname{Span}(\hat{z})$ (i.e., the orientation axis of $\hat{z}$). We say that $z$ and $\hat{z}$ share the same orientation (resp. direction), if $|\cos(z,\hat{z})|=1$ (resp. $\cos(z,\hat{z})=1$).

Thanks to Definition 1, we can now reformulate the previous refitting models in terms of block penalties. The Hard-constrained refitting models preserving Direction (11) and Orientation (12) of [2] as well as the flexible Quadratic model of CLEAR [7] correspond to the following block penalties:

[TABLE]

where $\iota_{\mathcal{C}}$ is the $0/+\infty$ indicator function of a set $\mathcal{C}$. These block penalties are either insensitive to directions (HO and QO) or intolerant to small changes of orientations (HD and HO), hence not satisfying (*cf. *drawbacks visible in Fig. 1).

When $b=1$, the orientation-based penalties (QO and HO) have absolutely no effect while the direction-based penalty HD preserves the sign of $(\Gamma\hat{x})_{i}$. In this paper, we argue that the direction of $(\Gamma\hat{x})_{i}$, for any $b\geq 1$, carries important information that is worth preserving when refitting, at least to some extent.

3.1 Desired properties of refitting block penalties

To compute global optimum of the refitting model (15), we only consider convex refitting block penalties $z\mapsto\phi(z,\hat{z})$. Hence, we now introduce properties a block penalty $\phi$ should satisfy for refitting purposes:

• (P1)

$\phi$ is convex, non negative and $\phi(z,\hat{z})=0$, if $\cos(z,\hat{z})=1$ or $\|z\|=0$,

• (P2)

$\phi(z^{\prime},\hat{z})\geq\phi(z,\hat{z})$ if $\|z^{\prime}\|=\|z\|$ and $\cos(z,\hat{z})\geq\cos(z^{\prime},\hat{z})$,

• (P3)

$z\mapsto\phi(z,\hat{z})$ is continuous,

Property (P1) stipulates that no configuration can be more favorable than $z$ and $\hat{z}$ having the same direction. Hence, the direction of the refitted solution should be encouraged to follow the one of the biased solution. Property (P2) imposes that for a fixed amplitude, the penalty should be increasing w.r.t. the angle formed with $\hat{z}$. Property (P3) enforces refitting that can continuously adapt to the data and be robust to small perturbations.

3.2 A new flexible refitting block penalty

We now introduce our refitting block penalty designed to preserve the desired features of $\hat{z}=\Gamma\hat{x}$ in a simple way. The Soft-penalized Direction penalty reads

[TABLE]

The properties of the different studied block penalties are presented in Table 1. The proposed SD model is the only one satisfying all the desired properties. As illustrated in Fig. 2, it is a continuous penalization that increases continuously with respect to the absolute angle between $z$ and $\hat{z}$.

3.3 SD block penalty: the best of both Bregman worlds

Denoting $D_{i}(\Gamma x)=D_{\|\cdot\|}((\Gamma x)_{i},(\Gamma\hat{x})_{i})$ as introduced in Section 2.1, the refitting models given in (10) and (11) can be expressed as

[TABLE]

Observing that the new SD block penalty (21) may be rewritten as $\phi_{\mathrm{SD}}(z,\hat{z})=\lambda\left(\|z\|-\left\langle z,\,{\hat{z}}/{\|\hat{z}\|}\right\rangle\right)=\lambda D_{\|\cdot\|}(z,\hat{z}),$ the SD refitting model is

[TABLE]

With such reformulations, connections between refitting models (22), (23) and (24) can be clarified. The solution $\tilde{x}^{\mathrm{IB}}$ [10] is too relaxed, as it only penalizes the directions $(\Gamma x)_{i}$ using $(\Gamma\hat{x})_{i}$, without aiming at preserving the co-support of $\hat{x}$. The solution $\tilde{x}^{\mathrm{HD}}$ [2] is too constrained: the direction within the co-support is required to be preserved exactly. Our proposed refitting $\tilde{x}^{\mathrm{SD}}$ lies in-between: it preserves the co-support, while authorizing some directional flexibility, as illustrated by the sharper square edges in Figure 1(g).

With respect to the Hard-constrained approach [2], an important difference is that we consider local inclusions of subgradients of the function $\lambda\|\cdot\|_{1,2}$ at point $(\Gamma x)_{i}$ instead of the global inclusion of subgradients of the function $\lambda\|\Gamma\cdot\|_{1,2}$ at point $x$ as in HD (11) and HO (12). Such a change of paradigm allows to adapt the refitting locally by preserving the co-support while including the flexibility of the original Bregman approach [10].

4 Refitting in practice

We first describe how computing $\hat{x}$ and its refitting $\tilde{x}$ in two successive steps and then propose a new numerical scheme for the joint computation of $\hat{x}$ and $\tilde{x}$.

4.1 Biased problem and posterior refitting

To obtain $\hat{x}$ solution of (1), we consider the primal dual formulation that reads

[TABLE]

where $\iota_{B_{2}^{\lambda}}$ is the indicator function of the $\ell_{2}$ ball of radius $\lambda$ that is [math] if $\|z_{i}\|\leq\lambda$ for all $i\in[m]$ and $+\infty$ otherwise. This problem can be solved with the iterative primal-dual algorithm [4] presented in the left part of Algorithm (29). For positive parameters satisfying $\tau\sigma<\|\Gamma\|$ and $\theta\in[0,1]$, the iterates $(\hat{z}^{k},\hat{x}^{k})$ converge to a saddle point $(\hat{z},\hat{x})$ satisfying $\hat{z}_{i}\in\partial\|\cdot\|((\Gamma\hat{x})_{i})$, $\forall i$.

Assume that the co-support $\hat{\mathcal{I}}=\operatorname{supp}(\Gamma\hat{x})$ of the biased solution has been identified, a posterior refitting can be obtained by solving (15) for any refitting block penalty $\phi$. To that end, we write the characteristic function of co-support preservation as $\sum_{i\in\hat{\mathcal{I}}^{c}}\iota_{\{0\}}(z_{i}),$ where $\iota_{\{0\}}(z)=0$ if $z=0$ and $+\infty$ otherwise. By introducing the convex function

[TABLE]

the general refitting problem (15) can be expressed as

[TABLE]

Subsequently, we can consider its primal dual formulation

[TABLE]

where $\omega_{\phi}^{*}(z,\Gamma\hat{x},\hat{\mathcal{I}})=\sup_{x}\langle z,x\rangle-\omega_{\phi}(x,\Gamma\hat{x},\hat{\mathcal{I}})$ is the convex conjugate, with respect to the first argument, of $\omega_{\phi}(\cdot,\Gamma\hat{x},\hat{\mathcal{I}})$. Such problem can again be solved with the primal-dual algorithm [4]. The crucial point is to have an accurate estimation of the vector $\Gamma\hat{x}$ and its support $\hat{\mathcal{I}}$. Yet, it is well known that estimating $\operatorname{supp}(\Gamma\hat{x})$ from an estimation $\hat{x}^{k}$ is not stable numerically: the support $\operatorname{supp}(\Gamma\hat{x}^{k})$ can be far from $\operatorname{supp}(\Gamma\hat{x})$ even though $\hat{x}^{k}$ is arbitrarily close to $\hat{x}$.

4.2 Joint-refitting algorithm

We now introduce a general algorithm aiming to jointly solve the original problem (1) and the refitting one (15) for any refitting block penalty $\phi$. This framework has been developed for stable projection onto the support in [6] and later extended to refitting with the Quadratic Orientation penalty in [7]. The strategy consists in solving in parallel the two problems (25) and (28).

Two iterative primal-dual algorithms are used for the biased variables $(\hat{z}^{k},\hat{x}^{k})$ and the refitted ones $(\tilde{z}^{k},\tilde{x}^{k})$. Let us now present the whole algorithm:

[TABLE]

with the operator $\Phi_{\tau}^{-}=(\mathrm{Id}+\tau\Phi^{\top}\Phi)^{-1}$ and the auxiliary variables $\hat{v}^{k}$ and $\tilde{v}^{k}$ that accelerate the algorithm. Following [4], for any positive scalars $\tau$ and $\sigma$ satisfying $\tau\sigma\|\Gamma^{\top}\Gamma\|<1$ and $\theta\in[0,1]$, the estimates $(\hat{z}^{k},\hat{x}^{k},\hat{v}^{k})$ of the biased solution converge to $(\hat{z},\hat{x},\hat{x})$, where $(\hat{z},\hat{x})$ is a saddle point of (25).

In the right part of Algorithm (29), we rely on the proximal operator that is, for a convex function $\psi$ and at point $z_{0}$, $\mathrm{prox}_{\sigma\psi}(z)=\operatorname{argmin}_{z}\tfrac{1}{2\sigma}\|z-z_{0}\|^{2}+\psi(z)$. From the block structure of the function $\omega_{\phi}$ defined in (26), the computation of its proximal operator may be realized pointwise. Since $\iota_{\{0\}}(z)^{*}=0$, we have

[TABLE]

Table 2 gives the expressions of the dual functions $\phi^{*}$ with respect to their first variable and their related proximal operators $\mathrm{prox}_{\sigma\phi^{*}}$ for the refitting block penalties considered in this paper. All details are given in the Appendix.

The idea behind this joint-refitting algorithm is to perform online co-support detection using the dual variable $\hat{z}^{k}$ of the biased variable $\hat{x}^{k}$. From relations in (4), we expect at convergence $\hat{z}^{k}$ to saturate on the support of $\Gamma\hat{x}$ and to satisfy the optimality condition $\hat{z}_{i}^{k}=\lambda\tfrac{(\Gamma\hat{x})_{i}}{\|(\Gamma\hat{x})_{i}\|}$. In practice, the norm of the dual variable $\hat{z}^{k}_{i}$ saturates to $\lambda$ relatively fast onto $\hat{\mathcal{I}}$. As a consequence, it is far more stable to detect the support of $\Gamma\hat{x}$ with the dual variable $\hat{z}^{k}$ than with the vector $\Gamma\hat{x}^{k}$ itself. In the first step of Algorithm (29), the condition $\|\hat{z}^{k}_{i}+\sigma(\Gamma\hat{v}^{k})_{i}\|>\lambda$ is thus used to detect elements of the support $i\in\hat{\mathcal{I}}^{k}$ of $\Gamma\hat{x}^{k}$ along iterations***As in [2], extended support $\|\hat{z}_{i}\|=\lambda$ can be tackled by testing $\|(\hat{z}^{k}+\sigma\Gamma\hat{v}^{k})_{i}\|\geq\lambda$.. The function $\Psi(\hat{z}^{k},\hat{v}^{k})$ aims at approximating $\Gamma\hat{x}$ with the current values of the available variables $(\hat{z}^{k},\hat{v}^{k})$ of the biased problem that is solved simultaneously. Following [7], the function $\Psi$ can be chosen as

[TABLE]

that satisfies $(\Psi(\hat{z},\hat{v}))_{i}=(\Gamma\hat{x})_{i}$ at convergence, while appearing to give very stable online estimations of the direction of $\Gamma\hat{x}$ through $\hat{z}^{k}$.

This joint-estimation considers at every iteration $k$ different refitting functions $\omega_{\phi}^{*}(.,\Psi(\hat{z}^{k},\hat{v}^{k}),\hat{\mathcal{I}}^{k})$ in (27). For fixed values of $\Psi(\hat{z}^{k},\hat{v}^{k})$ and $\hat{\mathcal{I}}^{k}$, the refitted variables $(\tilde{z}^{k},\tilde{x}^{k})$ in the right part of the Algorithm (29) converges since it exactly corresponds to the primal-dual algorithm [4] applied to the problem (28). However, unless $b=1$ (see [6]), we do not have guarantee of convergence of the presented scheme with a varying $\omega_{\phi}^{*}$. As in [7], we nevertheless observe convergence and a very stable behavior for this algorithm.

In addition to its better numerical stability, the running time of joint-refitting is more interesting than the posterior approach. In Algorithm (29), the refitting variables at iteration $k$ require the biased variables at the same iteration and the whole process can be realized in parallel without significantly affecting the running time of the original biased process. On the other hand, posterior refitting is necessarily sequential and the running time is doubled in general.

5 Results

We considered TViso regularization of degraded color images. We defined blocks obtained by applying $\Gamma=[\nabla_{x}^{R},\nabla_{y}^{R},\nabla_{x}^{G},\nabla_{y}^{G},\nabla_{x}^{B},\nabla_{y}^{B}]$ where $m=n$, $b=6$, and $\nabla^{C}_{d}$ denotes the gradient in the direction $d\in\{x,y\}$ for the color channel $C\in\{R,G,B\}$. We first focused on a denoising problem $y=x+w$ where $x$ is an 8bit color image and $w$ is an additive white Gaussian noise with standard deviation $\sigma=20$. We next focused on a deblurring problem $y=\Phi x+w$ where $x$ is an 8bit color image, $\Phi$ is a convolution simulating a directional blur, and $w$ is an additive white Gaussian noise with standard deviation $\sigma=2$. We chose $\lambda=4.3\sigma$. We applied the iterative primal-dual algorithm with our joint-refitting (Algorithm (29)) for $1,000$ iterations, with $\tau=1/4$, $\sigma=1/6$ and $\theta=1$.

Results are provided on Fig. 3 and 4. Comparisons of refitting with our proposed SD block penalty, HD (only for denoising) and QO are provided. Using our proposed SD block penalty offers the best refitting performances in terms of both visual and quantitative measures. The loss of contrast of TViso is well-corrected, amplitudes are enhanced while smoothness and sharpness of TViso is preserved. Meanwhile, the approach does not create artifacts, invert contrasts, or reintroduce information that were not recovered by TViso.

6 Conclusion

In this work, we have reformulated the refitting problem of solutions promoted by $\ell_{12}$ structured sparse analysis in terms of block penalties. We have introduced a new block penalty that interpolates between Bregman iterations [10] and direction preservation [2].This framework easily allows the inclusion of additional desirable properties of refitted solutions as well as new penalties that may increase refitting performances. In order to take advantage of our efficient joint-refitting algorithm, it is important to consider simple block penalty functions, which proximal operator can be computed explicitly or at least easily. Refitting in the case of other regularizers and loss functions will be investigated in the future.

Acknowledgments This project has been carried out with support from the French State, managed by the French National Research Agency (ANR-16-CE33-0010-01). This work was supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 777826.

Appendix A Proximity operators of block penalties

A.1 Convex conjugates $\phi^{*}$

We here compute the convex conjugate $\phi^{*}$ of the different block penalties $\phi(z,\hat{z})$ that only depends on $z\in\mathbb{R}^{b}$ and where $\hat{z}\in\mathbb{R}^{b}$ is a given fixed non null vector. We consider the following representation of the vectors $z$ with respect to the $\hat{z}$ axis: $z=\alpha\frac{\hat{z}}{\|\hat{z}\|}+\beta\frac{\hat{z}^{\perp}}{\|\hat{z}\|}$. This expression is valid for the case $b=2$. If $b=1$ then $z$ is only parameterized by $\alpha$. When $b>2$, $\hat{z}^{\perp}$ must be understood as a subspace $S$ of dimension $b-1$ and $\beta$ as a vector of $b-1$ components corresponding to each dimension of $S$. With this change of variables, we have $\|z\|^{2}=\alpha^{2}+\|\beta\|^{2}$ (since $\beta$ is of dimension $b-1$), $\cos(z,\hat{z})=\alpha/\sqrt{\alpha^{2}+\|\beta\|^{2}}$, $P_{\hat{z}}(z)=\alpha\hat{z}/\|\hat{z}\|$ and $z-P_{\hat{z}}(z)=\beta\hat{z}^{\perp}/\|\hat{z}\|$. We also observe for instance that $|\cos(z,\hat{z})|=1\Leftrightarrow\|\beta\|=0$. All the block penalties $\phi(z,\hat{z})$ can thus be expressed as $\phi(\alpha,\beta)$. The convex conjugate reads

[TABLE]

The block penalty $\phi_{\mathrm{HO}}$

reads

[TABLE]

Hence, it gives

[TABLE]

The block penalty $\phi_{\mathrm{HD}}$

reads

[TABLE]

Hence, it gives

[TABLE]

The block penalty $\phi_{\mathrm{QO}}$

reads

[TABLE]

Hence, it gives

[TABLE]

The optimality condition on $\beta$ give

[TABLE]

so that

[TABLE]

The block penalty $\phi_{\mathrm{QD}}$

reads

[TABLE]

Hence, it gives

[TABLE]

We observe that if $\alpha_{0}>0$, taking $\beta=0$ and letting $\alpha\to\infty$ leads to $\phi^{*}(\alpha_{0},\beta_{0})=+\infty$. Next, the optimality conditions on $\beta$ and $\alpha$ give

[TABLE]

so that

[TABLE]

The block penalty $\phi_{\mathrm{SD}}$

reads $\phi(z,\hat{z})=\frac{\lambda}{2}\|z\|(1-\cos(z,\hat{z})).$ Hence, it gives

[TABLE]

We observe that if $\sqrt{\|\beta_{0}\|^{2}+(\alpha_{0}+\lambda/2)^{2}}>\lambda/2$, then letting $\alpha\to\operatorname{sign}{(\alpha_{0}+\lambda/2)}\times\infty$ and $\beta\to\operatorname{sign}{\beta_{0}}\times\infty$ leads to $\phi^{*}(\alpha_{0},\beta_{0})=+\infty$. As a consequence we find

[TABLE]

A.2 Computing $\mathrm{prox}_{\sigma\phi^{*}}$

We here give the computation of the proximal operator of the different $\phi^{*}$ that is given at point $(\alpha_{0},\beta_{0})$ by

[TABLE]

Block penalty $\phi_{\mathrm{HO}}$.

We have

[TABLE]

Block penalty $\phi_{\mathrm{HD}}$.

We have

[TABLE]

Block penalty $\phi_{\mathrm{QO}}$.

We have

[TABLE]

since the optimality condition with respect to $\beta$ gives $\lambda(\beta-\beta_{0})+\sigma\|\hat{z}\|\beta=0$.

Block penalty $\phi_{\mathrm{QD}}$.

We have

[TABLE]

Block penalty $\phi_{\mathrm{SD}}$.

We have

[TABLE]

which just corresponds to the projection of the $\ell_{2}$ ball of $\mathbb{R}^{b}$ of radius $\lambda/2$ and center $(-\lambda/2,0)$.