A Finite Element Nonoverlapping Domain Decomposition Method with Lagrange Multipliers for the Dual Total Variation Minimizations

TL;DR

This paper introduces a primal-dual domain decomposition method using Lagrange multipliers for total variation regularized image processing, demonstrating superior performance in denoising, inpainting, and segmentation tasks.

Contribution

It develops a novel finite element domain decomposition approach with Lagrange multipliers for dual total variation minimization, improving computational efficiency and accuracy.

Findings

Outperforms existing state-of-the-art methods in image denoising, inpainting, and segmentation.

Effective handling of interface constraints via Lagrange multipliers.

Applicable to problems with $L^2$- and $L^1$-fidelity.

Abstract

In this paper, we consider a primal-dual domain decomposition method for total variation regularized problems appearing in mathematical image processing. The model problem is transformed into an equivalent constrained minimization problem by tearing-and-interconnecting domain decomposition. Then, the continuity constraints on the subdomain interfaces are treated by introducing Lagrange multipliers. The resulting saddle point problem is solved by the first order primal-dual algorithm. We apply the proposed method to image denoising, inpainting, and segmentation problems with either $L^2$-fidelity or $L^1$-fidelity. Numerical results show that the proposed method outperforms the existing state-of-the-art methods.

Tables

Table 1. Table 1: Performance of the proposed method for the image denoising problem with the L 1 superscript 𝐿 1 L^{1} fidelity term ( 5.2 )

1 1 Test image	1 1 $𝒩$	ALG1				LNP
		PSNR	iter	max inner iter	wall-clock time (sec)	PSNR	iter	max inner iter	wall-clock time (sec)
1 Cameraman $2048\times 2048$	1	47.68	-	546	86.30
	$2\times 2$	48.26	39	266	155.85	46.75	33	1761	587.86
	$4\times 4$	48.26	45	277	50.40	46.75	61	1679	141.31
	$8\times 8$	48.25	49	274	14.79	46.78	134	1796	42.37
	$16\times 16$	48.24	54	279	2.19	46.79	175	1200	7.82
1 Boat $2048\times 3072$	1	33.82	-	2093	476.02
	$2\times 2$	33.92	61	283	531.97	33.60	43	1815	970.12
	$4\times 4$	33.92	63	276	165.45	33.62	71	1823	255.70
	$8\times 8$	33.92	66	258	53.68	33.62	114	1413	87.23
	$16\times 16$	33.92	71	246	8.18	33.63	161	844	14.07

Table 2. Table 2: Performance of the proposed method for the image inpainting problem with the L 2 superscript 𝐿 2 L^{2} fidelity term ( 5.5 )

Test image	$𝒩$	PSNR	iter	max inner iter	wall-clock time (sec)
1 Cameraman $2048\times 2048$	1	24.96	-	2324	389.88
	$2\times 2$	25.00	438	613	805.88
	$4\times 4$	25.00	440	713	260.06
	$8\times 8$	25.00	445	99	88.96
	$16\times 16$	25.00	460	99	13.40
1 Boat $2048\times 3072$	1	24.32	-	1861	436.32
	$2\times 2$	24.33	285	817	1085.96
	$4\times 4$	24.33	286	588	361.04
	$8\times 8$	24.33	289	111	132.14
	$16\times 16$	24.33	298	117	25.63

Table 3. Table 3: Performance of the proposed method for the image inpainting problem with the L 1 superscript 𝐿 1 L^{1} fidelity term ( 5.6 )

1 1 Test image	1 1 $𝒩$	ALG1				LNP
		PSNR	iter	max inner iter	wall-clock time (sec)	PSNR	iter	max inner iter	wall-clock time (sec)
1 Cameraman $2048\times 2048$	1	34.70	-	1693	304.05
	$2\times 2$	35.56	78	921	615.5	34.18	173	1530	1132.71
	$4\times 4$	35.53	86	958	244.21	34.24	190	1672	323.98
	$8\times 8$	35.53	90	274	68.44	34.26	205	1796	97.65
	$16\times 16$	35.53	101	279	16.10	34.37	239	1200	24.67
1 Boat $2048\times 3072$	1	30.83	-	2262	604.34
	$2\times 2$	31.08	80	860	1283.93	30.53	368	1703	2134.45
	$4\times 4$	31.08	84	765	432.68	30.55	374	1888	635.88
	$8\times 8$	31.08	88	258	128.09	30.56	384	1413	217.33
	$16\times 16$	31.07	97	246	31.95	30.58	402	844	60.68

Table 4. Table 4: Performance of the proposed method for the image segmentation problem ( 5.7 )

1 1 Test image	1 1 $𝒩$	ALG1			DCT
		iter	max inner iter	wall-clock time (sec)	iter	max inner iter	wall-clock time (sec)
1 Cameraman $2048\times 2048$	1	-	1070	172.47
	$2\times 2$	16	258	153.72	14	290	85.67
	$4\times 4$	16	63	37.58	15	52	18.21
	$8\times 8$	16	55	6.65	15	52	5.31
	$16\times 16$	16	56	1.10	16	55	1.43

Equations

$$\underset{u\in BV(\Omega )}{\min }\left\{\frac{\alpha }{2}{\int }_{\Omega }(u-f{)}^{2}dx+TV(u)\right\},$$

$$TV(u)=\sup \left\{{\int }_{\Omega }u\mathrm{div}\mathbf{q}dx:\mathbf{q}\in (C_0^1(\Omega ){)}^2,|\mathbf{q}(x)|\le 1\forall x\in \Omega \right\},$$

$$\underset{\mathbf{p}\in H_0(\mathrm{div};\Omega )}{\min }\frac{1}{2\alpha }{\int }_{\Omega }(\mathrm{div}\mathbf{p}+\alpha f{)}^{2}dx\text{subject to }|\mathbf{p}(x)|\le 1\forall x\in \Omega$$

$$\underset{u\in BV(\Omega )}{\min }\left\{\alpha F(u)+TV(u)\right\},$$

$$\underset{\mathbf{p}\in H_0(\mathrm{div};\Omega )}{\min }\frac{1}{\alpha }{F}^{\ast }(\mathrm{div}\mathbf{p})\text{subject to }|\mathbf{p}(x)|\le 1\forall x\in \Omega ,$$

$$\Vert v{\Vert }_{p,\mathcal{H}}={\left(\sum _{i=1}^n|v_i{|}^p\right)}^{\frac{1}{p}},$$

$${⟨v,w⟩}_{\mathcal{H}}=\sum _{i=1}^n{v}_i{w}_i.$$

$$F^{\ast }(v^{\ast })=\underset{v\in \mathcal{H}}{\sup }\left\{⟨v,v^{\ast }⟩-F(v)\right\}$$

$$⟨v,Aw⟩=⟨A^{\ast }v,w⟩$$

$$\mathrm{dom}F=\left\{u\in \mathcal{H}:F(u)<\infty \right\}.$$

