Fast Nonoverlapping Block Jacobi Method for the Dual Rudin--Osher--Fatemi Model

TL;DR

This paper introduces a fast nonoverlapping block Jacobi method for the dual ROF model in image processing, achieving accelerated convergence rates and outperforming existing domain decomposition methods.

Contribution

It proposes a novel accelerated nonoverlapping block Jacobi method with proven $O(1/n^2)$ convergence for the dual ROF model, improving speed over existing methods.

Findings

Convergence rate of $O(1/n)$ for the basic method.

Accelerated method achieves $O(1/n^2)$ convergence.

Faster practical performance compared to existing methods.

Abstract

We consider nonoverlapping domain decomposition methods for the Rudin--Osher--Fatemi~(ROF) model, which is one of the standard models in mathematical image processing. The image domain is partitioned into rectangular subdomains and local problems in subdomains are solved in parallel. Local problems can adopt existing state-of-the-art solvers for the ROF model. We show that the nonoverlapping relaxed block Jacobi method for a dual formulation of the ROF model has the $O(1/n)$ convergence rate of the energy functional, where $n$ is the number of iterations. Moreover, by exploiting the forward-backward splitting structure of the method, we propose an accelerated version whose convergence rate is $O(1/n^2)$. The proposed method converges faster than existing domain decomposition methods both theoretically and practically, while the main computational cost of each iteration remains the same. We also provide the dependence of the convergence rates of the block Jacobi methods on the image size and the number of subdomains. Numerical results for comparisons with existing methods are presented.

Tables

Table 1. Table 1: List of notations used in Sections   2 and  3

notation	description	section
$\Omega$	image domain of size $M\times N$	2.1
$V$	a collection of functions: $\Omega \to ℝ$	2.1
$W$	a collection of functions: $\Omega \to {ℝ}^2$	2.1
$C$	$\{𝐩\in W:\left|{𝐩}_{ij}\right|\le 1\forall (i,j)\in W\}$	2.1
$\nabla$	discrete gradient operator	2.1
$\mathrm{div}$	discrete divergence operator, $\mathrm{div}=-{\nabla }^{\ast }$	2.1
$F$	$F\left(𝐩\right)=\frac{1}{2}{\Vert \mathrm{div}𝐩+\alpha f\Vert }_2^2$, dual energy functional	2.1
${\chi }_C$	characteristic function of $C\subset W$	2.1
$𝒩$	$M_s\times N_s$, number of subdomains in the decomposition ${\left\{{\Omega }_s\right\}}_{s=1}^{𝒩}$	2.2
$N_c$	number of colors on the decomposition ${\left\{{\Omega }_s\right\}}_{s=1}^{𝒩}$	2.2
$S_k$	union of all subdomains with color $k$	2.2
$W_k$	a collection of functions: $S_k\to {ℝ}^2$	2.2
$R_k$	$R_k𝐩={𝐩|}_{S_k}$, restriction operator: $W\to W_k$	2.2
$C_k$	$R_{k}C$	2.2
$F_k$	$F_k({𝐩}_k;𝐫)=F\left(R_k^{\ast }{𝐩}_k+\left(I-R_k^{\ast }{R}_k\right)𝐫\right)$, local energy functional	2.2
$D_k$	Bregman distance associated with $F_k(\cdot ;𝐫)$	2.2
$D$	$D(𝐩,𝐪)={\sum }_{k=1}^{N_c}{D}_k(R_k𝐩,R_k𝐪)$	2.2
${𝒮}_k$	${𝒮}_k\left(𝐩\right)={\arg \min }_{{𝐩}_k\in W_k}\left\{F_k({𝐩}_k;𝐩)+{\chi }_{C_k}\left({𝐩}_k\right)\right\}$	2.2
${\overline{𝒮}}_k$	${\overline{𝒮}}_k\left(𝐩\right)=R_k^{\ast }{𝒮}_k\left(𝐩\right)+\left(I-R_k^{\ast }{R}_k\right)𝐩$	2.2
${𝒫}_k$	${𝒫}_k\left(𝐩\right)={\arg \min }_{{𝐩}_k\in W_k}\left\{F_k(N_c{𝐩}_k-\left(N_c-1\right)R_k𝐩;𝐩)+{\chi }_{C_k}\left({𝐩}_k\right)\right\}$	3.2

Table 2. Table 2: Performance of the fast pre-relaxed block Jacobi method with various values of α 𝛼 \alpha ( 𝒩 = 8 × 8 𝒩 8 8 \mathcal{N}=8\times 8 )

Test image	$\alpha$	PSNR	iter	wall-clock time (sec)
Peppers $512\times 512$	5	23.88	26	0.34
	10	24.55	11	0.15
	20	23.91	8	0.12
Boat $2048\times 3072$	5	24.27	94	46.23
	10	24.86	10	7.08
	20	24.12	6	4.64

Table 3. Table 3: Performance of the fast pre-relaxed block Jacobi method with various numbers of subdomains  𝒩 𝒩 \mathcal{N}  ( α = 10 𝛼 10 \alpha=10 )

Test image	$𝒩$	PSNR	iter	wall-clock time (sec)
1 Peppers $512\times 512$	1	24.55	-	0.71
	$2\times 2$	24.55	9	0.70
	$4\times 4$	24.55	10	0.23
	$8\times 8$	24.55	11	0.15
	$16\times 16$	24.55	14	0.20
1 Boat $2048\times 3072$	1	24.86	-	53.24
	$2\times 2$	24.86	10	50.14
	$4\times 4$	24.86	10	14.47
	$8\times 8$	24.86	10	7.08
	$16\times 16$	24.86	10	4.70

Table 4. Table 4: Performance of the fast pre-relaxed block Jacobi method with various numbers of subdomains  𝒩 𝒩 \mathcal{N} , the stripe-shaped case ( α = 10 𝛼 10 \alpha=10 )

Test image	$𝒩$	PSNR	iter	wall-clock time (sec)
1 Peppers $512\times 512$	1	24.55	-	0.71
	4	24.55	7	0.57
	16	24.55	8	0.22
	64	24.55	13	0.19
	256	24.55	23	1.05
1 Boat $2048\times 3072$	1	24.86	-	53.24
	4	24.86	7	37.54
	16	24.86	7	10.53
	64	24.86	9	6.86
	256	24.86	14	7.01

Equations

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

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

$$\Omega =\left\{(i,j):1\le i\le M,1\le j\le N\right\}.$$

$$C=\left\{\mathbf{p}\in W:|{\mathbf{p}}_{ij}|\le 1\forall (i,j)\in \Omega \right\},$$

$$(\nabla u{)}_{ij}^1$$

$$(\nabla u{)}_{ij}^2$$

$$(\mathrm{div}\mathbf{p}{)}_{ij}$$

$$\underset{u\in V}{\min }\left\{\frac{\alpha }{2}\Vert u-f{\Vert }_2^2+\Vert \nabla u{\Vert }_1\right\},$$

$$\Vert \nabla u{\Vert }_1=\sum _{(i,j)\in \Omega }|(\nabla u{)}_{ij}|.$$

$$F(\mathbf{p})=\frac{1}{2}\Vert \mathrm{div}\mathbf{p}+\alpha f{\Vert }_2^2,\mathbf{p}\in W.$$

$$\underset{\mathbf{p}\in W}{\min }\left\{F(\mathbf{p})+{\chi }_C(\mathbf{p})\right\}$$

$$u^{\ast }=f+\frac{1}{\alpha }\mathrm{div}{\mathbf{p}}^{\ast }.$$

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

$$\Omega =\bigcup _{k=1}^{N_c}{S}_k,S_k\cap S_j=\varnothing \forall k,j.$$

$$R_k\mathbf{p}=\mathbf{p}{|}_{S_k},\mathbf{p}\in W.$$