$${\chi }_C(v)=\{\begin{array}{ll}0& \text{if }v\in C,\\ \infty & \text{if }v\notin C.\end{array}$$

$$\underset{u\in {\mathcal{H}}_1}{\min }\left\{F(u)+G(Ku)\right\},$$

$$\underset{u\in {\mathcal{H}}_1}{\min }\left\{F(u)+G(Ku)\right\}$$

$$=-\left[\underset{q\in {\mathcal{H}}_2}{\min }\left\{F^{\ast }(-K^{\ast }q)+G^{\ast }(q)\right\}\right].$$

$$X=\left\{u\in L^2(\Omega ):u{|}_T\text{ is constant}\forall T\in \mathcal{T}\right\}.$$

$${\varphi }_T(x)=\{\begin{array}{ll}1& \text{ if }x\in T,\\ 0& \text{ if }x\notin T,\end{array}T\in \mathcal{T}$$

$$u=\sum _{T\in \mathcal{T}}(u{)}_T{\varphi }_T.$$

$$Y=\left\{\mathbf{q}\in H_0(\mathrm{div};\Omega ):\mathbf{q}{|}_T\in {\mathcal{R}\mathcal{T}}_0(T)\forall T\in \mathcal{T}\right\},$$

$$\mathbf{q}(x_1,x_2)=\left[\begin{array}{c}{a}_1+b_1{x}_1\\ a_2+b_2{x}_2\end{array}\right].$$

$$\mathbf{p}=\sum _{i\in \mathcal{I}}(\mathbf{p}{)}_i{\mathbf{\psi }}_i.$$

$$C=\left\{\mathbf{p}\in Y:|(\mathbf{p}{)}_i|\le 1\forall i\in \mathcal{I}\right\}.$$

$$({\mathrm{proj}}_C\mathbf{p}{)}_i=\frac{(\mathbf{p}{)}_i}{\max \left\{1,|(\mathbf{p}{)}_i|\right\}}\forall i\in \mathcal{I}.$$

$$X_s=\left\{u\in L^2({\Omega }_s):u{|}_T\text{ is constant}\forall T\in {\mathcal{T}}_s\right\}.$$

$$R_{s}u=u{|}_{{\Omega }_s}.$$

$$(R_s^{\ast }{u}_s{)}_T=\{\begin{array}{ll}(u_s{)}_T& \text{ if }T\subset {\Omega }_s,\\ 0& \text{ if }T\not\subset {\Omega }_s,\end{array}T\in \mathcal{T}.$$

$${\tilde{Y}}_s=\left\{{\tilde{\mathbf{q}}}_s\in H(\mathrm{div};{\Omega }_s):{\tilde{\mathbf{q}}}_s\cdot {\mathbf{n}}_s=0\text{ on }\partial {\Omega }_s\setminus \Gamma \text{, }{\tilde{\mathbf{q}}}_s{|}_T\in {\mathcal{R}\mathcal{T}}_0(T)\forall T\in {\mathcal{T}}_s\right\}.$$

$${\tilde{C}}_s=\left\{{\tilde{\mathbf{p}}}_s\in {\tilde{Y}}_s:|({\tilde{\mathbf{p}}}_s{)}_i|\le 1\forall i\in {\tilde{\mathcal{I}}}_s\right\}.$$

$$({\mathrm{proj}}_{{\tilde{C}}_s}{\tilde{\mathbf{p}}}_s{)}_i=\frac{({\tilde{\mathbf{p}}}_s{)}_i}{\max \left\{1,|({\tilde{\mathbf{p}}}_s{)}_i|\right\}}\forall i\in {\tilde{\mathcal{I}}}_s.$$

$$B\tilde{\mathbf{p}}{|}_{{\Gamma }_{st}}={\tilde{\mathbf{p}}}_s\cdot {\mathbf{n}}_{st}-{\tilde{\mathbf{p}}}_t\cdot {\mathbf{n}}_{st},s<t,$$

$$\underset{u\in X}{\min }\left\{\mathcal{J}(u):=\alpha F(u)+TV(u)\right\},$$

$$TV(u)=\underset{\mathbf{p}\in C}{\max }{⟨u,\mathrm{div}\mathbf{p}⟩}_X,$$

Full text

††journal: arXiv

1 Introduction

After a pioneering work of Rudin et al. ROF:1992 , total variation minimization has been widely used in image processing. In their work, authors proposed an image denoising model with the total variation regularizer, which is called Rudin–Osher–Fatemi (ROF) model, as follows:

[TABLE]

where $\Omega$ is an image domain, $f$ is a corrupted image, $\alpha>0$ is a positive denoising parameter, and $TV(u)$ is the total variation of $u$ defined as

[TABLE]

where $|\cdot|$ denotes the Euclidean norm in $\mathbb{R}^{2}$. The solution space $BV(\Omega)$ is the collection of $L^{1}$ functions with finite total variation. Thanks to the anisotropic diffusion property of the total variation term, the model (1.1) effectively removes Gaussian noise and preserves edges and discontinuities of the image SC:2003 .

Total variation minimization can be applied to not only image denoising, but also other interesting problems in image processing. In CS:2002 , an image inpainting model with the total variation regularizer was proposed. The model in CS:2002 is a simple modification of (1.1), which excludes the inpainting domain from the domain of integration of the fidelity term in (1.1). In addition, (1.1) can be generalized to the image deconvolution problem by replacing $u$ in the fidelity term by $Au$, where $A$ is either a specific convolution operator CW:1998 ; YK:1996 . In order to treat impulse noise, the $TV$-$L^{1}$ model, which uses $L^{1}$ fidelity instead of $L^{2}$ fidelity, was introduced in CE:2005 ; Nikolova:2004 . It is well-known that the $TV$-$L^{1}$ model preserves contrast of the image, while the conventional ROF model does not. We also note that variational image segmentation problem can be represented as the total variation minimization by appropriate change of variables CEN:2006 . One can refer CP:2016 for further results of total variation minimization.

This paper is concerned with domain decomposition methods (DDMs) for such total variation regularized minimization problems. DDMs are suitable for parallel computation since they solve a large scale problem by dividing it into smaller problems and treating them in parallel. While DDMs for elliptic partial differential equations have been successfully developed over past decades, there have been relatively modest achievements in total variation minimization problems due to their own difficulties. At first, the total variation term is nonsmooth and nonseparable. Thus, the energy functional cannot be expressed as the sum of the local energy functionals in the subdomains in general. Even more, the solution space $BV(\Omega)$ allows discontinuities of a solution on the subdomain interfaces, so that it is difficult to impose appropriate boundary conditions to the local problems in the subdomains.

There have been several researches to overcome such difficulties and develop efficient DDMs for the total variation minimization CTWY:2015 ; DCT:2016 ; Fornasier:2007 ; FLS:2010 ; FS:2009 ; HL:2013 ; HL:2015 ; LLWY:2016 ; LN:2017 ; LNP:2018 ; LPP:2019 . The Gauss-Seidel and Jacobi type subspace correction methods for the ROF model were proposed in Fornasier:2007 ; FLS:2010 ; FS:2009 ; LLWY:2016 , and they were extended to the case of mixed $L^{1}/L^{2}$ fidelity in HL:2013 . However, in LN:2017 , a counterexample was provided for convergence of subspace correction methods. In DCT:2016 , a convergent overlapping DDM for the convex Chan–Vese model CEN:2006 was proposed. Recently, a domain decomposition framework for the case of $L^{1}$ fidelity was introduced in LNP:2018 . One of the effective ideas in DDMs for the total variation minimization is to consider the Fenchel–Rockafellar dual formulation of the model:

[TABLE]

instead of the original one. Here, $H_{0}(\operatorname{div};\Omega)$ denotes the space of vector fields $\mathbf{q}$: $\Omega\rightarrow\mathbb{R}^{2}$ such that $\operatorname{div}\mathbf{q}\in L^{2}(\Omega)$ and $\mathbf{q}\cdot\mathbf{n}=0$ on $\partial\Omega$. Consideration of the dual formulation resolves some difficulties mentioned above; the dual energy functional is separable, and the solution space $H_{0}(\operatorname{div};\Omega)$ requires some regularity on the subdomain interfaces. Even if there arises another difficulty of treating the inequality constraint $|\mathbf{p}(x)|\leq 1\hskip 5.69046pt\forall x\in\Omega$, several successful researches have been done CTWY:2015 ; HL:2015 ; LN:2017 ; LPP:2019 .

In this paper, we generalize the primal-dual DDM proposed in LPP:2019 to a wider class of the total variation minimizations. We consider a general model problem

[TABLE]

where $F$: $BV(\Omega)\rightarrow\bar{\mathbb{R}}$ is a proper, convex, and lower semicontinuous functional satisfying additional properties such that $F$ is separable and simple in the sense that (1.2) can be efficiently solved by the first order primal-dual algorithm CP:2011 . Such a class of total variation minimizations contains the ROF model, the $TV$-$L^{1}$ model, their inpainting variants, and the convex Chan–Vese model for image segmentation. As we noted above, we treat the Fenchel–Rockafellar dual problem of (1.2):

[TABLE]

where $F^{*}$ is the Legendre–Fenchel conjugate of $F$. The image domain is decomposed into a number of disjoint subdomains $\left\{\Omega_{s}\right\}_{s=1}^{\mathcal{N}}$, and the continuity along the subdomain interfaces is enforced. By treating the continuity constraints on the subdomain interfaces by the method of Lagrange multipliers, we obtain an equivalent saddle point problem. Application of the Chambolle–Pock primal-dual algorithm CP:2011 to the resulting saddle point problem yields our proposed method, that shows good performance for various total variation minimization problems such as image denoising, inpainting, and segmentation.

The rest of the paper is organized as follows. We present the basic settings for design of DDM in Sect. 2. In Sect. 3, we state the abstract model problem which generalizes various problems in image processing such as image denoising, inpainting, and segmentation. A convergent nonoverlapping DDM for the model problem is proposed in Sect. 4. We apply the proposed method to several image processing problems and compare with the existing state-of-the-art methods in Sect. 5. We conclude the paper with remarks in Sect. 6.

2 Preliminaries

In this section, we present the basic setting for design of DDM. At first, we introduce notations that will be used throughout the paper. Then, the discrete setting for the dual total variation minimization based on the finite element framework is provided.

2.1 Notations

Let $\mathcal{H}$ be the generic $n$-dimensional Hilbert space. For $v=[v_{1},...,v_{n}]^{T},w=[w_{1},...,w_{n}]^{T}\in\mathcal{H}$ and $1\leq p<\infty$, the $p$-norm of $v$ is denoted by

[TABLE]

and the Euclidean inner product of $v$ and $w$ is denoted by

[TABLE]

We may drop the subscript $\mathcal{H}$ if there is no ambiguity.

In this paper, we use the symbol ∗ as superscripts with two different meanings. At first, for a convex functional $F$: $\mathcal{H}\rightarrow\bar{\mathbb{R}}$, $F^{*}$: $\mathcal{H}\rightarrow\bar{\mathbb{R}}$ denotes the Legendre–Fenchel conjugate of $F$, which is defined as

[TABLE]

for every $v^{*}\in\mathcal{H}$. On the other hand, when $A$: $\mathcal{H}\rightarrow\mathcal{H}$ is a linear operator on $\mathcal{H}$, $A^{*}$ denotes the adjoint of $A$, that is,

[TABLE]

for every $v,w\in\mathcal{H}$. It is well-known that the matrix representation of $A^{*}$ is the transpose of the matrix representation of $A$.

For a convex functional $F$: $\mathcal{H}\rightarrow\bar{\mathbb{R}}$, the effective domain of $F$, denoted by $\operatorname{dom}F$, is defined as

[TABLE]

Also, for a convex subset $C$ of $\mathcal{H}$, we denote the relative interior of $C$ by $\operatorname{ri}C$, which is defined as the interior of $C$ when $C$ is regarded as a subset of its affine hull. One may refer Sect. 6 of Rockafellar:2015 for detailed topological properties of $\operatorname{ri}C$.

For a subset $C$ of $\mathcal{H}$, we define the characteristic functional $\chi_{C}$: $\mathcal{H}\rightarrow\bar{\mathbb{R}}$ of $C$ by

[TABLE]

It is clear that the functional $\chi_{C}$ is convex if and only if $C$ is a convex subset of $\mathcal{H}$.

2.2 Fenchel–Rockafellar duality

Throughout this paper, we use the notion of Fenchel–Rockafellar duality frequently. For the sake of completeness, we provide its key features. For more rigorous texts, readers may refer CP:2016 ; Rockafellar:2015 , for instance.

Let $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ be finite-dimensional Hilbert spaces. Consider the minimization problem of the form

[TABLE]

where $K$: $\mathcal{H}_{1}\rightarrow\mathcal{H}_{2}$ is a linear operator and $F$: $\mathcal{H}_{1}\rightarrow\bar{\mathbb{R}}$, $G$: $\mathcal{H}_{2}\rightarrow\bar{\mathbb{R}}$ are proper, convex, lower semicontinuous functionals. We use the assumptions on $F$ and $G$ presented in Corollary 31.2.1 of Rockafellar:2015 .

Assumption 2.1

There exists $u\in\mathcal{H}_{1}$ such that $u\in\operatorname{ri}(\operatorname{dom}F)$ and $Ku\in\operatorname{ri}(\operatorname{dom}G)$.

Under Assumption 2.1, the following relations hold Rockafellar:2015 :

[TABLE]

The saddle point problem in the right hand side of (2.2a) is called the primal-dual formulation of (2.1), and the minimization problem in (2.2b) is called the dual formulation of (2.1). Hence, it is enough to solve (2.2a) or (2.2b) instead of (2.1) in many cases.

2.3 Discrete setting

Let $\Omega\subset\mathbb{R}^{2}$ be a rectangular image domain consisting of a number of rows and columns of pixels. Since each pixel holds a value representing the intensity at a point, we may regard the image as a pixelwise constant function on $\Omega$. We write the collection of all pixels in $\Omega$ as $\mathcal{T}$ and define the space $X$ for an image as

[TABLE]

We note that $X\subset BV(\Omega)$. We regard each pixel in $\mathcal{T}$ as a square element so that $X$ becomes a piecewise constant square finite element space whose side length equals to 1. It is clear that the functions

[TABLE]

form a basis for $X$. For $u\in X$ and $T\in\mathcal{T}$, let $(u)_{T}\in\mathbb{R}$ denote the degree of freedom (dof) of $u$ associated with the basis function $\phi_{T}$. Then, $u$ is represented as

[TABLE]

To obtain a finite element discretization of the dual problem, the space $Y$ for the dual variables is defined by the lowest order Raviart–Thomas elements:

[TABLE]

where $\mathcal{RT}_{0}(T)$ is the collection of the vector functions $\mathbf{q}$: $T\rightarrow\mathbb{R}^{2}$ of the form

[TABLE]

We notice that the divergence operator in the continuous setting is well-defined on $Y$ and $\operatorname{div}Y\subset X$. Each dof of $Y$ is the value of the normal component on a pixel edge. Let $\mathcal{I}$ be the set of indices of the basis functions for $Y$ and $\left\{\bm{\psi}_{i}\right\}_{i\in\mathcal{I}}$ be the basis. For $\mathbf{p}\in Y$ and $i\in\mathcal{I}$, we denote the dof of $\mathbf{p}$ associated with the basis function $\bm{\psi}_{i}$ as $(\mathbf{p})_{i}\in\mathbb{R}$; we have

[TABLE]

Throughout the paper, the $p$-norm of $\mathbf{p}\in Y$ is defined as the $p$-norm of the vector of dofs of $\mathbf{p}$ for $1\leq p<\infty$. In this case, one may regard $\|\cdot\|_{2,Y}$ as a lumped $(L^{2}(\Omega))^{2}$-norm with proper scaling; see Remark 2.2 of LPP:2019 .

In order to treat the inequality constraints appearing in the dual formulation (1.3), we define the convex subset $C$ of $Y$ by

[TABLE]

The projection of $\mathbf{p}\in Y$ onto $C$ can be easily computed by

[TABLE]

2.4 Domain decomposition setting

We decompose the image domain $\Omega$ into $\mathcal{N}=N\times N$ disjoint rectangular subdomains $\left\{\Omega_{s}\right\}_{s=1}^{\mathcal{N}}$. For two subdomains $\Omega_{s}$ and $\Omega_{t}$ ($s<t$) sharing a subdomain edge (interface), let $\Gamma_{st}=\partial\Omega_{s}\cap\partial\Omega_{t}$ be the shared subdomain edge between them. Also, we define the union of the subdomain interfaces $\Gamma$ by $\Gamma=\bigcup_{s<t}\Gamma_{st}$. We assume that $\Gamma$ does not cut through any elements in $\mathcal{T}$.

Now, we define the local function spaces in the subdomains. For $s=1,...,\mathcal{N}$, let $\mathcal{T}_{s}$ be the collection of all pixels in $\Omega_{s}$. The local primal function space $X_{s}$ is defined by

[TABLE]

Obviously, we have $X=\bigoplus_{s=1}^{\mathcal{N}}X_{s}$. Furthermore, we define the natural restriction operator $R_{s}$: $X\rightarrow X_{s}$ by

[TABLE]

Then, the adjoint $R_{s}^{*}$: $X_{s}\rightarrow X$ of $R_{s}$ becomes the extension-by-zero operator

[TABLE]

The local dual function spaces are defined in the tearing-and-interconnectingsense FR:1991 . More precisely, we define the local dual function space $\tilde{Y}_{s}$ as

[TABLE]

Note that the boundary condition $\tilde{\mathbf{q}}_{s}\cdot\mathbf{n}_{s}=0$ is not imposed on $\Gamma\cap\partial\Omega_{s}$ for $\tilde{Y}_{s}$. Thus, $\tilde{Y}_{s}$ has dofs on the pixel edges contained in $\partial\Omega_{s}\cap\Gamma$. Let $\tilde{\mathcal{I}}_{s}$ be the set of indices of the basis functions for $\tilde{Y}_{s}$. Similarly to (2.3), the inequality-constrained subset $\tilde{C}_{s}$ of $\tilde{Y}_{s}$ is defined by

[TABLE]

The orthogonal projection onto $\tilde{C}_{s}$ is computed as

[TABLE]

Let $\tilde{Y}$ and $\tilde{C}$ be the direct sums of the sets $\tilde{Y}_{s}$’s and $\tilde{C}_{s}$’s, respectively. By definition, functions in $\tilde{Y}$ may have discontinuities on $\Gamma$. Let $\mathcal{I}_{\Gamma}$ be the collection of dofs of $Y$ on $\Gamma$. The jump operator $B$: $\tilde{Y}\rightarrow\mathbb{R}^{|\mathcal{I}_{\Gamma}|}$ measures the magnitude of such discontinuities, that is, $B$ is defined as

[TABLE]

where $\tilde{\mathbf{p}}_{s}=\tilde{\mathbf{p}}|_{\Omega_{s}}$. Clearly, there is a natural isomorphism between $Y$ and $\ker B\subset\tilde{Y}$. For later use, we provide an upper bound of the operator norm of $B$ LPP:2019 .

Proposition 2.2

The operator norm of $B$:* $\tilde{Y}\rightarrow\mathbb{R}^{|\mathcal{I}_{\Gamma}|}$ has a bound $\|B\|_{2}^{2}\leq 2$.*

Proposition 2.2 will be used for the estimation of the range of the parameters in the proposed method.

3 The model problem

The model problem we consider in this paper is the total variation regularized convex minimization problem defined on the image domain $\Omega$:

[TABLE]

where $F$: $X\rightarrow\bar{\mathbb{R}}$ is a proper, convex, lower semicontinuous functional, $TV(u)$ is the discrete total variation of $u$ given by

[TABLE]

and $\alpha>0$ is a positive parameter. In addition, we assume that $F$ satisfies the following three assumptions.

Assumption 3.1

$F$ is separable in the sense that there exist proper, convex, lower semicontinuous local energy functionals $F_{s}$: $X_{s}\rightarrow\bar{\mathbb{R}}$ such that

[TABLE]

where $R_{s}$ is the restriction operator defined in (2.5).

Assumption 3.2

For any $u\in X$ and $\sigma>0$, the proximity operator $\operatorname{prox}_{\sigma F}(u)$ defined by

[TABLE]

has a closed-form formula.

Assumption 3.3

$\operatorname{ri}(\operatorname{dom}F)\neq\emptyset$.

Assumption 3.1 makes (3.1) more suitable for designing DDMs. It will be explained in Sect. 4. Also, by Assumption 3.2, we are able to adopt the primal-dual algorithm CP:2011 to solve (3.1). Indeed, if Assumption 3.2 holds, we can solve an equivalent primal-dual form of (3.1),

[TABLE]

by the primal-dual algorithm. For the model problem (3.1), Assumption 3.3 is equivalent to Assumption 2.1 when $G(Ku)=TV(u)$ since $\operatorname{ri}(\operatorname{dom}(\|\cdot\|_{1,Y}))=Y$. That is, Assumption 3.3 is essential for the primal-dual equivalence. It is trivial that all the problems we will deal with in this paper satisfy Assumption 3.3.

Next, we consider the dual formulation of (3.1):

[TABLE]

From (3.3), it is possible to deduce a relationship between solutions of the primal problem (3.1) and the dual problem (3.4):

[TABLE]

As the standard discretization for the total variation minimization is the finite difference method, we provide a relation between the finite difference discretization and our finite element discretization. The following proposition means that a solution of the finite difference discretization of (1.2) can be recovered from a solution of (3.4).

Proposition 3.4

Assume that the image domain $\Omega$ consists of $M\times N$ pixels. Let $\mathbf{p}^{*}\in Y$ be a solution of (3.4). If $u^{*}\in X$ and $\mathbf{p}^{*}$ satisfy the primal-dual relation (3.5), then $u^{*}$ is a solution of the minimization problem

[TABLE]

where $Du$ is the forward finite difference operator

[TABLE]

and $(|Du|)_{ij}=|(Du)_{ij}^{1}|+|(Du)_{ij}^{2}|$.

Proof

It is straightforward by the same argument as Theorem 2.3 of LPP:2019 . ∎

Remark 3.5

In Proposition 3.4, $TV(u)$ equals to the finite difference anisotropic total variation $\||Du|\|_{1}$ due to the structure of the constraint set $C$ in (2.3). On the other hand, one can make $TV(u)$ equal to the finite difference isotropic total variation by replacing $C$ by an appropriate convex subset of $Y$; see LPP:2019 for details.

Remark 3.6

We give a remark on the relation between the continuous problem (1.2) and the discrete problem (3.1). For simplicity, we assume that the resolution of the image is $N\times N$. We fix the side length of $\Omega$ by 1 and introduce a side length parameter $h=1/N$. The space $X^{h}$ is defined in the same manner as in Sect. 2.3. Then the question we have is whether a solution $u^{h*}\in X^{h}$ of the minimization problem

[TABLE]

accumulates at a solution $u^{*}\in BV(\Omega)$ of (1.2) in $BV(\Omega)$ as $h\rightarrow 0$. However, since $TV(u^{h})$ does not $\Gamma$-converge to the $BV$-seminorm in $BV(\Omega)$ as $h\rightarrow 0$ (see Example 4.1 of Bartels:2012 ), this is not true in general. It means that the space $X^{h}$ may not be considered as a correct space for solving (1.2). To avoid this uncomfortable situation, one may use higher order Raviart–Thomas elements to discretize the dual formulation Bartels:2012 ; HHSNW:2018 .

Nevertheless, thanks to Proposition 3.4, one can construct a sequence accumulating at $u^{*}$ in the $L^{1}(\Omega)$-topology from $\{u^{h*}\}_{h>0}$. For any $u^{h}\in X^{h}$, we have

[TABLE]

where $\mathrm{dof}^{h}(u^{h})\in\mathbb{R}^{N\times N}$ is the dofs of $u^{h}$ and $\mathrm{D}^{h}$: $\mathbb{R}^{N\times N}\rightarrow\mathbb{R}^{N\times N}\times\mathbb{R}^{N\times N}$ is the forward finite difference operator defined similarly to (3.6). As it can be shown without much difficulty the right hand side of (3.7) $\Gamma$-converges to the $BV$-seminorm in $L^{1}(\Omega)$ as $h\rightarrow 0$ CLL:2011 ; WL:2011 , there exists an interpolation operator $I^{h}$: $\mathbb{R}^{N\times N}\rightarrow L^{1}(\Omega)$ such that $I^{h}\mathrm{dof}^{h}(u^{h*})$ accumulates at $u^{*}$ in $L^{1}(\Omega)$ as $h\rightarrow 0$. We note that for $u^{h}\in X^{h}$, $I^{h}\mathrm{dof}^{h}(u^{h})\not\in X^{h}$ in general.

The model problem (3.1) occurs in various areas of mathematical image processing. One of the typical examples is the image denoising problem. In ROF:1992 , authors proposed the well-known ROF model which consists of the $L^{2}$-fidelity term and the total variation regularizer. In the ROF model, $F$ is given by

[TABLE]

and Assumptions 3.1 and 3.2 are satisfied with

[TABLE]

In order to preserve the contrast of an image, the $TV$-$L^{1}$ model which uses the $L^{1}$-fidelity term was introduced in CE:2005 ; Nikolova:2004 :

[TABLE]

The local functionals $F_{s}$ and the proximity operator $\operatorname{prox}_{\sigma F}$ are readily obtained as follows:

[TABLE]

where the elementwise shrinkage operator $\operatorname{shr}_{S}$: $X\times\mathbb{R}_{>0}\rightarrow X$ on $S\subset\Omega$ is defined as follows WYYZ:2008 :

[TABLE]

with the convention $\frac{0}{|0|}=0$. For simplicity, $S$ is assumed to be a union of elements of $\mathcal{T}$ in (3.8).

The models for the image denoising problem are easily extended to the image inpainting problem CS:2002 . Let $D\subset\Omega$ be the inpainting domain and $f\in X$ be the known part of an image. We assume that $D$ does not cut through any elements in $\mathcal{T}$. We set $f=0$ on $D$ for simplicity. Also, let $A$: $X\rightarrow X$ be the restriction operator onto $\Omega\setminus D$, that is, $Au=0$ on $D$ for all $u\in X$. Then, $F$ is given by

[TABLE]

Note that Assumptions 3.1 and 3.2 are ensured with

[TABLE]

where $A_{s}u_{s}=R_{s}AR_{s}^{*}u_{s}$ for $s=1,...,\mathcal{N}$. One can easily check that $A=\bigoplus_{s=1}^{\mathcal{N}}A_{s}$. Similarly to the denoising model, we may use the $L^{1}$ fidelity term instead of the $L^{2}$ fidelity term as follows:

[TABLE]

In this case, $F_{s}$ and $\operatorname{prox}_{\sigma F}$ are given by

[TABLE]

Another typical example is the image segmentation problem. In CEN:2006 , authors proposed a convex image segmentation model with the total variation regularizer as follows:

[TABLE]

where $f$ is a given image, $c_{1}$ and $c_{2}$ are predetermined intensity values. Writing $g=(f-c_{1})^{2}-(f-c_{2})^{2}$ in (3.9) yields the following simpler form:

[TABLE]

Then, (3.10) is of the form (3.1) with

[TABLE]

and satisfies Assumptions 3.1 and 3.2 with

[TABLE]

Here, $\operatorname{proj}_{\left\{0\leq\cdot\leq 1\right\}}$ can be computed pointwise like (2.4).

We close this section by mentioning that the image deconvolution problem is not a case of (3.1) in general. For an operator $A$: $X\rightarrow X$ defined by the matrix convolution with some matrix kernel, $A\neq\bigoplus_{s=1}^{\mathcal{N}}R_{s}AR_{s}^{*}$ since the computation of $R_{s}Au$ needs values of $u$ in $\Omega_{s}$ as well as adjacent subdomains. Consequently, Assumption 3.1 cannot be satisfied for the deconvolution problem.

4 Proposed method

In this section, we extend the primal-dual DDM for the ROF model introduced in LPP:2019 to the more general model problem (3.1). The continuity of a solution on the subdomain interfaces is imposed in the dual sense, that is, it is imposed by the method of Lagrange multipliers. As a result, we obtain an equivalent saddle point problem, which is solved by the first order primal-dual algorithm CP:2011 .

We start the section by stating the following simple proposition, which means that the Legendre–Fenchel conjugate is separable if the original functional is separable.

Proposition 4.1

Let $F$: $X\rightarrow\bar{\mathbb{R}}$ be a proper, convex, lower semicontinuous functional satisfying Assumption 3.1. Then, its Legendre–Fenchel conjugate $F^{*}$ also satisfies Assumption 3.1.

Proof

For $u^{*}=\bigoplus_{s=1}^{\mathcal{N}}u_{s}^{*}\in X$, we have

[TABLE]

Letting $(F^{*})_{s}=F_{s}^{*}$ for $s=1,...,\mathcal{N}$ completes the proof. ∎

Thanks to Proposition 4.1, we can transform the dual model problem (3.4) to an equivalent constrained minimization problem

[TABLE]

In order to treat the continuity constraint $B\tilde{\mathbf{p}}=0$, the method of Lagrange multipliers for (4.1) yields the saddle point formulation

[TABLE]

where

[TABLE]

The following proposition summarizes the equivalence between the dual model problem (3.4) and the resulting saddle point problem (4.2).

Proposition 4.2

If $\mathbf{p}^{*}\in C$ is a solution of (3.4), then $\tilde{\mathbf{p}}^{*}=\bigoplus_{s=1}^{\mathcal{N}}\mathbf{p}^{*}|_{\Omega_{s}}\in\tilde{C}$ is a primal solution of (4.2). Conversely, if $\tilde{\mathbf{p}}^{*}$ is a primal solution of (4.2), then $\tilde{\mathbf{p}}^{*}\in\ker B$. Hence $\tilde{\mathbf{p}}^{*}\in C$ and $\tilde{\mathbf{p}}^{*}$ is a solution of (3.4).

Now, we are ready to propose the main algorithm of this paper. In the case of ROF model, recovering a primal solution $u^{*}$ from the computed dual solution $\mathbf{p}^{*}$ can be easily done by the primal-dual relation $u^{*}=f+\operatorname{div}\mathbf{p}^{*}/\alpha$. However, in the general case, the primal solution may not be obtained as in the ROF case since the primal-dual relation (3.5) does not always give an explicit formula for $u^{*}$. Instead, we consider an algorithm to obtain a primal solution $u^{*}$ and a dual solution $\mathbf{p}^{*}$ simultaneously. We begin with the primal-dual algorithm CP:2011 applied to (4.2). In each iteration of the primal-dual algorithm, we need to solve the local problems of the following form:

[TABLE]

for $\hat{\mathbf{p}}_{s}^{(n+1)}\in\tilde{Y}_{s}$. The term $\frac{1}{2\tau}\|\tilde{\mathbf{p}}_{s}-\hat{\mathbf{p}}_{s}^{(n+1)}\|_{2,\tilde{Y}_{s}}^{2}$ appears because we compute proximal descent/ascent in each step of the primal-dual algorithm. To obtain a primal solution $u^{(n+1)}$ and a dual solution $\tilde{\mathbf{p}}^{(n+1)}$ simultaneously, we replace (4.3) by the following primal-dual formulation of (4.3):

[TABLE]

There is another advantage to solve (4.4) instead of (4.3). Differently from the ROF model, it is sometimes cumbersome to get an explicit formula for $F_{s}^{*}$, which makes the design of a local solver difficult. We note that an explicit formula for $F^{*}$ in the case of $F(u)=\|Au-f\|_{1,X}$ with nonsingular $AA^{*}$ is given in DHN:2009 , but it is somewhat complicated. However, considering (4.4) does not require an explicit formula for $F^{*}$.

Similarly to the ROF case, the solution pair $(u^{(n+1)},\tilde{\mathbf{p}}^{(n+1)})$ can be constructed by assembling the local solution pairs $(u_{s}^{(n+1)},\tilde{\mathbf{p}}_{s}^{(n+1)})$ in the subdomains, i.e., $u^{(n+1)}=\bigoplus_{s=1}^{\mathcal{N}}u_{s}^{(n+1)}$ and $\tilde{\mathbf{p}}^{(n+1)}=\bigoplus_{s=1}^{\mathcal{N}}\tilde{\mathbf{p}}_{s}^{(n+1)}$. The local solution pair $(u_{s}^{(n+1)},\tilde{\mathbf{p}}_{s}^{(n+1)})$ in the subdomain $\Omega_{s}$ is obtained by solving the local problem

[TABLE]

and each local problem can be solved in parallel. We will address how to solve (4.5) in Sect. 5 in detail.

In summary, the proposed algorithm is presented in Algorithm 1. As in the ROF case, the range $\tau\sigma<1/2$ comes from Proposition 2.2.

Next, we analyze convergence of the proposed method. Since the sequence $\{(\tilde{\mathbf{p}}^{(n)},\lambda^{(n)})\}$ generated by Algorithm 1 agrees with the one generated by the standard primal-dual algorithm, $O(1/n)$ ergodic convergence is guaranteed.

Theorem 4.3

Let $\{(\tilde{\mathbf{p}}^{(n)},\lambda^{(n)})\}$ be the sequence generated by Algorithm 1. Then, $(\tilde{\mathbf{p}}^{(n)},\lambda^{(n)})$ converges to a saddle point $(\tilde{\mathbf{p}}^{*},\lambda^{*})$ of (4.2) and it satisfies

[TABLE]

for any $\tilde{\mathbf{p}}\in\tilde{Y}$ and $\lambda\in\mathbb{R}^{|\mathcal{I}_{\Gamma}|}$, where $\tilde{\mathbf{p}}_{n}=\frac{1}{n}\sum_{i=1}^{n}\tilde{\mathbf{p}}^{(i)}$ and $\lambda_{n}=\frac{1}{n}\sum_{i=1}^{n}\lambda^{(i)}$.

Proof

Since the term

[TABLE]

in (4.2) is convex, by Theorem 5.1 in CP:2016 , we get the desired result. ∎

By Theorem 4.3, we ensure the convergence of the sequences $\{\tilde{\mathbf{p}}^{(n)}\}$ and $\{\lambda^{(n)}\}$. On the other hand, since the sequence $\{u^{(n)}\}$ is a kind of byproduct of the primal-dual algorithm, the convergence theory for the primal-dual algorithm developed in CP:2011 ; CP:2016 does not ensure global convergence of $\{u^{(n)}\}$. Thus, we will prove that it tends to a solution of (3.1). For the sake of convenience, we rewrite (4.4) as the following more compact form:

[TABLE]

where $\operatorname{\widetilde{\operatorname{div}}}$: $\tilde{Y}\rightarrow X$ is defined as $\operatorname{\widetilde{\operatorname{div}}}\tilde{\mathbf{p}}=\bigoplus_{s=1}^{\mathcal{N}}\operatorname{div}\tilde{\mathbf{p}}_{s}$. Note that if $\tilde{\mathbf{p}}\in\ker B$, then $\tilde{\mathbf{p}}\in Y$ so that $\operatorname{\widetilde{\operatorname{div}}}\tilde{\mathbf{p}}=\operatorname{div}\tilde{\mathbf{p}}$. Furthermore, since $\tilde{Y}$ has dofs on $\Gamma$, we have $\ker\operatorname{\widetilde{\operatorname{div}}}^{*}=\left\{0\right\}$. Then, we readily see that $u^{(n+1)}$ is a solution of the minimization problem

[TABLE]

where $G_{n}$ is defined as

[TABLE]

the Legendre–Fenchel conjugate of

[TABLE]

At first, we verify the boundedness of $\left\{G_{n}\right\}$.

Lemma 4.4

There exist a finite functional $\overline{G}$:* $\tilde{Y}\rightarrow\mathbb{R}$ and a coercive functional $\underline{G}$: $\tilde{Y}\rightarrow\mathbb{R}$ such that*

[TABLE]

Proof

We first see that for any $\tilde{\mathbf{q}}\in\tilde{Y}$, we have $\frac{\tilde{\mathbf{q}}}{|\tilde{\mathbf{q}}|}\in\tilde{C}$, where $|\tilde{\mathbf{q}}|$ is the dof-wise absolute value of $\tilde{\mathbf{q}}$ and the division is done dof-wise with convention $\frac{0}{|0|}=0$. An upper bound of $\left\{G_{n}\right\}$ is obtained as follows:

[TABLE]

Take $\overline{G}(\tilde{\mathbf{q}})=\|\tilde{\mathbf{q}}\|_{1,\tilde{Y}}$, which is clearly finite.

Now, we find a lower bound of $\left\{G_{n}\right\}$. Note that $\hat{\mathbf{p}}^{(n)}=\tilde{\mathbf{p}}^{(n-1)}-\tau B^{*}\lambda^{(n)}$ in Algorithm 1. Since $\{\tilde{\mathbf{p}}^{(n)}\}$ and $\{\lambda^{(n)}\}$ are convergent by Theorem 4.3, $\{\hat{\mathbf{p}}^{(n)}\}$ is also convergent. Hence, $\{\hat{\mathbf{p}}^{(n)}\}$ is bounded, i.e. there exists a constant $L>0$ such that

[TABLE]

Then, we have

[TABLE]

Note that $\tilde{\mathcal{I}}_{s}$ denotes the set of indices of the basis functions for $\tilde{Y}_{s}$ and the value of $\sum_{s=1}^{\mathcal{N}}|\tilde{\mathcal{I}}_{s}|$ depends on the image size and the number of subdomains $\mathcal{N}$. Take

[TABLE]

which is coercive due to the term $\|\tilde{\mathbf{q}}\|_{1,\tilde{Y}}$. ∎

Next lemma provides a criterion for equi-coercivity of a collection of functionals.

Lemma 4.5

Let $\left\{F_{i}\right\}$, $i=1,2,\dots$, be a collection of functionals from $X$ into $\bar{\mathbb{R}}$. If there exists a coercive functional $F_{0}$:* $X\rightarrow\bar{\mathbb{R}}$ such that $F_{i}\geq F_{0}$ for all $i$, then $\left\{F_{i}\right\}$ is equi-coercive, that is, for every $t\geq 0$, there exists a compact subset $K_{t}$ of $X$ such that*

[TABLE]

Proof

Let $t\geq 0$. Since $F_{0}$ is coercive, there exists a compact subset $K_{t}$ of $X$ such that

[TABLE]

On the other hand, for every $i$, $F_{i}\geq F_{0}$ implies that

[TABLE]

Therefore, $\left\{F_{i}\right\}$ is equi-coercive. ∎

Now, with the help of the lemmas above, we state the main theorem, which ensures that the sequence $\{u^{(n)}\}$ generated by Algorithm 1 approaches to solutions of (3.1).

Theorem 4.6

Let $\{u^{(n)}\}$ be the sequence generated by Algorithm 1. Then, $\{u^{(n)}\}$ is bounded and every limit point of $\{u^{(n)}\}$ is a solution of (3.1).

Proof

Recall that $u^{(n)}$ is a solution of the minimization problem

[TABLE]

Since $F$ is proper, we may choose $u_{0}\in X$ with $F(u_{0})<\infty$. By Lemma 4.4, we have

[TABLE]

Thanks to the minimization property of $u^{(n)}$, we get

[TABLE]

By Lemmas 4.4 and 4.5, $\left\{G_{n}\right\}$ is equi-coercive, that is, there exists a compact subset $K_{t}$ of $\tilde{Y}$ independent of $n$ such that

[TABLE]

Thus, $G_{n}(-\operatorname{\widetilde{\operatorname{div}}}^{*}u^{(n)})\leq t$ implies that $-\operatorname{\widetilde{\operatorname{div}}}^{*}u^{(n)}\in K_{t}\cap\operatorname{ran}\operatorname{\widetilde{\operatorname{div}}}^{*}$. Since $\ker\operatorname{\widetilde{\operatorname{div}}}^{*}=\left\{0\right\}$, the map $-\operatorname{\widetilde{\operatorname{div}}}^{*}$ is a continuous isomorphism between $X$ and $\operatorname{ran}\operatorname{\widetilde{\operatorname{div}}}^{*}$. Therefore, we can deduce that

[TABLE]

which is a compact subset of $X$ independent of $n$. This implies that $\{u^{(n)}\}$ is bounded.

Now, we may refine $\{u^{(n)}\}$ so that it converges to its limit point $u^{*}$. Since $(u^{(n)},\tilde{\mathbf{p}}^{(n)})$ is a solution of a saddle point problem (4.6), it satisfies

[TABLE]

By Theorem 4.3, $\{\tilde{\mathbf{p}}^{(n)}\}$ converges to $\tilde{\mathbf{p}}^{*}$. As $u^{(n)}\rightarrow u^{*}$ and $\frac{\operatorname{\widetilde{\operatorname{div}}}\tilde{\mathbf{p}}^{(n)}}{\alpha}\rightarrow\frac{\operatorname{\widetilde{\operatorname{div}}}\tilde{\mathbf{p}}^{*}}{\alpha}$, by closedness of the graph of $\partial F$ (See Theorem 24.4 in Rockafellar:2015 ), we get

[TABLE]

By Proposition 4.2, $\tilde{\mathbf{p}}^{*}\in\ker B$, and hence $\tilde{\mathbf{p}}^{*}\in Y$. We obtain the relation

[TABLE]

From the facts that $\tilde{\mathbf{p}}^{*}$ is a solution of (3.4) (See Proposition 4.2) and $(u^{*},\tilde{\mathbf{p}}^{*})$ satisfies the relation (3.5), we conclude that $u^{*}$ is a solution of (3.1). ∎

As a direct consequence of of Theorem 4.6, we get the following result.

Corollary 4.7

Let $\{u^{(n)}\}$ be the sequence generated by Algorithm 1. Then, $\{\mathcal{J}(u^{(n)})\}$ converges to $\mathcal{J}(u^{*})$, where $u^{*}$ is a solution of (3.1).

5 Applications

In this section, we apply our proposed DDM to various total variation based image processing problems mentioned in Sect. 3. The proposed method was implemented in ANSI C with OpenMPI and compiled by Intel Parallel Studio XE. All the computations were done on a cluster composed of seven machines, where each machine has two Intel Xeon SP-6148 CPUs (2.4GHz, 20C), 192GB memory, and the operating system CentOS 7.4 64bit.

To emphasize efficiency of the proposed method as a parallel algorithm, we compare the wall-clock time of the proposed method with the primal-dual algorithm CP:2011 for the full dimension problem (3.4). The wall-clock time measures the total elapsed time including the communication time.

At each iteration of Algorithm 1, we solve the local saddle point problems of the form

[TABLE]

where $F_{s}$ is given in Assumption 3.1 and

[TABLE]

The primal-dual algorithm for (5.1) consists of computation of the proximity operators of $\alpha F_{s}(u_{s})$ and $G_{s}^{*}(\tilde{\mathbf{p}}_{s})$, that is,

[TABLE]

for some $\sigma_{0},\tau_{0}>0$. For more details, we refer readers to see CP:2011 . The proximity operator of $G_{s}^{*}(\tilde{\mathbf{p}}_{s})$ is computed easily as follows:

[TABLE]

Computation of the proximity operator of $\alpha F_{s}(u_{s})$ depends on the problem to solve. Thus, we will give details in each subsection.

Since $G_{s}^{*}$ is uniformly convex with parameter $1/\tau$, we are able to adopt the $O(1/n^{2})$ convergent primal-dual algorithm (CP:2011, , Algorithm 2) for the local problems. Such acceleration of the local solvers for DDMs was discussed in LNP:2018 ; LPP:2019 .

Next, we provide the setting of the used parameters. We set the parameters for the outer iterations by $\tau=50$ and $\sigma\tau=1/2$ for the image denoising and inpainting problems, and $\tau=1,\sigma\tau=1/2$ for the image segmentation problem. For the $O(1/n^{2})$ convergent algorithms for the local problems, we set $\gamma=1/8\tau$, $\tau_{0}=10$, and $\sigma_{0}\tau_{0}=1/8$. (The same notations for the parameters are used as in CP:2011 ). The full dimension problems ($\mathcal{N}=1$) are solved by the primal-dual algorithm with the same parameters as the corresponding local problems.

Finally, we note that the local solutions $(u_{s}^{(n-1)},\tilde{\mathbf{p}}_{s}^{(n-1)})$ from the previous outer iteration were chosen as initial guesses for the local problems to reduce the number of inner iterations.

5.1 Image denoising

We present the results of numerical experiments for the $TV$-$L^{1}$ model for image denoising:

[TABLE]

We note that numerical results of the proposed method for ROF model were given in LPP:2019 . In (5.2), $\alpha F_{s}(u_{s})$ is given by

[TABLE]

and its proximity operator can be computed as

[TABLE]

where the shrinkage operator $\operatorname{shr}_{s,S}$ on $S\subset\Omega_{s}$ is defined by replacing $\Omega$ and $\mathcal{T}$ by $\Omega_{s}$ and $\mathcal{T}_{s}$ in (3.8), respectively.

We use two test images “Cameraman $2048\times 2048$” and “Boat $2048\times 3072$” with $20\%$ salt-and-pepper noise (See Fig. 1). PSNR denotes peak signal-to-noise ratio. We set $\alpha=1$.

Thanks to Proposition 3.4, it is able to compare the proposed method with the existing methods for (5.2) based on the finite difference discretization. The following algorithms are used for our performance evaluation:

• •

ALG1: Proposed method described in Algorithm 1, $\mathcal{N}=8\times 8$.

• •

LNP: DDM proposed by Lee, Nam, and Park LNP:2018 with the anisotropic total variation, $\mathcal{N}=8\times 8$, $\tau=0.05$, $\sigma\tau=1/8$.

• •

CP: Primal-dual algorithm proposed by Chambolle and Pock CP:2011 with the anisotropic total variation, $\tau=10$, $\sigma\tau=1/8$.

In the numerical experiments, the following stopping criterion is used for outer iterations:

[TABLE]

and the following ones are used for inner iterations of ALG1 and LNP, respectively:

[TABLE]

Fig. 1 shows the denoised images obtained by ALG1 with $\mathcal{N}=16\times 16$ and CP ($\mathcal{N}=1$). In order to highlight the efficiency of the proposed method as a parallel solver, Table 1 shows the performance of the methods with the varying number of subdomains. In addition, to compare the convergence rate of the proposed method with existing algorithms, we present Fig. 2 which shows decay of $\frac{\mathcal{J}(u^{(n)})-\mathcal{J}(u^{*})}{\mathcal{J}(u^{*})}$ for 1000 iterations of three algorithms ALG1, LNP, and CP, where the minimum primal energy $\mathcal{J}(u^{*})$ is computed approximately by $10^{6}$ iterations of the primal-dual algorithm.

5.2 Image inpainting

We first consider the inpainting model with the $L^{2}$ fidelity term:

[TABLE]

Here, $A$: $X\rightarrow X$ is the restriction operator onto $\Omega\setminus D$, so that its matrix representation is a diagonal matrix whose diagonal entries are either [math] or $1$. Thus, computation of the proximity operator of $\alpha F_{s}(u_{s})=\frac{\alpha}{2}\|A_{s}u_{s}-f\|_{2,X_{s}}^{2}$ is as easy as the case of the denoising problem. Indeed, we have

[TABLE]

Two test images “Cameraman $2048\times 2048$” and “Boat $2048\times 3072$” corrupted by additive Gaussian noise with mean [math] and variance $0.05$ and a text mask are used for numerical experiments (See Fig. 3). The model parameter $\alpha$ is set as $\alpha=10$. Table 2 shows the computation results of the proposed method for the varying number of subdomains. Stop conditions (5.3) and (5.4a) are used for outer iterations and inner iterations of ALG1, respectively. Fig. 3 shows the resulting images of ALG1 with $\mathcal{N}=16\times 16$ and CP. Fig. 4 shows decay of $\frac{\mathcal{J}(u^{(n)})-\mathcal{J}(u^{*})}{\mathcal{J}(u^{*})}$ for two algorithms ALG1 and CP. Here, $\mathcal{J}(u^{*})$ is computed by $10^{6}$ iterations of the primal-dual algorithm.

Now, we consider the following $L^{1}$ inpainting model:

[TABLE]

Similarly to the $L^{2}$ inpainting problem, the proximity operator of $\alpha F_{s}(u_{s})=\alpha\|A_{s}u_{s}-f\|_{1,X}$ is given by

[TABLE]

Test images are corrupted by $20\%$ salt-and-pepper noise and the same text mask as the $L^{2}$ inpainting problem  (See Fig. 5). We use $\alpha=1$ as the model parameter. Numerical results of ALG1, LNP, and CP for (5.6) are given in Table 3. Both ALG1 and LNP uses (5.3) as the stopping criterion for outer iterations. Stopping criteria for local problems are given in (5.4a) and (5.4b) for ALG1 and LNP, respectively. The recovered images obtained by ALG1 and CP are given in Fig. 5. Fig. 6 shows decay of the value of $\frac{\mathcal{J}(u^{(n)})-\mathcal{J}(u^{*})}{\mathcal{J}(u^{*})}$ for three algorithms ALG1, LNP, and CP.

5.3 Image segmentation

As we mentioned in Sect. 3, the convex Chan–Vese model for the image segmentation model is represented as

[TABLE]

where $g=(f-c_{1})^{2}-(f-c_{2})^{2}$. The proximity operator of $\alpha F_{s}(u_{s})=\left<u,g\right>_{X}+\chi_{\left\{0\leq u\leq 1\right\}}(u)$ is computed as

[TABLE]

As shown in Fig. 7, we use a test image “Cameraman $2048\times 2048$.” We set the model parameters $\alpha=10$, $c_{1}=0.6$, and $c_{2}=0.1$ heuristically. Also, to convert the results to binary functions, we use a threshold parameter $1/2$. Numerical results of ALG1, DCT, and CP for (5.7) are presented in Table 4, where DCT denotes the following algorithm:

• •

DCT: DDM proposed by Duan, Chang, and Tai DCT:2016 with the anisotropic total variation, $\mathcal{N}=8\times 8$, $\tau=1$, $\sigma\tau=1/8$.

Both ALG1 and DCT uses (5.3) as a stop condition for outer iterations. For stop conditions for local problems, ALG1 uses (5.4a) while DCT uses (5.4b). The segmentation results obtained by ALG1 and CP are provided in Fig. 7. Fig. 8 presents the energy decay of several algorithms containing the proposed method for (5.7). Since the minimum primal energy $\mathcal{J}(u^{*})$ is negative in this case, we plot the value of $\frac{\mathcal{J}(u^{(n)})-\mathcal{J}(u^{*})}{|\mathcal{J}(u^{*})|}$, where $\mathcal{J}(u^{*})$ is computed by $10^{6}$ iterations of the primal-dual algorithm.

5.4 Discussion

As shown in Figs. 1, 3, 5, and 7, the results of the full dimension problem and the proposed method are not visually distinguishable. Note that the resulting images of the proposed method show no trace of the subdomain interfaces. In Tables 1–4, one can observe that solutions for different values of $\mathcal{N}$ have different PSNRs. It means that the algorithms converge to different solutions. This is because (1.2) admits nonunique solutions in general. However, as the stopping criterion (5.3) indicates, their primal energies tend to the minimum.

In Tables 1–4, one can see that the number of outer iterations monotonically increases as the number of subdomains $\mathcal{N}$ increases. We also observe that the numbers of maximum inner iterations of the proposed method are much smaller than the numbers of iterations of the full dimension problem. The reason is that we utilize more accelerated solvers for the local problems than standard ones. Reduction of the numbers of inner iterations makes the proposed method faster.

For every problem, we see that the wall-clock time is decreasing as the number of subdomains $\mathcal{N}$ grows. In particular, large scale images such as “Cameraman $2048\times 2048$” and “Boat $2048\times 3072$” can be processed in a minute with sufficiently many subdomains, while it takes quite a long time with a single domain. This shows efficiency of the proposed method as a parallel solver for image processing.

Finally, as shown in Figs. 2, 4, 6, and 8, the proposed method shows good performance in terms of the decay rate of the primal energy $\mathcal{J}(u^{(n)})$. To be more precise, the primal energy of the proposed method decreases much faster than the one of CP for the full-dimension problem. The proposed method also outperforms LNP which is a recently developed DDM for $TV$-$L^{1}$ problems. Even more, for the segmentation problem, the convergence rate of the proposed method seems to be faster than $O(1/n)$. However, a theoretical evidence for such fast convergence is still missing. It is observed in Figs. 2, 4, 6, and 8 that the primal energy of the proposed method becomes stagnant when the relative error is sufficiently small. This is due to that local problems are solved inexactly by iterative methods in each iteration. Such stagnation of the primal energy is not problematic in practice because the quality of the recovered image becomes acceptable enough before the stagnation starts. Meanwhile, in terms of wall-clock time, the proposed method outperforms LNP for both denoising and inpainting problems while it shows the similar performance to DCT for the segmentation problem.

6 Conclusion

In this paper, we generalized the primal-dual DDM for the ROF model proposed in LPP:2019 to more general total variation minimization problem. The Fenchel–Rockafellar dual of the model problem was considered. We constructed the constrained minimization problem which has domain decomposition structure and is equivalent to the dual model problem. The constrained minimization problem was converted to the equivalent saddle point problem by the method of Lagrange multipliers. The resulting saddle point problem was solved by the first order primal-dual algorithm and convergence of the dual solution was guaranteed. We also proved convergence of the primal solution. Numerical results showed that the proposed method is superior to the existing methods in the sense of convergence rate, and is much faster with sufficiently many subdomains than the full dimension problem.

Even though the proposed DDM is applicable for various total variation regularized problems, we point out that the proposed method is not appropriate for the image deconvolution problem. We will develop a DDM with similar strategy which is applicable for the image deconvolution problem later.