$$(R_k^{\ast }{\mathbf{p}}_k{)}_{ij}=\{\begin{array}{ll}({\mathbf{p}}_k{)}_{ij}& \text{if}(i,j)\in S_k,\\ 0& \text{if}(i,j)\notin S_k,\end{array}{\mathbf{p}}_k\in W_k.$$

$$\sum _{k=1}^{N_c}{R}_k^{\ast }{R}_k=I.$$

$${\chi }_C(\mathbf{p})=\sum _{k=1}^{N_c}{\chi }_{C_k}(R_k\mathbf{p}),\mathbf{p}\in W.$$

$$F_k({\mathbf{p}}_k;\mathbf{r})=F(R_k^{\ast }{\mathbf{p}}_k+(I-R_k^{\ast }{R}_k)\mathbf{r}),{\mathbf{p}}_k\in W_k.$$

$$D_k({\mathbf{p}}_k,{\mathbf{q}}_k)=F_k({\mathbf{p}}_k;\mathbf{r})-F_k({\mathbf{q}}_k;\mathbf{r})-⟨\nabla F_k({\mathbf{q}}_k;\mathbf{r}),{\mathbf{p}}_k-{\mathbf{q}}_k⟩,{\mathbf{p}}_k,{\mathbf{q}}_k\in W_k,$$

$$D_k({\mathbf{p}}_k,{\mathbf{q}}_k)=\frac{1}{2}\Vert \mathrm{div}{R}_k^{\ast }({\mathbf{p}}_k-{\mathbf{q}}_k){\Vert }_2^2.$$

$$D(\mathbf{p},\mathbf{q})=\sum _{k=1}^{N_c}{D}_k(R_k\mathbf{p},R_k\mathbf{q}),\mathbf{p},\mathbf{q}\in W.$$

$$F(\mathbf{p})\le F(\mathbf{q})+⟨\nabla F(\mathbf{q}),\mathbf{p}-\mathbf{q}⟩+N_{c}D(\mathbf{p},\mathbf{q}).$$

$$F(\mathbf{p})-F(\mathbf{q})$$

$$=\frac{1}{2}{\Vert \mathrm{div}\left(\sum _{k=1}^{N_c}{R}_k^{\ast }{R}_k(\mathbf{p}-\mathbf{q})\right)\Vert }_2^2$$

$$\le N_c\sum _{k=1}^{N_c}\frac{1}{2}\Vert \mathrm{div}{R}_k^{\ast }{R}_k(\mathbf{p}-\mathbf{q}){\Vert }_2^2=N_{c}D(\mathbf{p},\mathbf{q}).$$

$${\mathcal{S}}_k(\mathbf{p})=\underset{{\mathbf{p}}_k\in W_k}{\mathrm{argmin}}\left\{F_k({\mathbf{p}}_k;\mathbf{p})+{\chi }_{C_k}({\mathbf{p}}_k)\right\},\mathbf{p}\in W,$$

$${\bar{\mathcal{S}}}_k(\mathbf{p})=R_k^{\ast }{\mathcal{S}}_k(\mathbf{p})+(I-R_k^{\ast }{R}_k)\mathbf{p},\mathbf{p}\in W.$$

$$F_k({\mathbf{p}}_k;\mathbf{q})-F({\bar{\mathcal{S}}}_k(\mathbf{q}))\ge \frac{1}{2}\Vert \mathrm{div}{R}_k^{\ast }({\mathbf{p}}_k-R_k{\bar{\mathcal{S}}}_k(\mathbf{q})){\Vert }_2^2.$$

$${\mathbf{p}}^{(n+1)}=\frac{1}{N_c}\sum _{k=1}^{N_c}{\bar{\mathcal{S}}}_k({\mathbf{p}}^{(n)})=\left(1-\frac{1}{N_c}\right){\mathbf{p}}^{(n)}+\frac{1}{N_c}\sum _{k=1}^{N_c}{R}_k^{\ast }{\mathcal{S}}_k({\mathbf{p}}^{(n)})$$

$$\frac{1}{N_c}F(\mathbf{p})+\left(1-\frac{1}{N_c}\right)F(\mathbf{q})-F\left(\frac{1}{N_c}\sum _{k=1}^{N_c}{\bar{\mathcal{S}}}_k(\mathbf{q})\right)\ge D\left(\mathbf{p},\frac{1}{N_c}\sum _{k=1}^{N_c}{\bar{\mathcal{S}}}_k(\mathbf{q})\right)-D(\mathbf{p},\mathbf{q}).$$

$$\mathcal{S}(\mathbf{q})=\frac{1}{N_c}\sum _{k=1}^{N_c}{\bar{\mathcal{S}}}_k(\mathbf{q}).$$

Full text

{AMS}

65N55, 65Y05, 65B99, 65K10, 68U10

1 Introduction

As large-scale images have become available in these days, it takes a long time to process such images and demands to reduce the time, especially in real time, are required. To accomplish such goal, there has been arisen the necessity of mathematical methods to make efficient use of distributed memory computers. In this perspective, we consider domain decomposition methods (DDMs), which have been already used successfully for elliptic partial differential equations [20, 24].

In DDMs, a large-scale problem defined on a domain $\Omega$ is splitted into smaller ones called local problems defined on subdomains $\left\{\Omega_{s}\right\}_{s=1}^{\mathcal{N}}$ satisfying $\bar{\Omega}=\bigcup_{s=1}^{\mathcal{N}}\bar{\Omega}_{s}$. We assign one local problem per processor so that all local problems can be solved in parallel in distributed memory computers. A solution of the full-dimension problem is obtained by the assembly of solutions of local problems. Thus, DDMs are suitable for distributed memory computers. Block Jacobi and Gauss–Seidel methods are typical examples of DDMs.

In this paper, we consider block Jacobi methods for the Rudin–Osher–Fatemi (ROF) model [22], which is one of the standard models in variational image processing given by

[TABLE]

where $\alpha$ is a positive parameter, $f$ is a corrupted image, $TV(u)$ is the total variation of $u$, and $BV(\Omega)$ is the space of functions of bounded variation on an image domain $\Omega\subset\mathbb{R}^{2}$. There have been numerous approaches to solve Eq. 1; see [2, 5, 6, 12, 26, 27] for instance.

There have been several remarkable preceding researches on DDMs for variational image processing. Subspace correction methods for Eq. 1 were proposed in [10, 11], and various techniques to deal with local problems were considered in [15, 16]. In [13], subspace correction methods for the case of mixed $L^{1}$/$L^{2}$-fidelity were considered. However, it was shown in [17] that such subspace correction methods for Eq. 1 may not converge to a minimizer. To ensure convergence to a minimizer, in [8, 14, 17], subspace correction methods for the Fenchel–Rockafellar dual formulation of Eq. 1 given by

[TABLE]

were proposed, where $C_{0}^{1}(\Omega;\mathbb{R}^{2})$ is the set of continuously differentiable vector fields $\mathbf{p}$: $\Omega\rightarrow\mathbb{R}^{2}$ with compact support. Recently, interface-based domain decomposition methods for total variation minimization were proposed [18, 19].

Over the past decade, block coordinate descent type algorithms for nonsmooth convex optimization have been successfully developed [3, 25]. They solve a minimization problem by successively minimizing the functional with respect to one of the coordinate blocks. In the context of domain decomposition methods, they correspond to nonoverlapping block Gauss–Seidel methods. Indeed, we observe that convergence rate analyses for block coordinate descent methods given in [7, 23] can be directly applied to the nonoverlapping block Gauss–Seidel method for Eq. 2 in [14]. However, to the best of our knowledge, there are no available analysis of the convergence rate for the nonoverlapping block Jacobi method for Eq. 2.

In this paper, we prove the $O(1/n)$ convergence rate of the nonoverlapping relaxed block Jacobi method for Eq. 2, which was proposed in [14]. We note that it was shown in [8] that the overlapping counterpart is $O(1/n)$ convergent and the authors gave a remark that their proof cannot be directly applied to the nonoverlapping case. In addition to the analysis of the convergence rate, we propose an accelerated nonoverlapping block Jacobi method for Eq. 2 which has the $O(1/n^{2})$ convergence rate. First, we modify the relaxed block Jacobi method so that it has the forward-backward splitting structure [9]. We call such a modified method as pre-relaxed block Jacobi method. Thanks to its forward-backward structure, we obtain an $O(1/n^{2})$ convergent algorithm by adding the momentum technique introduced in [2]. Numerical results ensure faster convergence of the proposed method compared to the existing ones.

The rest of the paper is organized as follows. Basic settings for the paper are presented in Section 2. We analyze the convergence rate of the relaxed block Jacobi method and propose its variants including an $O(1/n^{2})$ convergent method in Section 3. We consider how to deal with local problems of the proposed methods in Section 4. We compare the proposed methods to the existing ones by numerical experiments in Section 5. We conclude the paper with remarks in Section 6.

2 Preliminaries

In this section, we review the standard finite difference framework for the ROF model. Then, we formulate the Fenchel–Rockafellar dual problem of the discrete ROF model. Finally, a nonoverlapping domain decomposition setting for the dual problem is presented. A list of notations used in Sections 2 and 3 is given in Table 1.

2.1 Discrete setting

We assume that the grayscale image domain $\Omega$ is composed of $M\times N$ pixels; i.e.,

[TABLE]

Let $V$ and $W$ be the collections of functions from $\Omega$ into $\mathbb{R}$ and $\mathbb{R}^{2}$, respectively. Both spaces are equipped with the usual Euclidean inner products $\left<\cdot,\cdot\right>$ and their induced norms $\|\cdot\|_{2}$. In addition, let $C$ be the convex subset of $W$ given by

[TABLE]

where $|\mathbf{p}_{ij}|=(|p_{ij}^{1}|^{2}+|p_{ij}^{2}|^{2})^{1/2}$ for $\mathbf{p}=(p^{1},p^{2})$.

The discrete gradient operator $\nabla$: $V\rightarrow W$ is defined by

[TABLE]

We define the discrete divergence operator $\mathrm{div}$: $W\rightarrow V$ as the minus adjoint of $\nabla$, i.e.,

[TABLE]

With the operators defined above, a discrete version of Eq. 1 is stated as

[TABLE]

where $\|\nabla u\|_{1}$ is the discrete total variation of $u$ defined as the 1-norm of $\nabla u$:

[TABLE]

Now, we consider the Fenchel–Rockafellar dual problem of Eq. 5. We define the energy functional $F$: $W\rightarrow\mathbb{R}$ by

[TABLE]

It is well-known that a solution $u^{*}$ of Eq. 5 can be recovered from a solution $\mathbf{p}^{*}$ of

[TABLE]

by a simple algebraic formula

[TABLE]

Here $\chi_{C}$: $W\rightarrow\bar{\mathbb{R}}$ is the characteristic function of $C\subset W$ defined as

[TABLE]

We call Eq. 6 the Fenchel–Rockafellar dual problem of Eq. 5. For more details, readers may refer [5, 21].

2.2 Domain decomposition setting

First, we partition the image domain $\Omega$ into $\mathcal{N}=M_{s}\times N_{s}$ nonoverlapping rectangular subdomains $\left\{\Omega_{s}\right\}_{s=1}^{\mathcal{N}}$. All subdomains can be classified into $N_{c}$ classes (colors) by the usual coloring technique (see section 2.5.1 of [24]); local problems on subdomains in the same class are solved independently. Details of the coloring technique will be given in Section 5. Let $S_{k}$ be the union of all subdomains with color $k$ for $k=1,\dots,N_{c}$. Then, we have

[TABLE]

For each $k=1,\dots,N_{c}$, we define the local function space $W_{k}$ as the collection of functions from $S_{k}$ to $\mathbb{R}^{2}$. Also, we define the restriction operator $R_{k}$: $W\rightarrow W_{k}$ as

[TABLE]

Then, its adjoint $R_{k}^{*}$: $W_{k}\rightarrow W$ becomes the natural extension operator:

[TABLE]

Note that $R_{k}^{*}R_{k}$ is the orthogonal projection from $W$ onto $W_{k}$ and

[TABLE]

We clearly have $W=\bigoplus_{k=1}^{N_{c}}R_{k}^{*}W_{k}$. Similarly, by defining $C_{k}=R_{k}C$, we have $C=\bigoplus_{k=1}^{N_{c}}R_{k}^{*}C_{k}$ and

[TABLE]

Given $\mathbf{r}\in W$, the local energy functional $F_{k}$: $W_{k}\rightarrow\mathbb{R}$ is defined as

[TABLE]

The Bregman distance [4] associated with $F_{k}(\cdot;\mathbf{r})$ is denoted by $D_{k}$, i.e.,

[TABLE]

and one can observe that it is independent of $\mathbf{r}$. Indeed, we have

[TABLE]

We define the function $D$: $W\times W\rightarrow\mathbb{R}$ as

[TABLE]

We provide the descent lemma corresponding to $D$.

Lemma 2.1**.**

For any $\mathbf{p},\mathbf{q}\in W$, we have

[TABLE]

Proof 2.2**.**

By direct calculation, we have

[TABLE]

Finally, we define two operators $\mathcal{S}_{k}$: $W\rightarrow W_{k}$ and $\bar{\mathcal{S}}_{k}$: $W\rightarrow W$ by

[TABLE]

and

[TABLE]

As the local energy functional $F_{k}$ is not strongly convex, the minimization problem in Eq. 11 may admit nonunique minimizers. In this case, we take $\mathcal{S}_{k}(\mathbf{p})$ as any one among them. By the optimality of $\mathcal{S}_{k}$, we have the following lemma.

Lemma 2.3**.**

For any $\mathbf{p}_{k}\in C_{k}$ and $\mathbf{q}\in C$, we have

[TABLE]

Proof 2.4**.**

*It is clear from the optimality condition of $\mathcal{S}_{k}(\mathbf{q})$ (see the proof of Lemma 3.3 in [8]). *

Lemma 2.3 shall be useful in the convergence analysis of block Jacobi methods in the next section. We note that it was also used in the case of overlapping domain decomposition [8].

3 Block Jacobi methods

In this section, we analyze the convergence rate of the relaxed block Jacobi method (Algorithm 2 in [14]) for Eq. 6. Then, by observing a resemblance between the relaxed block Jacobi method and the forward-backward splitting [2, 9], we propose a variant of the method called the pre-relaxed block Jacobi method, which has the exact forward-backward splitting structure. Based on the pre-relaxed block Jacobi method, we also propose an accelerated block Jacobi method which has the $O(1/n^{2})$ convergence rate.

3.1 Relaxed Block Jacobi method

The relaxed block Jacobi method for Eq. 6 is presented in the following, which was first proposed in [14].

Note that

[TABLE]

in Algorithm 1. That is, $\mathbf{p}^{(n+1)}$ is obtained by relaxation of $\sum_{k=1}^{N_{c}}R_{k}^{*}\mathcal{S}_{k}(\mathbf{p}^{(n)})$ with $\mathbf{p}^{(n)}$, and it is the reason why we call Algorithm 1 the relaxed block Jacobi method.

In [14], the convergence of Algorithm 1 was shown without the convergence rate. It is summarized in the following proposition.

Proposition 3.1**.**

*Let $\left\{\mathbf{p}^{(n)}\right\}$ be the sequence generated by Algorithm 1. Then, it satisfies the followings:

(i) The sequence $\left\{F(\mathbf{p}^{(n)})\right\}$ is decreasing, so that it converges.

(ii) The sequence $\left\{\mathbf{p}^{(n)}\right\}$ is bounded and has a limit point which is a solution of Eq. 6. *

First, we prove that the convergence rate of Algorithm 1 is $O(1/n)$. The following lemma is a main ingredient for our proof.

Lemma 3.2**.**

For any $\mathbf{p},\mathbf{q}\in C$, we have

[TABLE]

Proof 3.3**.**

For the sake of convenience, let

[TABLE]

By the convexity of $F$, we have

[TABLE]

Invoking Lemma 2.3 with $\mathbf{p}_{k}=R_{k}\mathbf{p}$ and summing over $k=1,\dots,N_{c}$ give

[TABLE]

Also, by the relation

[TABLE]

and the convexity of the functional $\frac{1}{2}\|\mathrm{div}\cdot\|_{2}^{2}$, we obtain

[TABLE]

Summation of Eq. 14 over $k=1,\dots,N_{c}$ yields

[TABLE]

On the other hand, we have

[TABLE]

*where the last inequality is due to the convexity of $F$. Summation of Eqs. 12, 13, 15 and 16 yields the desired result. *

Now, we are ready to prove $O(1/n)$ convergence of Algorithm 1. The proof presented in this paper uses the same argument as the proof of Theorem 3.1 in [2].

Theorem 3.4**.**

Let $\left\{\mathbf{p}^{(n)}\right\}$ be the sequence generated by Algorithm 1 and $\mathbf{p}^{*}$ be a solution of Eq. 6. Then, for any $n\geq 1$, we have

[TABLE]

Proof 3.5**.**

Using Lemma 3.2 with $\mathbf{p}=\mathbf{p}^{*}$ and $\mathbf{q}=\mathbf{p}^{(j)}$, we get

[TABLE]

Summation of this inequality over $j=0,1,\dots,n-1$ yields

[TABLE]

Also, by Lemma 3.2 with $\mathbf{p}=\mathbf{q}=\mathbf{p}^{(j)}$, we obtain

[TABLE]

Multiplying this inequality by $j$ and summing over $j=1,\dots,n-1$ yields

[TABLE]

Adding Eq. 17 times $N_{c}$ and Eq. 18, we obtain

[TABLE]

That is,

[TABLE]

Therefore, we conclude that

[TABLE]

Remark 3.6**.**

*In [8], the relaxed block Jacobi method with overlapping domain decomposition has $O(1/n)$ convergence rate with the constant that tends to $\infty$ as the overlapping size tends to 0. Therefore, the proof in [8] does not guarantee the convergence of the method in the nonoverlapping case. *

By the definition of $D(\mathbf{p}^{*},\mathbf{p}^{(0)})$ in Eq. 10, the estimate given in Theorem 3.4 depends on the number of subdomains $\mathcal{N}$. The following lemma describes such dependency in detail.

Lemma 3.7**.**

Assume that $\Omega$ with $M\times N$ pixels is partitioned into $\mathcal{N}=M_{s}\times N_{s}$ nonoverlapping rectangular subdomains of the same size. For any $\mathbf{p}\in W$, we have

[TABLE]

where $c_{1}$ is given by

[TABLE]

Proof 3.8**.**

For $\mathbf{p}\in W$, we define $R^{s}\mathbf{p}=\mathbf{p}|_{\Omega_{s}}$ for $s=1,\dots,\mathcal{N}$. Then, we clearly have

[TABLE]

For $(i,j)\in\Omega$, by the definition of $\mathrm{div}$, one requires $\mathbf{p}_{ij}$, $\mathbf{p}_{i-1,j}$, and $\mathbf{p}_{i,j-1}$ to compute $(\mathrm{div}\mathbf{p})_{ij}$. Let $\tilde{\Omega}_{s}$ and $\Omega_{s}^{\circ}$ be subsets of $\Omega$ defined by

[TABLE]

and

[TABLE]

respectively; see Fig. 1 for graphical description. Since

[TABLE]

one can observe that

[TABLE]

and

[TABLE]

Thus, we obtain that

[TABLE]

where $E_{1}$, $E_{2}$, $E_{3}$, $E_{4}$, and $E_{5}$ are regions indicated in Fig. 1(c). Note that the cardinalities of these regions are $\frac{N}{N_{s}}-1$, $\frac{N}{N_{s}}$, $\frac{M}{M_{s}}-1$, $\frac{M}{M_{s}}$, and $1$, respectively. By the Cauchy–Schwarz inequality, we have

[TABLE]

Similarly, one obtains the followings:

[TABLE]

Combining all inequalities stated above, we conclude that

[TABLE]

Summing the last inequality for $s=1,\dots,\mathcal{N}$ yields

[TABLE]

*This completes the proof. *

Remark 3.9**.**

In the case of the stripe-shaped domain decomposition, say $M_{s}=1$ and $N_{s}=\mathcal{N}$, one may obtain a sharper estimate than Lemma 3.7 since we do not need to deal with the terms related to $E_{1}$ and $E_{2}$ in LABEL:edges. As a result, we obtain

[TABLE]

As a direct consequence of Lemma 3.7, we get the following upper bound for $D(\mathbf{p}^{*},\mathbf{p})$ for any $\mathbf{p}\in C$.

Corollary 3.10**.**

Let $\mathbf{p}^{*}$ be a solution of Eq. 6. Then, for any $\mathbf{p}\in C$, we have

[TABLE]

*where $c_{1}$ is given in Eq. 19. *

Proof 3.11**.**

From the optimality condition of $\mathbf{p}^{*}$, we get

[TABLE]

Also, since both $\mathbf{p}$ and $\mathbf{p}^{*}$ are in $C$, we have

[TABLE]

With Lemmas 3.7, 21 and 22, we readily obtain

[TABLE]

Combining Theorems 3.4 and 3.10, we obtain the following estimate of convergence rate, in which dependence on the size of the image and the number of subdomains is revealed.

Corollary 3.12**.**

Let $\left\{\mathbf{p}^{(n)}\right\}$ be the sequence generated by Algorithm 1 and $\mathbf{p}^{*}$ be a solution of Eq. 6. Then, for any $n\geq 1$, we have

[TABLE]

*where $c_{1}$ is given in Eq. 19. *

Remark 3.13**.**

One can obtain a relation between the convergence rate and the total length of the subdomain interfaces from Corollary 3.12. Let $\Gamma_{st}=\partial\Omega_{s}\cap\partial\Omega_{t}$ for $s<t$ and $\Gamma=\bigcup_{s<t}\Gamma_{st}$. If the decomposition $\{\Omega_{s}\}_{s=1}^{\mathcal{N}}$ is of $M_{s}\times N_{s}$, it is clear that the length of $\Gamma$ is given by

[TABLE]

For fixed $M$ and $N$, $c_{1}$ satisfies

[TABLE]

*where $c_{2}$ is a constant independent of $M_{s}$ and $N_{s}$. Therefore, it is expected that the convergence rate of Algorithm 1 may deteriorate as $|\Gamma|$ increases. Using Remark 3.9, the same bound as Eq. 23 can be yielded for the stripe-shaped domain decomposition. *

3.2 Pre-relaxed Block Jacobi method

First, we observe that Lemma 3.2 has the similar form to one in the fundamental lemma of the forward-backward splitting (see Lemma 2.3 of [2]):

[TABLE]

where $\mathcal{T}$: $W\rightarrow W$ is a certain operator which computes the next iterate of the algorithm, and $L$ is a constant. If we replace $F(\mathbf{p})$ in Eq. 24 by $\frac{1}{N_{c}}F(\mathbf{p})+\frac{N_{c}-1}{N_{c}}F(\mathbf{q})$, then we obtain the formula in Lemma 3.2. Note that an $O(1/n)$ convergent algorithm with the forward-backward splitting Eq. 24 has the FISTA acceleration [2]. However, even though Algorithm 1 has the $O(1/n)$ convergence rate, it is difficult to apply the FISTA acceleration because of such differences. Since the relaxation step is crucial to guarantee the convergence of the algorithm [8], to apply the FISTA acceleration, we modify the algorithm so that a formula of the form Eq. 24 can be yielded while it still has the relaxation step.

We define an alternative local solution operator $\mathcal{P}_{k}$: $W\rightarrow W_{k}$ by

[TABLE]

The difference between two operators $\mathcal{S}_{k}$ and $\mathcal{P}_{k}$ is that, $\mathbf{p}_{k}$ in $\mathcal{S}_{k}$ is replaced by the relaxed term $N_{c}\mathbf{p}_{k}-(N_{c}-1)R_{k}\mathbf{p}$ in $\mathcal{P}_{k}$. Now, we propose a variant of Algorithm 1 which uses $\mathcal{P}_{k}$ as the local solution operator; see Algorithm 2.

Compared to Algorithm 1, the relaxation step is inside the local solution operator in Algorithm 2. That is, the relaxation step is done before computing the local solution operator. This is why we call Algorithm 2 the pre-relaxed block Jacobi method. The following lemma says that an application of the operator $\sum_{k=1}^{N_{c}}R_{k}^{*}\mathcal{P}_{k}$ is in fact a proximal descent step with respect to a pseudometric

[TABLE]

Lemma 3.14**.**

For $\mathbf{q}\in C$, $\sum_{k=1}^{N_{c}}R_{k}^{*}\mathcal{P}_{k}(\mathbf{q})$ is a solution of the minimization problem

[TABLE]

Proof 3.15**.**

For $\mathbf{p}_{k}\in W_{k}$, let

[TABLE]

Let us define

[TABLE]

Then, we have

[TABLE]

Thus, the optimality condition of Eq. 25 reads as

[TABLE]

Assembly of Eq. 28 for $k=1,\dots,N_{c}$ gives

[TABLE]

*which is the optimality condition of Eq. 26. Therefore, $\sum_{k=1}^{N_{c}}R_{k}^{*}\mathcal{P}_{k}(\mathbf{q})$ is a solution of Eq. 26. *

Together with Lemmas 2.1 and 3.14, we obtain a formula for Algorithm 2 which has the exactly same form as Eq. 24. The fundamental lemma for Algorithm 2 is presented in the following.

Lemma 3.16**.**

For any $\mathbf{p},\mathbf{q}\in C$, we have

[TABLE]

Proof 3.17**.**

For the sake of convenience, we define $A$ as in Eq. 27 and

[TABLE]

Invoking Lemma 2.1 yields

[TABLE]

By Lemma 3.14, $\mathcal{P}(\mathbf{q})$ satisfies Eq. 29, that is,

[TABLE]

Also, by the convexity of $F$, we have

[TABLE]

*Then, we get the desired result by summation of Eqs. 30, 31 and 32. *

Since Lemma 3.16 has exactly the same form as Lemma 2.3 of [2], analysis of the $O(1/n)$ convergence rate is straightforward.

Theorem 3.18**.**

Let $\left\{\mathbf{p}^{(n)}\right\}$ be the sequence generated by Algorithm 2 and $\mathbf{p}^{*}$ be a solution of Eq. 6. Then, for any $n\geq 1$, we have

[TABLE]

*where $c_{1}$ is given in Eq. 19. *

Proof 3.19**.**

In Theorem 3.1 of [2], we replace $L(f)$ and $\|\cdot\|$ by $N_{c}$ and

[TABLE]

respectively, the notations we use. Consequently, we have the following estimate:

[TABLE]

*Application of Corollary 3.10 to Eq. 33 yields the desired result. *

Remark 3.20**.**

*In Theorem 3.18, we gave a little sharper bound than Theorem 3.4. Indeed, our numerical experiments presented in Section 5 will show that Theorem 3.18 has a faster energy decay than Theorem 3.4 at initial dozens of iterations. *

In addition, an accelerated version of Algorithm 2 with the $O(1/n^{2})$ convergence rate can be designed by the same way as FISTA [2]. It is summarized in Algorithm 3.

Finally, we state the convergence theorem for Algorithm 3. Theorem 3.21 says that Algorithm 3 is superior to Algorithms 1 and 2. In Section 5, we will present some numerical examples where Algorithm 3 shows better performance than others.

Theorem 3.21**.**

Let $\left\{\mathbf{p}^{(n)}\right\}$ be the sequence generated by Algorithm 3 and $\mathbf{p}^{*}$ be a solution of Eq. 6. Then, for any $n\geq 1$, we have

[TABLE]

*where $c_{1}$ is given in Eq. 19. *

Proof 3.22**.**

With the same argument as in Theorem 4.4 of [2], we obtain

[TABLE]

*Then, application of Corollary 3.10 to Eq. 34 completes the proof. *

Remark 3.23**.**

*In Theorem 3.21, the $O(1/n^{2})$ convergence is satisfied regardless of $f$ and $\alpha$ since they only affect on $F$. That is, neither the noise level of a given image nor the value of the weight parameter $\alpha$ interferes with the $O(1/n^{2})$ convergence of Algorithm 3. Similar arguments apply to Algorithms 1 and 2. *

4 Local problems

In this section, we discuss how to deal with local problems of the block Jacobi methods introduced in this paper. The coloring technique used in the proposed methods is also considered.

First, we assume that the domain decomposition $\left\{\Omega_{s}\right\}_{s=1}^{\mathcal{N}}$ is of the window shape as shown in Fig. 2(a). All subdomains $\Omega_{s}$ have the same size. Let $\Omega_{ij}$ be the subdomain on the $i$-th row and the $j$-th column. We define the local function space $W^{s}$ as the collection of functions from $\Omega_{s}$ to $\mathbb{R}^{2}$. Also, we define the restriction operator $R^{s}$: $W\rightarrow W^{s}$ similarly to Eq. 7 and $C^{s}=R^{s}C$.

Local problems in $\Omega_{s}$ have the following general form:

[TABLE]

for some $g\in V$. Let $\tilde{\Omega}_{s}$ be the subset of $\Omega$ containing $\Omega_{s}$ and its adjacent lines of pixels on the bottom and the right sides; see Fig. 1(a). Also, we define $\tilde{V}^{s}$ as the collection of functions from $\tilde{\Omega}_{s}$ to $\mathbb{R}$. One can observe that $\mathrm{div}(R^{s})^{*}\mathbf{p}_{s}\in\tilde{V}_{s}$ for all $\mathbf{p}_{s}\in W^{s}$. Thus, it is natural to define the local divergence operator $\mathrm{div}_{s}$: $W^{s}\rightarrow\tilde{V}^{s}$ by

[TABLE]

and the local gradient operator $\nabla_{s}$: $\tilde{V}^{s}\rightarrow W^{s}$ by $\nabla_{s}=-\mathrm{div}_{s}^{*}$. With these local operators, it is straightforward to obtain the equivalent primal and primal-dual form of Eq. 35 as follows (see, for example, [21] for details):

[TABLE]

Consequently, existing state-of-the-art solvers for the ROF model can be adopted for local problems of the proposed methods. For instance, one may solve Eq. 35 by FISTA [2], Eq. 37a by the augmented Lagrangian method [27], or Eq. 37b by the first order primal-dual algorithm [6].

Next, we consider the coloring technique used in the proposed methods. As we mentioned above, we assign the same color to subdomains if local problems associated to them are solved independently. To compute $g\in\tilde{V}^{s}$ in Eq. 35, we require the values of $\mathbf{q}\in W$ in the area marked in Fig. 3, i.e., $\tilde{\Omega}_{s}$ plus its adjacent lines of pixels on the top and the left sides. It means that the local problem Eq. 35 on $\Omega_{ij}$ is dependent on $\Omega_{i-1,j}$, $\Omega_{i-1,j+1}$, $\Omega_{i,j-1}$, $\Omega_{i,j+1}$, $\Omega_{i+1,j-1}$, and $\Omega_{i+1,j}$. Only $\Omega_{i-1,j-1}$ and $\Omega_{i+1,j+1}$ can have the same color as $\Omega_{ij}$ among its neighboring subdomains. It can be accomplished with 3 colors; each subdomain $\Omega_{ij}$ is colored with the color $(i-j)\textrm{ mod }3$. See Fig. 4 for visual description. We notice that the same coloring technique was introduced in [17] for primal DDMs.

We conclude the section with a remark on domain decomposition shapes. Since DDMs are designed for use on distributed memory computers, it is important to reduce the amount of communications among processors. In this sense, the window-shaped domain decomposition shown in Fig. 2(a) is most efficient because the total length of the subdomain interfaces is minimized. Moreover, we observed in Eq. 23 that the convergence rate of block Jacobi methods rely on the length of the subdomain interfaces. However, in the case when we have only a few number of processors, the stripe-shaped domain decomposition shown in Fig. 2(b) can be a good alternative. Since the stripe-shaped domain decomposition can be colored with only 2 colors, one may expect better convergence rate; recall that all the convergence theorems presented in this paper are dependent on the number of colors $N_{c}$. Numerical comparison between two types of domain decomposition will be presented in the next section.

5 Numerical experiments

In this section, we present numerical results for the block Jacobi methods introduced in this paper. All codes were programmed using ANSI C with OpenMPI and compiled by Intel Parallel Studio XE. All computations were performed on a computer cluster composed of seven machines, where each machine is equipped with two Intel Xeon SP-6148 CPUs (2.4GHz, 20C), 192GB RAM, and the operating system CentOS 7.4 64bit. Two test images “Peppers $512\times 512$” and “Boat $2048\times 3072$” are used in our numerical experiments. Both are corrupted with additive Gaussian noise with mean [math] and variance $0.05$; see Fig. 5. As a measurement of the quality of denoising, we use the peak-signal-to-noise ratio (PSNR) defined by

[TABLE]

where $\mathrm{MAX}=1$ is the possible maximum value of pixel intensity and $f_{\textrm{orig}}$ is the original clean image. The weight parameter $\alpha$ in Eq. 6 is chosen as 5, 10, and 20 in our experiments.

For all experiments, we set the initial guess as $\mathbf{p}^{(0)}=\mathbf{0}$. To reduce the time elapsed in solving local problems, the local solutions $\mathbf{p}_{s}^{(n-1)}$ from the previous iteration are chosen as initial guesses for the local problems at each iteration.

First, we compare the energy decay of several block methods: relaxed block Jacobi (Algorithm 1, RJ), block Gauss–Seidel (Algorithm 1 in [14] , GS), pre-relaxed block Jacobi (Algorithm 2, PJ), and fast pre-relaxed block Jacobi (Algorithm 3, FPJ) methods. Figure 6 plots $\frac{F(\mathbf{p}^{(n)})-F(\mathbf{p}^{*})}{F(\mathbf{p}^{*})}$ of the block methods per iteration in $\log$-$\log$ scale, where $\mathbf{p}^{*}$ is a solution of the full-dimension problem Eq. 6 obtained by $10^{6}$ FISTA [2] iterations. The number of subdomains $\mathcal{N}$ is fixed as $8\times 8$. Local problems were solved by FISTA with the stop criterion

[TABLE]

where the operator $\mathrm{div}_{s}$ was defined in Eq. 36. The above criterion was chosen so that stagnations of the energy do not occur (see Figure 5 in [19]). As we justified theoretically in Section 3, the energy functional of FPJ converges to the minimum much faster than other methods. As noted in Remark 3.23, such fast convergence occurs regardless of the given image $f$ and the weight parameter $\alpha$. In particular, the convergence rate of FPJ is faster than that of GS, while other block Jacobi methods converge more slowly than GS. We notice that the $O(1/n)$ convergence of GS can be deduced as a direct consequence of [7, 23]. It is observed that the energy decay of FPJ is not monotone. This is due to the nature of FISTA, and one may use MFISTA [1] instead of FISTA acceleration to ensure monotone decay of the energy. Note that the monotonicity of MFISTA given in Appendix A of [1] can be applied to our framework directly. Interestingly, as we mentioned in Remark 3.20, PJ seems to be faster than RJ at first 20 iterations, but eventually shows similar convergence results compared to RJ in our numerical experiments. Meanwhile, as shown in Fig. 7, all the results from three block Jacobi methods with $\alpha=10$ have the same PSNR and are not visually distinguishable. The same applies to the cases $\alpha=5,20$, and we omit the results for those cases.

For practical use in image restoration, the stop criteria need not be too strict; see [8] for details. In the following experiments, we use

[TABLE]

for outer iterations, and

[TABLE]

for local problems.

We compare the results of FPJ with respect to the weight parameter $\alpha$. Figure 8 shows the results of FPJ with $\alpha=5,10,$ and $20$ with $\mathcal{N}=8\times 8$. As it is well-known for the ROF model, the staircase effect becomes strong as $\alpha$ decreases. Among three results, $\alpha=10$ gives the highest PSNR. We observe from Table 2 that both the number of FPJ iterations and the wall-clock time increase as $\alpha$ decreases.

Table 3 shows the performance of FPJ with various numbers of subdomains $\mathcal{N}$ when $\alpha=10$. The case $\mathcal{N}=1$ presents the results of FISTA applied to the full-dimension problem for Eq. 6 using the stopping criterion Eq. 38. In the case of the test image “Peppers $512\times 512$,” the wall-clock time of FPJ decreases as $\mathcal{N}$ grows from 1 to $8\times 8$, but increases when $\mathcal{N}$ becomes $16\times 16$. It is because the size of “Peppers $512\times 512$” is so small that the portion of communication time between processors in the wall-clock time is not negligible compared to computation time if we use $\mathcal{N}=16\times 16$ subdomains. For a sufficiently large test image “Boat $2048\times 3072$,” the wall-clock time of FPJ decreases as $\mathcal{N}$ increases. It shows the efficiency of FPJ as a parallel solver for large-scale images. We also observe numerically that the number of FPJ iterations increases as $\mathcal{N}$ grows. Figure 9 shows the denoised test images by FPJ with $\mathcal{N}=16\times 16$. As they show no trace of the subdomain interfaces and have the same PSNR as the case $\mathcal{N}=1$, we can say that FPJ correctly solves Eq. 6 and recovers images.

Finally, we compare the convergence rate of FPJ for the window-shaped domain decomposition and the stripe-shaped one. Figure 10 shows the decay of $\frac{F(\mathbf{p}^{(n)})-F(\mathbf{p}^{*})}{F(\mathbf{p}^{*})}$ for both domain decompositions when $\mathcal{N}=4$ and $256$. Since $N_{c}=2$ for the stripe shape and $N_{c}=3$ for the window shape, the stripe-shaped domain decomposition converges faster than the window-shaped one when $\mathcal{N}=4$. However, since the total lengths of the subdomain interfaces are $O(\mathcal{N}^{1/2})$ and $O(\mathcal{N})$ for the window and stripe shape, respectively, as we observed in Remark 3.13, the window-shaped domain decomposition is faster than the stripe-shaped one when $\mathcal{N}$ is large enough; see Fig. 10(c) and (d). Table 4 provides the performance of FPJ with the stripe-shaped domain decomposition with various values of $\mathcal{N}$ when $\alpha=10$. Comparing Tables 3 and 4, it can be inferred that if $\mathcal{N}$ is small, the stripe-shaped domain decomposition is more efficient than the window-shaped one, and if $\mathcal{N}$ is large enough, the opposite holds.

6 Conclusion

We proved the $O(1/n)$ convergence rate of the nonoverlapping relaxed block Jacobi method for the dual ROF model proposed in [14]. Then, incorporating the FISTA acceleration technique [2] into this method, we proposed the $O(1/n^{2})$ convergent fast pre-relaxed block Jacobi method. We implemented various block decomposition methods and showed that the fast pre-relaxed block Jacobi method outperforms others numerically. The fast pre-relaxed block Jacobi method also showed good performance as a parallel solver for large-scale images.

Since our convergence analysis cannot be applied directly to overlapping cases, we plan to investigate ways to extend our acceleration methodology to overlapping methods [8] in our future work. In addition, to design scalable DDMs, construction of the coarse grid correction will be considered.