# A Preconditioned Alternating Minimization Framework for Nonconvex and Half
Quadratic Optimization

Shengxiang Deng Ismail Ben Ayed Hongpeng Sun Institute for Mathematical
Sciences, Renmin University of China, China. Email<EMAIL_ADDRESS>de Technolgie Superieure de Montréal, Canada. Email:
<EMAIL_ADDRESS>for Mathematical Sciences, Renmin
University of China, China. Email<EMAIL_ADDRESS>
###### Abstract

For some typical and widely used non-convex half-quadratic regularization
models and the Ambrosio-Tortorelli approximate Mumford-Shah model, based on
the Kurdyka-Łojasiewicz analysis and the recent nonconvex proximal algorithms,
we developed an efficient preconditioned framework aiming at the linear
subproblems that appeared in the nonlinear alternating minimization procedure.
Solving large-scale linear subproblems is always important and challenging for
lots of alternating minimization algorithms. By cooperating the efficient and
classical preconditioned iterations into the nonlinear and nonconvex
optimization, we prove that only one or any finite times preconditioned
iterations are needed for the linear subproblems without controlling the error
as the usual inexact solvers. The proposed preconditioned framework can
provide great flexibility and efficiency for dealing with linear subproblems
and guarantee the global convergence of the nonlinear alternating minimization
method simultaneously.

#### Key words.

Alternating minimization, half quadratic, nonconvex optimization, Kurdyka-
Łojasiewicz property, linear preconditioner

#### AMS subject classifications.

65K10, 90C25, 90C26, 65F08

##  1 Introduction

The aim of this paper is to develop a preconditioned framework to deal with
linear subproblems for the nonlinear and nonconvex alternating minimization
algorithms, while applying to some nonconvex half quadratic regularized
problems or the Ambrosio-Tortorelli approximate Mumford-Shah model [2, 42].
For the half quadratic models, we mainly focus on the truncated quadratic
model, the Geman-McClure model, and the Hebert-Leahy model [28, 29, 27, 31].
The truncated quadratic models including the Geman-Reynolds type model [28]
and Geman-Yang type model [29] are rooted from the Markov random fields. All
the half quadratic models addressed here have statistical interpretations and
we refer to [9, 30, 53] for more details. These nonconvex regularizations are
widely used for image restorations, segmentation, stereo, optical flow, and so
on, which have vast applications in medical imaging, computer vision, and
inverse problems [8, 9, 10, 53, 32, 42].

Due to the importance of these models, there are many theoretical and
algorithmic studies on each of these models [6, 10, 21, 22, 20, 23, 33, 37,
43, 44, 48]. For the truncated quadratic model, the graduated non-convexity
algorithm was proposed in [10] and see [18, 43] for recent developments.
Difference of convex algorithm (DCA) [37], preconditioned DCA [25], and the
first-order primal-dual algorithm [48] are also developed for the truncated
quadratic models. For the German-McClure model, Hebert-Leahy model, and
Ambrosio-Tortorelli model, there are also lots of analysis and algorithmic
developments; see [6, 8, 9, 21, 32] and so on.

Now, let’s restrict our attention on the alternating minimization (AM)
optimization. The alternating minimization algorithms are essentially the same
as the majorize-minimize type algorithm including the expectation–maximization
(EM) algorithm and the bound optimization [1, 23, 26, 19, 35, 45, 49]. The
general alternating minimization method for $L(\boldsymbol{u},\boldsymbol{y})$
reads as follows

$\inf_{\boldsymbol{u}\in X,\boldsymbol{y}\in
Y}L(\boldsymbol{x},\boldsymbol{y}):=f(\boldsymbol{u})+Q(\boldsymbol{u},\boldsymbol{y})+g(\boldsymbol{y}).$
(1.1)

Given initial value $\boldsymbol{y}^{0}$ (or $\boldsymbol{x}^{0}$), iterate
for $k=0,1,\cdots$, until some stopping criterion is fulfilled

$\displaystyle\boldsymbol{u}^{k+1}$
$\displaystyle\in\operatorname*{arg\,min}_{\boldsymbol{u}\in
X}L(\boldsymbol{u},\boldsymbol{y}^{k}),$ (1.2a)
$\displaystyle\boldsymbol{y}^{k+1}$
$\displaystyle\in\operatorname*{arg\,min}_{\boldsymbol{y}\in
Y}L(\boldsymbol{u}^{k+1},\boldsymbol{y}),$ (1.2b)

where $f$ and $g$ are proper lower semicontinuous functions, $Q$ is a $C^{1}$
function with local Lipschitz continuous gradient $\nabla Q$, and $X$, $Y$ are
the corresponding finite dimensional Hilbert spaces. Throughout this paper, we
assume all the variables, the spaces, and the operators including the
integrals setting are all finite dimensional. For the half quadratic models,
(1.1) can be designed for the following equivalent minimization problem, where
the pioneering work can be find in [28, 29]:

$\inf_{\boldsymbol{u}\in X}F(\boldsymbol{u}),\quad
F(\boldsymbol{u}):=D(\boldsymbol{u})+R(\boldsymbol{u})=\inf_{\boldsymbol{y}\in
Y}L(\boldsymbol{u},\boldsymbol{y}).$ (1.3)

Here $D(\boldsymbol{u})$ is the data term and $R(\boldsymbol{u})$ is the
regularization term. Instead of minimizing the original functional (1.3), it
is convenient to minimizing (1.1) with the auxiliary variable
$\boldsymbol{y}$. In discrete optimization including the binary optimization,
the auxiliary variable $\boldsymbol{y}$ is called “line process” in [10, 28,
29] or more general “outlier process” [9]. These line or outlier processes are
popular and widely used, not only because they have physical or statistic
intuition and the ability to model spatial properties of discontinuities [9]
but also they can make the optimization procedure more efficient and stable.
Actually, one can get (1.1) from the original (1.3) by the Fenchel-Rockafellar
duality theory. We refer to the Geman-Reynolds model [28] and Geman-Yang model
[29] which are two typical and different duality processes to get (1.1) from
(1.3). In [6], the alternating minimization focused on convex regularization
term $R(\boldsymbol{u})$ is discussed. In [21, 22, 1] and the thesis [20], two
alternating minimization algorithms including “ARTUR” and “LEGEND” algorithms
were developed with global convergence analysis for convex models with convex
regularizations. Although there are rich studies on alternating minimization
for the half quadratic models or the Ambrosio-Tortorelli model, however,
solving the large scale linear subproblems for $\boldsymbol{u}$ (or
$\boldsymbol{y}$) is still very challenging. Recent attempts can be found in
[23, 35] where conjugate gradient method with line search as an inexact solver
was developed. Based on the weighted quadratic approximations on
$F(\boldsymbol{u})$, truncated conjugate gradient method with error control
was also developed for $C^{1}$ function $F(\boldsymbol{u})$ in [45].

Inspired by the novel preconditioning techniques developed for linear
subproblems that appeared in nonlinear convex optimizations [15, 14], the new
development of Kurdyka-Łojasiewicz (KL) analysis, and the proximal nonconvex
optimization [3, 5, 4, 13], we proposed a preconditioned framework for the
nonconvex alternating minimization aiming at dealing with the linear problems
especially large-scale problems efficiently for (1.1). Our contributions
belong to the following folds. First, we proposed a preconditioned framework
that can deal with any linear subequation of $\boldsymbol{u}$ or
$\boldsymbol{y}$. Any finite number of preconditioned iterations can guarantee
the convergence of the whole nonlinear alternating minimization algorithm
without solving the linear subequations with middle or high accuracy, which is
different from the inexact solvers through error control. Especially, for the
Ambrosio-Tortorelli model, preconditioned iterations can be employed for both
linear equations involving $\boldsymbol{u}$ and $\boldsymbol{y}$. Second, with
the analysis of the semialgebraic sets and o-minimal structure, we prove that
all the functions $L(\boldsymbol{u},\boldsymbol{y})$ of discussed models are
KL functions. Together with the boundedness of the iteration sequence, we
obtain the global convergence of the proposed preconditioned alternating
minimization by the recent developments of proximal minimization algorithms
[3, 5, 4, 13]. Third, our preconditioned framework can deal with the nonsmooth
and nonconvex truncated quadratic problems efficiently with global
convergence, while they are precluded by convex assumptions as in [6, 1, 22]
or $C^{1}$ conditions as in [45]. Fourth, we developed several efficient red-
black Gauss-Seidel preconditioners for both isotropic and anisotropic
equations in divergence form by finite difference method [39, 51]. Numerical
tests show that one can get lower energy and better reconstructions more
efficiently with the proposed preconditioned framework compared with solving
the linear system of (1.2) with moderate accuracy and without any proximal
terms.

The rest of this paper is organized as follows. In section 2, we give a brief
introduction to the half quadratic models, the Ambrosio-Tortorelli model, and
the basic KL analysis. In section 3, we give an illustration of the motivation
of our preconditioned framework and prove that all
$L(\boldsymbol{u},\boldsymbol{y})$ of discussed models are KL functions. In
section 4, we first prove the iterative sequence
$\\{(\boldsymbol{u}^{k},\boldsymbol{y}^{k})\\}$ is bounded for each model. We
then get the global convergence with the KL properties by [3, 5, 4, 13]. In
section 5, we first give the detailed five-point stencils of the symmetric
Gauss-Seidel iterations along with the preconditioners and the numerical
experiments then follows. In section 6, we give some discussions and
conclusions.

##  2 Nonconvex half quadratic models and KL functions

For the half quadratic models, let’s begin with the following truncated
quadratic regularization with gradient operator
$\nabla:=[\nabla_{1},\nabla_{2}]$ for image restoration [10, 9, 42, 43]

$\displaystyle\inf_{\boldsymbol{u}\in
X}F(\boldsymbol{u})=D(\boldsymbol{u})+\int_{\Omega}P^{I}(\boldsymbol{u})\mathrm{d}\sigma,\quad
P^{I}(\boldsymbol{u}):=\frac{\mu}{2}\min(|\nabla\boldsymbol{u}|^{2},\frac{\lambda}{\mu}),$
(2.1) $\displaystyle\inf_{\boldsymbol{u}\in
X}F(\boldsymbol{u})=D(\boldsymbol{u})+\int_{\Omega}P^{A}(\boldsymbol{u})\mathrm{d}\sigma,\quad
P^{A}(\boldsymbol{u}):=\frac{\mu}{2}\sum_{i=1}^{2}\min(|\nabla_{i}\boldsymbol{u}|^{2},\frac{\lambda}{\mu}),$
(2.2)

where $\lambda$ and $\mu$ are positive constants, $X$ is finite dimensional
discrete image space and henceforth $D(\boldsymbol{u})$ is the following data
term with $A$ being a linear and bounded operator

$D(\boldsymbol{u}):={\|A\boldsymbol{u}-\boldsymbol{u}_{0}\|_{2}^{2}}/{2}.$

$P^{I}$ or $P^{A}$ is the isotropic or anisotropic truncated regularizations.
Let’s introduce the image space $X$ and the dual spaces $Y$ and $Z$, i.e.,

$X=\\{\boldsymbol{u}:\Omega\rightarrow\mathbb{R}\\},\quad
Y=\\{\boldsymbol{u}:\Omega\rightarrow\mathbb{R}^{2}\\},\quad
Z=\\{\boldsymbol{u}:\Omega\rightarrow\mathbb{R}\\},$ (2.3)

where $\Omega$ is the image domain. The truncated quadratic model (2.1) can be
reformulated as the following Geman-Reynolds model [28],

$\inf_{\boldsymbol{u}\in X,\boldsymbol{b}\in
Z}L_{GR}(\boldsymbol{u},\boldsymbol{b}):=D(\boldsymbol{u})+\frac{\lambda}{2}\int_{\Omega}\left[\boldsymbol{b}(\frac{\mu}{\lambda}|\nabla\boldsymbol{u}|^{2}-1)+I_{\\{0\leq\boldsymbol{b}\leq
1\\}}(\boldsymbol{b})+1\right]\mathrm{d}\sigma,$ (2.4)

or the following Geman-Yang model [29, 42]: $\inf_{\boldsymbol{u}\in
X,\boldsymbol{l}\in Y}L_{GY}(\boldsymbol{u},\boldsymbol{l})$,
$\boldsymbol{l}=(\boldsymbol{l}_{1},\boldsymbol{l}_{2})^{T}$ with
$L_{GY}(\boldsymbol{u},\boldsymbol{l})$ defined by,

$D(\boldsymbol{u})+\frac{\mu}{2}\|\nabla\boldsymbol{u}-\boldsymbol{l}\|_{2}^{2}+\mu\int_{\Omega}H(\boldsymbol{l};\frac{\lambda}{\mu})\mathrm{d}\sigma,\
H(\boldsymbol{l};\frac{\lambda}{\mu}):=\begin{cases}\sqrt{\frac{\lambda}{\mu}}|\boldsymbol{l}|-\frac{|\boldsymbol{l}|^{2}}{2},\
|\boldsymbol{l}|\leq\sqrt{\frac{\lambda}{\mu}},\\\ \frac{\lambda}{2\mu},\
\quad|\boldsymbol{l}|>\sqrt{\frac{\lambda}{\mu}}.\end{cases}$ (2.5)

For the anisotropic case, similarly, we can reformulate (2.2) as

$\inf_{\boldsymbol{u}\in X,\boldsymbol{b}\in
Y}L_{GR}^{A}(\boldsymbol{u},\boldsymbol{b}):=D(\boldsymbol{u})+\frac{\lambda}{2}\sum_{i=1}^{2}\int_{\Omega}\left[\boldsymbol{b}_{i}(\frac{\mu}{\lambda}|\nabla_{i}\boldsymbol{u}|^{2}-1)+I_{\\{0\leq\boldsymbol{b}_{i}\leq
1\\}}(\boldsymbol{b}_{i})+1\right]\mathrm{d}\sigma,$ (2.6)

and

$\inf_{\boldsymbol{u}\in X,\boldsymbol{l}\in
Y}L_{GY}^{A}(\boldsymbol{u},\boldsymbol{l}):=D(\boldsymbol{u})+\sum_{i=1}^{2}\big{(}\frac{\mu}{2}\|\nabla_{i}\boldsymbol{u}-\boldsymbol{l}_{i}\|_{2}^{2}+\mu\int_{\Omega}H(\boldsymbol{l}_{i};\frac{\lambda}{\mu})\mathrm{d}\sigma\big{)}.$
(2.7)

The following German-McClure model is also widely used [27, 28]

$\inf_{\boldsymbol{u}\in
X}F_{GM}(\boldsymbol{u}):=D(\boldsymbol{u})+\frac{\mu}{2}\int_{\Omega}\frac{|\nabla\boldsymbol{u}|^{2}/\lambda}{1+|\nabla\boldsymbol{u}|^{2}/\lambda}\mathrm{d}\sigma,$
(2.8)

which is equivalent to the following minimization problem [9, 20]

$\inf_{\boldsymbol{u}\in X,\boldsymbol{b}\in
Z}L_{GM}(\boldsymbol{u},\boldsymbol{b}):=D(\boldsymbol{u})+\frac{\mu}{2}\int_{\Omega}\left[\boldsymbol{b}\frac{|\nabla\boldsymbol{u}|^{2}}{\lambda}+\boldsymbol{b}-2\sqrt{\boldsymbol{b}}+1\right]\mathrm{d}\sigma,\quad\boldsymbol{b}\geq
0.$ (2.9)

Similarly, the anisotropic German-McClure model can be

$\inf_{\boldsymbol{u}\in X,\boldsymbol{b}\in
Y}L_{GM}(\boldsymbol{u},\boldsymbol{b}):=D(\boldsymbol{u})+\sum_{i=1}^{2}\frac{\mu}{2}\int_{\Omega}\left[\boldsymbol{b}_{i}\frac{|\nabla_{i}\boldsymbol{u}|^{2}}{\lambda}+\boldsymbol{b}_{i}-2\sqrt{\boldsymbol{b}_{i}}+1\right]\mathrm{d}\sigma,$
(2.10)

where $\boldsymbol{b}=(\boldsymbol{b}_{1},\boldsymbol{b}_{2})$ and
$\boldsymbol{b}_{i}\geq 0$, $i=1,2$. Now, let’s turn to the Hebert-Leahy model
which reads as follows [31]

$\inf_{\boldsymbol{u}\in
X}F_{HL}(\boldsymbol{u}):=D(\boldsymbol{u})+\frac{\mu}{2}\int_{\Omega}\log(1+\frac{|\nabla\boldsymbol{u}|^{2}}{\lambda})\mathrm{d}\sigma,$
(2.11)

which is equivalent to the following minimization problem [8, 9, 20]

$\inf_{\boldsymbol{u}\in X,\boldsymbol{b}\in
Z}L_{HL}(\boldsymbol{u},\boldsymbol{b}):=D(\boldsymbol{u})+\frac{\mu}{2}\int_{\Omega}(\boldsymbol{b}\cdot\frac{|\nabla\boldsymbol{u}|^{2}}{\lambda}+\boldsymbol{b}-\log(\boldsymbol{b})-1)\mathrm{d}\sigma,\quad\boldsymbol{b}>0.$
(2.12)

Similarly, the anisotropic Hebert-Leahy model can be

$\inf_{\boldsymbol{u}\in X,\boldsymbol{b}\in
Y}L_{HL}(\boldsymbol{u},\boldsymbol{b}):=D(\boldsymbol{u})+\sum_{i=1}^{2}\frac{\mu}{2}\int_{\Omega}(\boldsymbol{b}_{i}\cdot\frac{|\nabla_{i}\boldsymbol{u}|^{2}}{\lambda}+\boldsymbol{b}_{i}-\log(\boldsymbol{b}_{i})-1)\mathrm{d}\sigma,$
(2.13)

where $\boldsymbol{b}=(\boldsymbol{b}_{1},\boldsymbol{b}_{2})$ and
$\boldsymbol{b}_{i}\geq 0$, $i=1,2$. We will also discuss the following
Ambrosio–Tortorelli approximation of the Mumford-Shah model [2]

$\displaystyle\inf_{\boldsymbol{u}\in X,\boldsymbol{s}\in
Z}L_{MS}(\boldsymbol{u},\boldsymbol{s}):=$
$\displaystyle\frac{\|\boldsymbol{u}-\boldsymbol{u}_{0}\|_{2}^{2}}{2}+\alpha\int_{\Omega}\boldsymbol{s}^{2}|\nabla\boldsymbol{u}|^{2}\mathrm{d}\sigma$
$\displaystyle+\lambda\left(\varepsilon\int_{\Omega}|\nabla\boldsymbol{s}|^{2}\mathrm{d}\sigma+\int_{\Omega}\frac{1}{4\varepsilon}|\boldsymbol{s}-1|^{2}\mathrm{d}\sigma\right).$
(2.14)

Based on difference of convex functions, we will give another interpretation
of the equivalence of the truncated model (2.1) and (2.5) compared to [29,
42]. The equivalence of the anisotropic cases (2.2) and (2.7) is similar and
omitted here. Let’s begin with the following lemma.

######  Lemma 1.

Denote $h_{a}(\boldsymbol{t}):=\max(|\boldsymbol{t}|^{2}/{2},{a}/{2})$ which
is a convex function with positive constant $a$ and $\boldsymbol{t}\in Y$. The
Fenchel dual function of $h_{a}(\boldsymbol{t})$ is

$h_{a}^{*}(\boldsymbol{l})=\begin{cases}\sqrt{a}|\boldsymbol{l}|-a/2,\quad|\boldsymbol{l}|\leq\sqrt{a},\\\
{|\boldsymbol{l}|^{2}}/{2},\quad|\boldsymbol{l}|>\sqrt{a}.\end{cases}$ (2.15)

###### Proof.

It can be checked that $|\boldsymbol{l}||\boldsymbol{t}|\leq
h_{a}(\boldsymbol{t})$ when $|\boldsymbol{l}|\leq\sqrt{a}$, since the curve
$|\boldsymbol{l}||\boldsymbol{t}|$ is below
$\max(|\boldsymbol{t}|^{2}/{2},{a}/{2})$ as functions of $|\boldsymbol{t}|$
for fixed $\boldsymbol{l}$. While $|\boldsymbol{l}|\leq\sqrt{a}$, we thus have

$\displaystyle\sup_{\boldsymbol{t}\in Y}$
$\displaystyle\left\\{\langle\boldsymbol{t},\boldsymbol{l}\rangle-\max({|\boldsymbol{t}|^{2}}/{2},{a}/{2})\right\\}=\sup_{|\boldsymbol{t}|\in[0,\sqrt{a}]}\left\\{\langle\boldsymbol{t},\boldsymbol{l}\rangle-\max({|\boldsymbol{t}|^{2}}/{2},{a}/{2})\right\\}$
$\displaystyle=\sup_{|\boldsymbol{t}|\in[0,\sqrt{a}]}\left\\{\langle\boldsymbol{t},\boldsymbol{l}\rangle-{a}/{2}\right\\}=\sup_{|\boldsymbol{t}|=\sqrt{a}}\left\\{\langle\boldsymbol{t},\boldsymbol{l}\rangle-{a}/{2}\right\\}=|\boldsymbol{l}|\sqrt{a}-{a}/{2}.$

While $|\boldsymbol{l}|>\sqrt{a}$, it can be verified that
$|\boldsymbol{l}||\boldsymbol{t}|=h_{a}(t)$ while
$|\boldsymbol{t}|=\frac{a}{2|\boldsymbol{l}|}$ or $2|\boldsymbol{l}|$. We thus
conclude that, while $|\boldsymbol{l}|>\sqrt{a}$,

$\displaystyle\sup_{\boldsymbol{t}\in Y}$
$\displaystyle\left\\{\langle\boldsymbol{t},\boldsymbol{l}\rangle-\max({|\boldsymbol{t}|^{2}}/{2},{a}/{2})\right\\}=\sup_{|\boldsymbol{t}|\in[0,2|\boldsymbol{l}|]}\left\\{\langle\boldsymbol{t},\boldsymbol{l}\rangle-\max({|\boldsymbol{t}|^{2}}/{2},{a}/{2})\right\\}$
$\displaystyle=\sup_{|\boldsymbol{t}|\in[\sqrt{a},2|\boldsymbol{l}|]}\left\\{\langle\boldsymbol{t},\boldsymbol{l}\rangle-{|\boldsymbol{t}|^{2}}/{2}\right\\}=\left\\{\langle\boldsymbol{t},\boldsymbol{l}\rangle-{|\boldsymbol{t}|^{2}}/{2}\right\\}|_{{|\boldsymbol{t}|=|\boldsymbol{l}|}}={|\boldsymbol{l}|^{2}}/{2}.$

With the discussion of the two cases above, we get (2.15). ∎

######  Lemma 2.

The truncated models (2.1) and (2.5) are equivalent.

###### Proof.

For (2.1), it can be written as the difference of convex functions (DCA):
$F(\boldsymbol{u})=F_{1}(\boldsymbol{u})-F_{2}(\boldsymbol{u})$ with convex
$F_{1}$ and $F_{2}$

$F_{1}(\boldsymbol{u}):=D(\boldsymbol{u})+\int_{\Omega}(\frac{\mu}{2}|\nabla\boldsymbol{u}|^{2}+\frac{\lambda}{\mu})d\sigma,\quad
F_{2}(\boldsymbol{u})=\int_{\Omega}\frac{\mu}{2}\max(|\nabla\boldsymbol{u}|^{2},\frac{\lambda}{\mu})\mathrm{d}\sigma.$
(2.16)

By Lemma 1, we have

$\frac{\mu}{2}\max(|\nabla\boldsymbol{u}|^{2},\frac{\lambda}{\mu})={\mu}\max(|\nabla\boldsymbol{u}|^{2}/2,\frac{\lambda}{\mu}/2)=\mu\sup_{\boldsymbol{l}\in
Y}\left\\{\langle\nabla\boldsymbol{u},\boldsymbol{l}\rangle-
h_{\frac{\lambda}{\mu}}^{*}(\boldsymbol{l})\right\\}.$ (2.17)

Substituting (2.17) into (2.16), we get this lemma with direct calculation. ∎

Henceforth, we will vectorize all the image variables $\boldsymbol{u}$,
$\boldsymbol{u}_{0}$ and the auxiliary variable $\boldsymbol{y}$ including
$\boldsymbol{l}$, $\boldsymbol{b}$ and $\boldsymbol{s}$ as
$u\in\mathbb{R}^{MN}$, $u_{0}\in\mathbb{R}^{M_{0}N_{0}}$ and $y$ including
${l}$, ${b}$ and ${s}$ in the corresponding spaces and for the corresponding
models. Besides, we will still use the linear operators $A$, $\nabla_{i}$,
$i=1,2$ as their corresponding matrix versions after vectorizing all the
variables. Now, we will touch some necessary tools from convex and variational
analysis [41, 24, 46]. The graph of a multivalued mapping $F:R^{n}\Rightarrow
R^{m}$ is defined by

$\operatorname{gph}F:=\\{(x,y)\in R^{n}\times R^{m}:y\in F(x)\\},$

whose domain is defined by $\operatorname{dom}F:=\\{x\ |\
F(x)\neq\emptyset\\}$. Similarly the graph of an extended real-valued function
$f:R^{n}\to R\cup\\{+\infty\\}$ is defined by

$\operatorname{gph}f:=\\{(x,s)\in R^{n}\times R:s=f(x)\\}.$

Let $h:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ be a proper lower
semicontinuous function. Denote $\operatorname{dom}h:=\\{x\in\mathbb{R}^{n}:\
h(x)<+\infty\\}$. For each $x\in\operatorname{dom}f$, the limiting-
subdifferential of $h$ at $x\in\mathbb{R}^{n}$, written $\partial f$, is
defined as follows [41, 46],

$\displaystyle\partial h(x):=\bigg{\\{}$
$\displaystyle\xi\in\mathbb{R}^{n}:\exists x_{n}\rightarrow
x,h(x_{n})\rightarrow h(x),\xi_{n}\rightarrow\xi,$
$\displaystyle\lim_{y\rightarrow x}\inf_{y\neq
x_{n}}\frac{h(y)-h(x_{n})-\langle\xi_{n},y-x_{n}\rangle}{|y-x_{n}|}\geq
0\bigg{\\}}.$

It is known that the above subdifferential $\partial h$ reduces to the
classical subdifferential in convex analysis when $h$ is convex. It can be
seen that a necessary condition for $x\in\mathbb{R}^{n}$ to be a minimizer of
$h$ is $0\in\partial h$ [3].

For the global and local convergence analysis, we also need the Kurdyka-
Łojasiewicz (KL) property and KL exponent. While the KL properties can help
obtain the global convergence of iterative sequences, the KL exponent can help
provide a local convergence rate.

######  Definition 1 (KL property and KL exponent).

A proper closed function $h$ is said to satisfy the KL property at
$\bar{x}\in\operatorname{dom}\partial h$ if there exists $a\in(0,+\infty]$, a
neighborhood $\mathcal{O}$ of $\bar{x}$, and a continuous concave function
$\psi:[0,a)\rightarrow(0,+\infty)$ with $\psi(0)=0$ such that:

  * _(i)_

$\psi$ is continuous differentiable on $(0,a)$ with $\psi^{\prime}>0$.

  * _(ii)_

For any $x\in\mathcal{O}$ with $h(\bar{x})<h(x)<h(\bar{x})+a$, one has

$\psi^{\prime}(h(x)-h(\bar{x}))\operatorname{dist}(0,\partial h(x))\geq 1.$
(2.18)

A proper closed function $h$ satisfying the KL property at all points in
$\operatorname{dom}\partial h$ is called a KL function. If $\psi$ in (2.18)
can be chosen as $\psi(s)=cs^{1-\theta}$ for some $\theta\in[0,1)$ and $c>0$,
we say that $h$ satisfies KL properties at $\bar{x}$ with exponent $\theta$.
This means that for some $\bar{c}>0$, we have

$\operatorname{dist}(0,\partial h(x))\geq\bar{c}(h(x)-h(\bar{x}))^{\theta}.$
(2.19)

If $h$ satisfies KL property with exponent $\theta\in[0,1)$ at all the points
of $\operatorname{dom}\partial h$, we call $h$ is a KL function with exponent
$\theta$.

For KL functions, the semialgebraic functions and definable functions in an
o-minimal structure provide a vast field of applications including the KL
analysis for our models to be discussed.

######  Definition 2 (Semialgebraic set and Semialgebraic function [4]).

A subset $S$ of $\mathbb{R}^{n}$ is called a real semialgebraic set if there
exists a finite number of real polynomial functions $P_{i,j},\
Q_{i,j}:\mathbb{R}^{n}\rightarrow\mathbb{R}$, such that

$S=\bigcup_{j=1}^{p}\bigcap_{j=1}^{q}\\{x\in\mathbb{R}^{n}:P_{i,j}=0,\
Q_{i,j}>0\\}.$

A function $f:\mathbb{R}^{n}\to\mathbb{R}\cup\\{+\infty\\}$ is semialgebraic
if its graph is a semialgebraic set of $\mathbb{R}^{n+1}$.

A very useful conclusion is that a semialgebraic function has the KL property
with $\psi(s)=cs^{1-\theta}$ for some
$\theta\in\left[0,1\right)\cap\mathbb{Q}$ and $c>0$, which can be seen as a
corollary of Theorem 3.2 of [11]. The following Tarski-Seidenberg theorem is
very useful for the analysis of KL properties.

######  Theorem 1 (Tarski-Seidenberg [4]).

Let $S$ be a semialgebraic set in $\mathbb{R}^{m+n}$, then

$\bar{S}:=\\{x\in\mathbb{R}^{m}:(x,y)\in S\ \text{for some}\
y\in\mathbb{R}^{n}\\}$

is a semialgebraic set.

For lots of cases, when the Tarski-Seidenberg theorem is not applicable, we
can turn to the o-minimal structure which is originated from real algebraic
geometry and can cover more complicated cases. We refer to [50, 5] for its
definition. Verifying the o-minimal structure directly with the definition is
much more complicated compared to verifying the semialgebraic set and we will
focus on the existed results that can be employed directly. Let $\Delta$ be an
o-minimal structure. A set $A$ is called definable if $A\in\Delta$. A map $f$
is said to be definable if its graph
$\operatorname{gph}f\subseteq\mathbb{R}^{m+n}$ is definable [50]. Due to their
dramatic impacts, these structures are being studied extensively. One of the
interests of such structures in optimization is due to the following nonsmooth
extension of KL property [5, 12] (see Theorem 11 of [12]).

######  Theorem 2 ([12]).

Any proper lower semicontinuous function
$f:\mathbb{R}^{n}\Rightarrow\mathbb{R}\cup\\{+\infty\\}$ that is definable in
an o-minimal structure $\Delta$ has the KL property at each point of
$\operatorname{dom}\partial f$. Moreover, the function $\psi$ is definable in
$\Delta$.

Let $R_{\exp}=(R,+,\cdot,\exp)$. Wilkie proved that $R_{\exp}$ is model
complete [52]. As a direct consequence of this theorem, each definable sets in
$R_{\exp}$ is the image of the zero set of a function in
$P(x,y,e^{x},e^{y})=0$ under a natural projection (see page 3 of [50] or
[36]). Then by a Khovanskii result on fewnomials [34], $R_{\exp}$ is an
o-minimal structure. An analytic proof of Wilkie’s theorem is given in [40].

######  Theorem 3 ([50, 36]).

The images in $\mathbb{R}^{n}$ for $n=0,1,2,\ldots$ under projection maps
$\mathbb{R}^{n+k}\Rightarrow\mathbb{R}^{n}$ of sets with the form
$\\{(x,y)\in\mathbb{R}^{n+k}:P(x,y,e^{x},e^{y})=0\\}$ is definable, where $P$
is a real polynomial in $2(n+k)$ variables with $x:=(x^{1},\ldots,x^{n})$,
$y:=(y^{1},\ldots,y^{k})$, $e^{x}:=(e^{x^{1}},\ldots,e^{x^{n}})$ and
$e^{y}:=(e^{y^{1}},\ldots,e^{y^{k}})$.

##  3 Preconditioned framework and KL properties

In this section, we will investigate the following preconditioned alternating
minimization framework (3.1) and give an illustration of our motivation
through the Lemma 3 to be discussed. The proposed framework includes the
updates of $u^{k+1}$ in (2.4), (2.5), (2.9), (2.12), (2.14) for the isotropic
cases and (2.6), (2.7), (2.10), (2.13) for the anisotropic cases along with
the update $s^{k+1}$ in (2.14). Our preconditioned framework for (1.1) is as
follows:

$\displaystyle u^{k+1}$
$\displaystyle\in\operatorname*{arg\,min}{L(u,y^{k})+\frac{1}{2}\|u-u^{k}\|_{M_{k}}^{2}},$
(3.1a) $\displaystyle y^{k+1}$
$\displaystyle\in\operatorname*{arg\,min}{L(u^{k+1},y)+\frac{1}{2}\|y-y^{k}\|_{N_{k}}^{2}},$
(3.1b)

where $M_{k},N_{k}$ are the proximal matrices satisfying [5, 4]

$0<\gamma_{-}I\leq M_{k}\leq\gamma_{+}I,\quad 0<\mu_{-}I\leq
N_{k}\leq\mu_{+}I,\quad 0<\gamma_{-},\ \mu_{-}<\gamma_{+},\ \mu_{+}<+\infty.$
(3.2)

Here $y^{k}$ denotes $b^{k}$, $l^{k}$ or $s^{k}$ in the corresponding models.
We will employ the metric induced by $M_{k}$ depending on the corresponding
model, which turns out to be the classical and powerful preconditioned
iterations. It can bring out flexibility and efficiency for linear
subproblems.

Henceforth we will focus on the denoising problems with $A=I$. We can
reformulate the original Euler-Lagrangian equation for $u^{k+1}$ (or $s^{k+1}$
for (2.14)) as the following general form

$\mathbb{T}_{k}u=\mathfrak{b}^{k}\ \text{in}\ \Omega,\quad\frac{\partial
u^{k+1}}{\partial\nu}|_{\partial\Omega}=0,\quad\mathbb{T}_{k}:=\gamma^{k}(x)I+\nabla^{*}\mathcal{B}^{k}\nabla,$
(3.3)

where $\gamma_{k}(x)\geq\bar{\gamma}$ for arbitrary $x\in\Omega$ with the
constant $\bar{\gamma}>0$ and the linear operator $\mathcal{B}^{k}$ is
positive semidefinite. We have $\mathbb{T}_{k}=I-\mu\Delta$,
$\mathfrak{b}^{k}=u_{0}+\nabla^{*}l^{k}$ with $\gamma^{k}(x)=1.0$ and
$\mathcal{B}^{k}=\text{Diag}[\mu I,\mu I]$ for (2.5),
$\mathbb{T}_{k}=I+\mu\nabla^{*}{b}^{k}\nabla$, $\mathfrak{b}^{k}=u_{0}$ with
$\gamma^{k}(x)=1.0$ and $\mathcal{B}^{k}=\text{Diag}[\mu{b}^{k},\mu{b}^{k}]$
for (2.4) and $\mathfrak{b}^{k}=u_{0}$ with $\gamma^{k}(x)=1.0$

$\mathbb{T}_{k}=I+\nabla^{*}(\mathcal{B}^{k}\nabla),\quad\mathcal{B}^{k}:=\text{Diag}[\mu/\lambda
b_{1}^{k},\mu/\lambda b_{2}^{k}],$ (3.4)

for (2.10) and (2.13). Similarly, for $s^{k+1}$ update in (2.14), we have
$\mathbb{T}_{k}=(\frac{\lambda}{2\varepsilon}+2\alpha|\nabla
u^{k}|^{2})I-\lambda\varepsilon\Delta$ with
$\gamma^{k}(x)=\frac{\lambda}{2\varepsilon}+2\alpha|\nabla u^{k}|^{2}$ and
$\mathcal{B}^{k}=\text{Diag}[\alpha\varepsilon,\alpha\varepsilon]$. All the
other cases are similar. Inspired by the recent development of preconditioning
technique for linear subproblems in the nonlinear convex [15, 14, 16] or
nonconvex iteration [25], we will introduce the classical preconditioned
iteration to deal with linear subproblems (3.3). Our motivation mainly comes
from the following lemma, i.e., Lemma 3.

######  Lemma 3.

With appropriately chosen linear operators $M_{k}$ and $N_{k}$ satisfying
(3.2), for the update of $u^{k+1}$ in the isotropic cases (2.4), (2.5), (2.9),
(2.12), (2.14) and the anisotropic cases (2.6), (2.7), (2.10), (2.13) along
with $s^{k+1}$ in (2.14), these updates can be written as the following
classical preconditioned iteration for the original equation (3.3),

$u^{k+1}:=u^{k}+\mathbb{M}_{k}^{-1}[\mathfrak{b}^{k}-\mathbb{T}_{k}u^{k}],\quad\mathbb{M}_{k}:=\mathbb{T}_{k}+M_{k}.$
(3.5)

###### Proof.

Suppose the original Euler-Lagrangian equation for the corresponding
functional is (3.3). With adding proximal term in (3.1a), the Euler-Lagrangian
equation for $u$ becomes

$\mathbb{T}_{k}u+M_{k}(u-u^{k})-\mathfrak{b}^{k}=0.$ (3.6)

Still with notation $u^{k+1}$ as the solved $u$ above, we thus arrive at

$\displaystyle u^{k+1}$
$\displaystyle=(M_{k}+\mathbb{T}_{k})^{-1}(\mathfrak{b}^{k}+M_{k}u^{k})$
(3.7a)
$\displaystyle=(M_{k}+\mathbb{T}_{k})^{-1}[(M_{k}+\mathbb{T}_{k})u^{k}+\mathfrak{b}^{k}-\mathbb{T}_{k}u^{k}]$
(3.7b)
$\displaystyle=u^{k}+\mathbb{M}_{k}^{-1}[\mathfrak{b}^{k}-\mathbb{T}_{k}u^{k}],$
(3.7c)

which is essentially the classical preconditioned iteration for solving the
linear equation $\mathbb{T}_{k}u=\mathfrak{b}^{k}$ [47]. We thus reformulate
the proximal iteration in (3.1) as the preconditioned iteration (3.5), which
will turn out very useful. ∎

Through Lemma 3, we introduce the classical preconditioned iterations from
numerical linear algebra and computation to the nonlinear alternating
minimization by specially designed positive definite proximal terms $M_{k}$
and $N_{k}$. The idea can also be found in [16, 25]. Here, for self-
completeness, we give some explanation through a concrete example.

######  Remark 1.

Suppose the discretization of the linear operator $\mathbb{T}_{k}$ in Lemma 3
is $D_{k}-E_{k}-E_{k}^{*}$ where $D_{k}$ is the diagonal part, $-E_{k}$
represents the strict lower triangular part and $E_{k}^{*}$ is the transpose
of $E_{k}$. Here we still use $\mathbb{T}_{k}$ as its corresponding discrete
matrix. If choosing $\mathbb{M}_{k}$ as the symmetric Gauss-Seidel
preconditioner involving $\mathbb{T}_{k}$, considering the positive
definiteness requirement of $M_{k}$ as in (3.2), we can choose [47] (chapter
4.1) (or [15])

$\mathbb{M}_{k}=\mathbb{T}_{k}+E_{k}^{*}D_{k}^{-1}E_{k}+\eta
I=D_{k}-E_{k}-E_{k}^{*}+E_{k}^{*}D_{k}^{-1}E_{k}+\eta I,$ (3.8)

where $\eta>0$ being a tiny positive constant and $I$ being the identity
matrix. The small perturbation with $\eta I$ is to guarantee the positive
definiteness of $M_{k}$. Actually, we have the proximal metric

$M_{k}=\mathbb{M}_{k}-\mathbb{T}_{k}=E_{k}^{*}D_{k}^{-1}E_{k}+\eta I\geq\eta
I,$ (3.9)

which is positive definite. However, we do not need to calculate the explicit
form of $\mathbb{M}_{k}$ or $M_{k}$. Instead, we can do it through rewriting
(3.7) as follows

$\displaystyle u^{k+1}$
$\displaystyle=(M_{k}+\mathbb{T}_{k})^{-1}[(M_{k}+\mathbb{T}_{k})u^{k}+\mathfrak{b}^{k}+\eta
u^{k}-(\mathbb{T}_{k}+\eta I)u^{k}]$ (3.10a)
$\displaystyle=u^{k}+\mathbb{M}_{k}^{-1}[\bar{\mathfrak{b}}^{k}-\bar{\mathbb{T}}_{k}u^{k}],$
(3.10b)

where $\bar{\mathfrak{b}}^{k}:=\mathfrak{b}^{k}+\eta u^{k}$ and
$\bar{\mathbb{T}}_{k}={\mathbb{T}}_{k}+\eta I$. This means the update (3.10)
is exactly the one time symmetric Gauss-Seidel iteration for the linear
equation $\bar{\mathbb{T}}_{k}u=\bar{\mathfrak{b}}^{k}$, which is also
equivalent to one time symmetric Gauss-Seidel iteration for the linear
equation ${\mathbb{T}}_{k}u={\mathfrak{b}}^{k}$ with preconditioner in (3.8).

Furthermore, it is proved that any finite preconditioned iterations still
provide a preconditioner [15] satisfying (3.2), i.e.,

$u^{k+(i+1)/n}=u^{k+i/n}+\mathbb{M}_{k}^{-1}(\bar{\mathfrak{b}}^{k}-\bar{\mathbb{T}}_{k}u^{k+i/n}),\quad
i=0,\ldots,n-1$ (3.11)

corresponds to
$u^{k+1}=u^{k}+\mathbb{M}_{k,n}^{-1}(\bar{\mathfrak{b}}^{k}-\bar{\mathbb{T}}_{k}u^{k})$
and we just need to choose large enough $\gamma_{+}$ according to fixed $n$
for meeting the requirement of $\mathbb{M}_{k,n}$ in (3.2). We thus built a
flexible framework for introducing the classical preconditioning techniques.
However, how to design efficient preconditioners for the corresponding
$\mathbb{T}_{k}$ especially the anisotropic cases is still very subtle and
challenging. We leave them to section 5.

Now, let’s turn to the discussion of the KL-properties of (2.4), (2.5), (2.9),
(2.12), and (2.14). The anisotropic cases (2.6), (2.7), (2.10), (2.13) are
completely similar and we omit them here. Let’s begin with the KL properties
of $L_{GY}(u,l)$.

######  Lemma 4.

The isotropic $L_{GY}(u,l)$ is KL-function.

###### Proof.

We see $\operatorname{gph}L_{GY}$ can be written as follows with
$l=(l_{1},l_{2})\in\mathbb{R}^{2MN}$. Denoting $x=(u,l_{1},l_{2},z)$ and
$|D(u)_{i}|:={|(Au-u_{0})_{i}|^{2}}/{2}$, we have

$\displaystyle\operatorname{gph}L_{GY}=\big{\\{}x:z=\sum_{i=1}^{M_{0}N_{0}}|D(u)_{i}|+\frac{\lambda}{2}\sum_{i=1}^{MN}\frac{|\nabla
u-l|_{i}^{2}}{2}+\sum_{i=1}^{MN}H(l_{i};\frac{\lambda}{\mu})\big{\\}}=$
$\displaystyle\big{\\{}x:z-(\sum_{i=1}^{M_{0}N_{0}}|D(u)_{i}|+\frac{\lambda}{2}\sum_{i=1}^{MN}\frac{|\nabla
u-l|_{i}^{2}}{2}+\sum_{i=1}^{M_{0}N_{0}}\sqrt{\frac{\lambda}{\mu}}|l_{i}|-\frac{|l_{i}|^{2}}{2})=0,l_{1i}^{2}+l_{2i}^{2}>\frac{\lambda}{\mu})\big{\\}}$
$\displaystyle\bigcup\big{\\{}x:z-(\sum_{i=1}^{M_{0}N_{0}}|D(u)_{i}|+\frac{\lambda}{2}\sum_{i=1}^{MN}\frac{|\nabla
u-l|_{i}^{2}}{2}+\frac{\lambda}{\mu}),0\leq
l_{1i}^{2}+l_{2i}^{2}\leq\frac{\lambda}{\mu})\big{\\}}.$

Denote $s:=\sqrt{l_{1}^{2}+l_{2}^{2}}$,
$U_{i}:=\\{(x,s):l_{1i}^{2}+l_{2i}^{2}=s_{i}^{2},s_{i}>\sqrt{\frac{\lambda}{\mu}}\\}$
and $V_{i}:=\\{(x,s):l_{1i}^{2}+l_{2i}^{2}=s_{i}^{2},0\leq
s_{i}\leq\sqrt{\frac{\lambda}{\mu}}\\}$ for $i=1,2,\cdots,M_{0}N_{0}$. We
found the above representation of $\operatorname{gph}L_{GY}$ can be formulated
as

$\displaystyle\big{(}\big{\\{}(x,s):z-(\sum_{i=1}^{M_{0}N_{0}}|D(u)_{i}|+\frac{\lambda}{2}\sum_{i=1}^{MN}\frac{|\nabla
u-l|_{i}^{2}}{2}+\sum_{i=1}^{M_{0}N_{0}}\sqrt{\frac{\lambda}{\mu}}s_{i}-\frac{s_{i}^{2}}{2})=0\big{\\}}\bigcap_{i=1}^{M_{0}N_{0}}U_{i}\big{)}$
$\displaystyle\bigcup\big{(}\big{\\{}(x,s):z-(\sum_{i=1}^{M_{0}N_{0}}|D(u)_{i}|+\frac{\lambda}{2}\sum_{i=1}^{MN}\frac{|\nabla
u-l|_{i}^{2}}{2}+\frac{\lambda}{\mu})=0\big{\\}}\bigcap_{i=1}^{M_{0}N_{0}}V_{i}\big{)}$

Since all the sets above are semialgebraic sets, with Tarski-Seidenberg
Theorem, i.e., Theorem 1, $\operatorname{gph}L_{GY}$ is semialgebraic sets.
$L_{GY}$ is semialgebraic function and is a KL-function. ∎

For the cases of the isotropic or anisotropic $L_{GR}(u,b)$, $L_{GM}(u,b)$,
and $L_{MS}(u,s)$, the proofs are quite similar to Lemma 5 and we omit here.

######  Lemma 5.

The functions $L_{GR}(u,b)$, $L_{GM}(u,b)$, and $L_{MS}(u,s)$ are semialgeraic
functions and thus are KL-functions.

For the KL property of $L_{HL}(u,b)$, we need to employ the o-minimal
structure. Actually, we have the following lemma.

######  Lemma 6.

The function $L_{HL}(u,b)$ is definable and is a KL-function.

###### Proof.

Denoting $b=e^{w}$, we have $\log(b)=w$. Since

$\operatorname{gph}L_{HL}(u,b)=\\{(u,b,z):z=\sum_{i=1}^{M_{0}N_{0}}\frac{|(Au-
u_{0})_{i}|^{2}}{2}+\sum_{i=1}^{MN}\frac{\mu}{2}(b_{i}\frac{|(\nabla
u)_{i}|}{\lambda}+b_{i}-\log(b_{i})-1)\\},$

it thus can be reformulated as

$\displaystyle\biggl{\\{}(u,b,z,w):z=\sum_{i=1}^{M_{0}N_{0}}\frac{|(Au-
u_{0})_{i}|^{2}}{2}+\sum_{i=1}^{MN}\frac{\mu}{2}(e^{w_{i}}\frac{|(\nabla
u)_{i}|}{\lambda}+e^{w_{i}}-w_{i}-1),b_{i}=e^{w_{i}}\biggl{\\}}$
$\displaystyle=\biggl{\\{}(u,b,z,w):\biggl{(}z-\sum_{i=1}^{M_{0}N_{0}}\frac{|(Au-
u_{0})_{i}|^{2}}{2}+\sum_{i=1}^{MN}\frac{\mu}{2}(e^{w_{i}}\frac{|(\nabla
u)_{i}|}{\lambda}+e^{w_{i}}-w_{i}-1)\biggl{)}^{2}$
$\displaystyle\qquad\qquad\qquad\qquad+\sum_{i=1}^{MN}(b_{i}-e^{w_{i}})^{2}=0\biggl{\\}}.$

The above set on $\\{(u,b,z,w)\\}$ is a zero set of
$P(u,b,z,w,e^{u},e^{b},e^{z},e^{w})=0$ with real polynomial function $P$. Here
$e^{u}=(e^{u_{1}},\cdots,e^{u_{MN}})$, $e^{b}=(e^{b_{1}},\cdots,e^{b_{MN}})$,
$e^{z}=(e^{z_{1}},\cdots,e^{z_{MN}})$, and $e^{w}=b$. We thus conclude the
$\operatorname{gph}L_{HL}(u,b)$ is definable due to Theorem 3. Thanks to
Theorem 2, $L_{HL}(u,b)$ is also definable and is a KL function. ∎

Furthermore, although the analysis of the KL exponent is very challenging,
however, for the KL exponent of the anisotropic $L_{GY}^{A}(u,l)$, we have the
following lemma.

######  Lemma 7.

The anisotropic $L_{GY}^{A}(u,l)$ is a KL-function with an exponent of
$\frac{1}{2}$.

###### Proof.

Suppose the vectorized $\boldsymbol{l}_{1}$ and $\boldsymbol{l}_{2}$ are
$(l_{1,1},\cdots,l_{1,MN})$ and $(l_{2,1},\cdots,l_{2,MN})$. Let’s introduce

$r_{i,j}(l_{i,j})=\begin{cases}\sqrt{\frac{\lambda}{\mu}}|l_{i,j}|-\frac{|l_{i,j}|^{2}}{2},&|l_{i,j}|\leq\sqrt{\frac{\lambda}{\mu}},\\\
\frac{\lambda}{2\mu},&|l_{i,j}|>\sqrt{\frac{\lambda}{\mu}},\end{cases}$

where $i=1,2$ and $j=1,\cdots,MN$. We thus rewrite $L_{GY}^{A}(u,l)$ as
follows

$\displaystyle\sum_{i=1}^{M_{0}N_{0}}\frac{|(Au-
u_{0})_{i}|^{2}}{2}+\sum_{i=1,2}\sum_{j=1}^{MN}\left(\frac{\mu}{2}{|(\nabla_{i}u)_{j}-l_{i,j}|^{2}}+\mu
r_{i,j}(l_{i,j})\right)$ (3.12)
$\displaystyle\sum_{i=1}^{M_{0}N_{0}}\frac{|(Au-
u_{0})_{i}|^{2}}{2}+\sum_{i=1,2}\sum_{j=1}^{MN}\left(\frac{\mu}{2}{|(\nabla_{i}u)_{j}-l_{i,j}|^{2}}+\mu\min_{k=1,2}\\{r_{i,j,k}(l_{i,j})+\zeta_{i,j,k}(l_{i,j})\\}\right)$

where
$r_{i,j,1}(l_{i,j})=\sqrt{\frac{\lambda}{\mu}}|l_{i,j}|-\frac{|l_{i,j}|^{2}}{2}$,
$r_{i,j,2}(l_{i,j})=\frac{\lambda}{2\mu}$, $i=1,2,j=1,\cdots,MN$, and

$\zeta_{i,j,1}(l_{i,j}):=I_{\\{l_{i,j}:|l_{i,j}|\leq\sqrt{\frac{\lambda}{\mu}}\\}}(l_{i,j}),\quad\zeta_{i,j,2}(l_{i,j}):=I_{\\{l_{i,j}:|l_{i,j}|>\sqrt{\frac{\lambda}{\mu}}\\}}(l_{i,j}).$
(3.13)

In order to put the piecewise polyhedral terms including $\zeta_{i,j,k}$ and
$|l_{i,j}|$ together, introducing

$\delta_{i,j,1}(l_{i,j}):=\zeta_{i,j,1}(l_{i,j})+\sqrt{{\lambda}{\mu}}|l_{i,j}|,\quad\delta_{i,j,2}(l_{i,j}):=\zeta_{i,j,2}(l_{i,j}),$

with (3.12), we can also reformulate $L_{GY}^{A}$ as follows

$L_{GY}^{A}(u,l)=\min_{k=1,2}\\{\mathfrak{F}_{i,j,k}(u,l)+\delta_{i,j,k}(l_{i,j})\\},\\\
$ (3.14)

where

$\displaystyle\mathfrak{F}_{i,j,1}(u,l):=\sum_{i=1}^{M_{0}N_{0}}\frac{|(Au-
u_{0})_{i}|^{2}}{2}+\frac{\mu}{2}\sum_{i=1,2}\sum_{j=1}^{MN}\left(\frac{|(\nabla_{i}u)_{j}-l_{i,j}|^{2}}{2}-\frac{|l_{i,j}|^{2}}{2}\right),$
(3.15)
$\displaystyle\mathfrak{F}_{i,j,2}(u,l):=\sum_{i=1}^{M_{0}N_{0}}\frac{|(Au-
u_{0})_{i}|^{2}}{2}+\frac{\mu}{2}\sum_{i=1,2}\sum_{j=1}^{MN}\left(\frac{|(\nabla_{i}u)_{j}-l_{i,j}|^{2}}{2}+\frac{\lambda}{2}\right).$
(3.16)

Denote $x:=(u,l)=(u,l_{1},l_{2})^{T}$. It can be seen that
$\delta_{i,j,1}(l_{i,j})$ and $\delta_{i,j,2}(l_{i,j})$ are polyhedral
functions, since both the constraint sets in (3.13) are polyhedral sets and
$|l_{i,j}|$ is also polyhedral function. Furthermore, it can be readily
checked that $\mathfrak{F}_{i,j,k}(u,l)$ can be written as follows,

$\mathfrak{F}_{i,j,k}(x)=\mathfrak{F}_{i,j,k}(u,l_{1},l_{2})=\frac{1}{2}x^{T}M_{k}x+x^{T}c_{k}+b_{k},\quad
k=1,2,$ (3.17)

where $M_{1}$ and $M_{2}$ are symmetric matrices,

$M_{1}=\Sigma^{T}\Sigma,\quad M_{2}=M_{1}-\text{Diag}[0,\mu I,\mu
I],\quad\Sigma:=\begin{bmatrix}A&0&0\\\ \sqrt{\mu}\nabla_{1}&-\sqrt{\mu}I&0\\\
\sqrt{\mu}\nabla_{2}&0&-\sqrt{\mu}I\end{bmatrix},$

and

$c_{k}=[-A^{T}u_{0},0,0]^{T},\quad k=1,2,\quad b_{1}=\|u_{0}\|_{2}^{2}/2,\quad
b_{2}=\|u_{0}\|_{2}^{2}/2+{\lambda}/{2}.$

Since $\mathfrak{F}_{i,j,k}$ can be written in (3.17) and $\delta_{i,j,k}$
with $k=1,2$ are polyhedral functions, by [38] (Corollary 5.2), we conclude
that $L_{GY}^{A}$ is a KL function of $(u,l_{1},l_{2})$ with KL exponent
$1/2$. ∎

##  4 Global Convergence of the Proposed Algorithms

In this section, we will discuss the convergence of (3.1) that can cover all
the models discussed in this paper. We will first prove the boundedness of the
iteration sequence $(u^{k},y^{k})$ in (3.1) for the corresponding models and
the convergence proof then follows. Actually, from the updates (3.1a) and
(3.1b), we can easily arrive at:

$\displaystyle L(u^{k},y^{k-1})+\frac{1}{2}\|u^{k}-u^{k-1}\|_{M_{k-1}}^{2}\leq
L(u^{k-1},y^{k-1}),$ (4.1a) $\displaystyle
L(u^{k},y^{k})+\frac{1}{2}\|y^{k}-y^{k-1}\|_{N_{k-1}}^{2}\leq
L(u^{k},y^{k-1}),$ (4.1b) $\displaystyle
L(u^{k},y^{k})+\frac{1}{2}\|u^{k}-u^{k-1}\|_{M_{k-1}}^{2}+\frac{1}{2}\|y^{k}-y^{k-1}\|_{N_{k-1}}^{2}\leq
L(u^{k-1},y^{k-1}),$ (4.1c)

which tells that $L(u^{k},y^{k})$ is bounded by $0$ and $L(u^{0},y^{0})$ and
is decreasing.

######  Lemma 8.

With $N_{k}=\mu I$, $M_{k}=M$ satisfying (3.2) and $A^{*}A$ has a bounded
inverse, the sequence $\\{(u^{k},l^{k}),k\in\mathbb{N}\\}$ of the iterations
(3.1) for $L_{GY}$ is bounded.

###### Proof.

By the condition on $A^{*}A$, we see $F(u)$ in (2.1) is coercive on $u$.
Since,

$F(u^{k})\leq\inf_{l\in Y}L_{GY}(u^{k},l)\leq L_{GY}(u^{k},l^{k})\leq
L_{GY}(u^{0},l^{0}),$

we conclude that $u^{k}$ must be bounded by the coercivity of $F(u)$. It can
be verified that $L_{GY}$ is also coercive on $l$ for any fixed $u$, we get
the boundedness of $l^{k}$. ∎

######  Lemma 9.

With $N_{k}=\lambda I/2$, $M_{k}$ satisfying (3.2) and $A^{*}A$ has a bound
inverse, the sequence $\\{(u^{k},b^{k}),k\in\mathbb{N}\\}$ of the iterations
(3.1) for $L_{GR}$ is bounded.

###### Proof.

The boundedness of $\\{b^{k}\\}$ follows from the constraint $I_{\\{0\leq
b\leq 1\\}}(b)$ and the projection

$b^{k+1}=\mathcal{P}_{[0,1]}(b^{k}-\frac{\lambda}{\mu}|\nabla
u^{k}|^{2}+1)\in[0,1],$ (4.2)

where $\mathcal{P}_{[0,1]}$ is the projection to $[0,1]$.

By the condition on $A^{*}A$, we get the the coercivity of $L_{GR}$ (2.1). We
thus get the boundedness of $u^{k}$ similarly to Lemma 8. ∎

Now, let’s turn to the boundedness of the iterative sequence of $L_{GM}$.

######  Lemma 10.

With $N_{k}=\mu I/2$, $0<b^{0}\leq 1$, $M_{k}$ satisfying (3.2) and $A^{*}A$
has a bounded inverse, the sequence $\\{(u^{k},b^{k}),k\in\mathbb{N}\\}$ of
the iterations (3.1) for $L_{GM}$ is bounded.

###### Proof.

We will first prove $0<b^{k+1}\leq 1$ by induction on assumption $0<b^{k}\leq
1$. The update for the $b_{k+1}$ is
$(\xi+1-\frac{1}{\sqrt{b_{k+1}}}+b_{k+1}-b_{k})=0$. With introducing
$x:=\sqrt{b_{k+1}}$, we arrive at

$x^{3}+px+q=0,\ \,\ p:=\xi+1-b_{k}\geq 0,\ q:=-1,\ \xi=\frac{|\nabla
u|^{2}}{\lambda}.$

Observing that $p>0$, $q<0$, we conclude
$\Delta:=\frac{q^{2}}{4}+\frac{p^{3}}{27}\geq\frac{1}{4}$ and
$-\frac{1}{2}+\sqrt{\Delta}\geq 0$. Using the celebrated Cardano’s formula for
the depressed cubic equation, we get

$x=\sqrt[3]{\frac{1}{2}+\sqrt{\Delta}}+\sqrt[3]{\frac{1}{2}-\sqrt{\Delta}}=\sqrt[3]{\frac{1}{2}+\sqrt{\Delta}}-\sqrt[3]{-\frac{1}{2}+\sqrt{\Delta}}.$
(4.3)

Besides, Noting that

$\frac{1}{2}+\sqrt{\Delta}=\sqrt[3]{-\frac{1}{2}+\sqrt{\Delta}}^{3}+\sqrt[3]{1}^{3}\leq(\sqrt[3]{-\frac{1}{2}+\sqrt{\Delta}}+\sqrt[3]{1})^{3}\rightarrow\sqrt[3]{\frac{1}{2}+\sqrt{\Delta}}\leq\sqrt[3]{-\frac{1}{2}+\sqrt{\Delta}}+{1}$

we obtain that $0<x\leq 1$, i.e. $0<b_{k+1}=x^{2}\leq 1$.

By the condition on $A^{*}A$, we get the the coercivity of $L_{GM}$ on $u$ in
(2.8). The boundedness of $\\{u^{k}\\}$ follows similarly to Lemma 8. ∎

######  Lemma 11.

With $N_{k}=\mu I/2$, $0\leq b^{0}\leq 1$, $M_{k}$ satisfying (3.2) and
$A^{*}A$ has a bounded inverse, the sequence
$\\{(u^{k},b^{k}),k\in\mathbb{N}\\}$ of the iterations (3.1) for $L_{HL}$ is
bounded.

###### Proof.

Note that $b_{k+1}=\operatorname*{arg\,min}_{b}\frac{\mu}{2}(\xi b+b-\log
b-1)+\frac{\|b-b_{k}\|_{\mu I}^{2}}{4}$ with $\xi:=\frac{|\nabla
u|^{2}}{\lambda}\geq 0$ leading to the quadratic equation
$b^{2}+(\xi+1-b^{k})b-1=0$. Using the quadratic formula and choosing the
positive root, we get

$b_{k+1}=(-a+\sqrt{a^{2}+4})/{2},\quad a=\xi+1-b_{k}\geq 0.$ (4.4)

By the assumption on $b^{0}$ and observing that $a^{2}+4\leq(a+2)^{2}$, we
have $0<b_{k+1}\leq 1$ by deduction. The sequence $\\{b^{k}\\}$ is thus
bounded.

By the condition on $A^{*}A$, we get the coercivity of $L_{HL}(u,b)$ for $u$
as in (2.11). The boundedness of $\\{u^{k}\\}$ follows similarly to Lemma 8. ∎

######  Lemma 12.

The sequence $\\{(u^{k},s^{k}),k\in\mathbb{N}\\}$ of the iterations (3.1) for
$L_{MS}$ is bounded.

###### Proof.

It can be seen that although $L_{MS}$ in (2.14) is not convex on $(u,s)$, it
is however coercive on $(u,s)$, since $L_{MS}(u,s)\rightarrow+\infty$ either
$|u|\rightarrow+\infty$ or $|s|\rightarrow+\infty$. Then $(u^{k},s^{k})$ must
be bounded by the boundednes of $L_{MS}(u^{k},s^{k})$ for all $k$. ∎

With (4.2), (4.3), and (4.4), we get the updates of $b^{k}$ for $L_{GR}$,
$L_{GM}$, and $L_{HL}$ correspondingly. Now, let’s discuss the update of
$l^{k}$ for $L_{GY}$.

######  Lemma 13.

Assuming $N_{k}=\kappa_{k}I$, for the $l_{i}$ updates of the anisotropic
$L_{GY}^{A}(u,l)$ with $i=1,2$ (or for the isotropic $L_{GY}(u,l)$ with $i=I$
case), we have

$l_{i}^{k+1}=\begin{cases}\hat{l}^{k}_{i}-\sqrt{a}\tau\frac{\hat{l}^{k}_{i}}{|\hat{l}^{k}_{i}|}\quad&\tau\sqrt{a}<|\hat{l}^{k}_{i}|<(1+\tau)\sqrt{a}\\\
0\quad&-\tau\sqrt{a}\leq|\hat{l}_{i}^{k}|\leq\tau\sqrt{a}\\\
\frac{1}{1+\tau}\hat{l}^{k}_{i}\quad&|\hat{l}^{k}_{i}|>(1+\tau)\sqrt{a}\end{cases},$
(4.5)

where
$a=\frac{\lambda}{\mu},\tau=\frac{\mu}{\kappa_{k}},\hat{l}^{k}_{i}=l_{i}^{k}+\frac{\mu}{\kappa_{k}}\nabla_{i}u^{k+1}$
with $i=1,2$ for the anisotropic case and the norm $|\cdot|$ denotes the
absolute value (or for the isotropic case with $i=I$,
$l^{k}_{I}=({l_{1}^{k}}_{I},{l_{2}^{k}}_{I})$, $\nabla_{I}=\nabla$, and
$|l_{I}|=|(l_{1},l_{2})|=\sqrt{l_{1}^{2}+l_{2}^{2}}$).

###### Proof.

Using (3.1b) and noting
$h_{\frac{\lambda}{\mu}}^{*}(l_{i})=H(l_{i};\frac{\lambda}{\mu})-{l_{i}^{2}}/{2}+\frac{1}{2}\frac{\lambda}{\mu}$,
we obtain that:

$\displaystyle l_{i}^{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{l_{i}}\frac{\mu}{2}\|\nabla_{i}u^{k+1}-l_{i}\|^{2}+\mu
H(l_{i};\frac{\lambda}{\mu})+\frac{1}{2}\|l_{i}-l_{i}^{k}\|_{\kappa_{k}I}^{2}$
$\displaystyle=\operatorname*{arg\,min}_{l_{i}}h_{a}^{*}(l_{i})+\frac{1}{2\mu}\|l_{i}-(l_{i}^{k}+\frac{\mu}{\kappa_{k}}\nabla_{i}u^{k+1})\|_{\kappa_{k}I}^{2}$
$\displaystyle=\operatorname*{arg\,min}_{l_{i}}h_{a}^{*}(l_{i})+\frac{1}{2\mu}\|l_{i}-\hat{l}^{k}_{i}\|_{\kappa_{k}I}^{2},$

where $h_{\tau}^{*}$ is the same as in (2.15). Using the first-order optimal
condition, we have

$0\in\partial
h_{a}^{*}(l_{i}^{k+1})+\frac{1}{\tau}(l_{i}^{k+1}-\hat{l}^{k}_{i})=\begin{cases}\sqrt{a}\frac{l_{i}^{k+1}}{|l_{i}^{k+1}|}+\frac{1}{\tau}(l^{k+1}_{i}-\hat{l}^{k}_{i})\quad&0<|l_{i}^{k+1}|<\sqrt{a}\\\
B_{\sqrt{a}}(0)-\frac{1}{\tau}\hat{l}^{k}_{i}\quad&l_{i}^{k+1}=0\\\
l^{k+1}_{i}+\frac{1}{\tau}(l^{k+1}_{i}-\hat{l}^{k}_{i})\quad&|l^{k+1}_{i}|\geq\sqrt{a}\end{cases},$
(4.6)

where $B_{\sqrt{a}}(0)$ is a closed ball with center 0 and radius $\sqrt{a}$
in $\mathbb{R}$ or $\mathbb{R}^{2}$. Solving this equation, we get (4.6). ∎

Now let’s turn to the global convergence of the iteration (3.1) for our
models. Actually, with the boundedness of $\\{(u^{k},y^{k})\\}$ proved above,
by [4] (Theorem 6.2), we can get the global convergence. The following
proposition will tell the convergence properties, i.e., Proposition 1, whose
proof can be obtained by slight modifications of the proof in [5] (Proposition
3.1) and is omitted here.

######  Proposition 1.

Let $\\{(u^{k},y^{k})\\}$ be a sequence generated by our algorithm (3.1).
Denoting $\Gamma$ as the set of accumulation points of the sequence
$\\{(u^{k},y^{k})\\}$, then the following statements hold:

  1. _(i)_

$\lim_{k\to\infty}(\|u^{k}-u^{k-1}\|_{M_{k-1}}^{2}+\|y^{k}-y^{k-1}\|_{N_{k-1}}^{2})=0$
and it follows from

$\sum_{k=1}^{\infty}(\|u^{k}-u^{k-1}\|_{M_{k-1}}^{2}+\|y^{k}-y^{k-1}\|_{N_{k-1}}^{2})<+\infty;$

  2. _(ii)_

if $(u^{k},y^{k})$ is bounded, then $\Gamma$ is a nonempty compact connected
set, and

  3. _(iii)_

$\Gamma\subset\emph{crit}\ L$. Moreover, $L(u,y)\equiv\zeta$ on $\Gamma$ is
constant, where $\zeta=\inf L(u^{k},y^{k})<+\infty$;

$\emph{dist}((u^{k},y^{k}),\Gamma)\to 0,\ \text{as}\ k\to\infty,$

where $\emph{crit}\ L:=\\{(\bar{x},\bar{y}):0\in\partial
L(\bar{x},\bar{y})\\}$.

The following Theorem 4 can be found in [5] (see Theorem 4.6 therein), which
tells that the iteration sequence $\\{(u^{k},y^{k})\\}$ has a finite length
and converges to a critical point.

######  Theorem 4.

Assuming $L$ in (3.1) is a KL-function and $\\{(u^{k},y^{k})\\}$ is bounded,
we have

$\sum_{k=1}^{+\infty}(\|u^{k}-u^{k-1}\|_{M_{k-1}}+\|y^{k}-y^{k-1}\|_{N_{k-1}})<+\infty,$
(4.7)

which means that $(u^{k},y^{k})$ converges to a critical point of $L$.

Actually, with the global convergence for convex and proximal alternating
optimization given in [7, 17], our preconditioned framework can also deal with
convex half quadratic models. We next consider the convergence rate of the
sequence $(u^{k},y^{k})$ for our models if the KL exponent of $L(u,y)$ is
known. This kind of convergence rate analysis is standard; see [3, 5, 38] for
more comprehensive analysis and is omitted here.

######  Theorem 5 (convergence rate).

Assume that $(u^{k},y^{k})$ converges to $(\bar{u},\bar{y})$ and $L(u,y)$ has
the KL property at $(\bar{u},\bar{y})$ with
$\psi(s)=cs^{1-\theta},\theta\in\left[0,1\right)$, $c>0$. Then the following
estimations hold:

  1. 1.

If $\theta=0$ then the sequence $(u^{k},y^{k})$ converges in a finite number
of steps.

  2. 2.

If $\theta\in\left(0,\frac{1}{2}\right]$ then there exist $c>0$ and
$\tau\in\left[0,1\right)$, such that $\|(u^{k},y^{k})-(\bar{u},\bar{y})\|\leq
c\tau^{k}$.

  3. 3.

If $\theta\in\left(\frac{1}{2},1\right]$ then there exists $c>0$, such that
$\|(u^{k},y^{k})-(\bar{u},\bar{y})\|\leq ck^{-\frac{1-\theta}{2\theta-1}}$.

With Lemma 7 and Theorem 5 we get the following convergence rate for (2.7).

######  Corollary 1.

Due to Lemma 7, the anisotropic $L_{GY}^{A}(u,l)$ is a KL function of $(u,l)$
with KL exponent $\frac{1}{2}$. The sequence $(u^{k},l^{k})$ generated by
preconditioned alternating iterations 3.1 thus has linear convergence rate due
to Theorem 5.

##  5 Numerical Part

In this part, we will first discuss the finite difference method and the
preconditioners. Although, the linear equations are in divergence form with
varying coefficients such as the linear subproblems for the $L_{GR}$,
$L_{GM}$, $L_{HL}$, and $L_{MS}$, however, we can still benefit from the
preconditioning framework by dealing with the linear subproblems with any
finite times symmetric Gauss-Seidel iterations. The numerical experiments
focused on image denoising and segmentation also show the efficiency of the
proconditioned framework.

###  5.1 The Preconditioners and the finite difference method

Now let’s turn to the discretization of the following equation with
$\mathbb{B}$ that is similar in (3.4) which covers all the linear equations
discussed in this paper

$\gamma(x)u+\nabla^{*}\mathbb{B}(x)\nabla
u=z(x),\quad\mathbb{B}(x)=\text{Diag}[d^{1}(x),d^{2}(x)],\ x\in\Omega.$ (5.1)

Here we use $d^{1}$, $d^{2}$ instead of $b_{1}$, $b_{2}$ as in (3.4) for the
convenience of lower indices to be employed. In order to discrete equation
(5.1), let’s first introduce the widely used forward and backward difference
operators $\nabla^{+}:=(\nabla_{x}^{+},\nabla_{y}^{+})^{T}$ and
$\nabla^{-}:=(\nabla_{x}^{-},\nabla_{y}^{-})^{T}$ [17]

$\nabla_{x}^{+}u_{i,j}=\begin{cases}u_{i+1,j}-u_{i,j}\quad&1\leq i<m,\\\
0\quad&i=m,\end{cases}\quad\nabla_{y}^{+}u_{i,j}=\begin{cases}u_{i,j+1}-u_{i,j}\quad&1\leq
j<n,\\\ 0\quad&i=n,\end{cases}$ (5.2)

and

$\nabla_{x}^{-}u_{i,j}=\begin{cases}u_{1,j}\quad&i=1,\\\
u_{i,j}-u_{i-1,j}\quad&1<i\leq m,\\\
-u_{m-1,j}\quad&i=m,\end{cases}\quad\nabla_{y}^{-}u_{i,j}=\begin{cases}u_{i,1}\quad&j=1,\\\
u_{i,j}-u_{i,j-1}\quad&1<j\leq n,\\\ -u_{i,n-1}\quad&j=n.\end{cases}$ (5.3)

If choosing $\nabla^{+}$ as the gradient and its adjoint as the divergence
operator, we get

$\gamma_{i,j}u_{i,j}-\left(\nabla_{x}^{-}(d^{1}_{i,j}\nabla_{x}^{+}u_{i,j})+\nabla_{y}^{-}(d^{2}_{i,j}\nabla_{y}^{+}u_{i,j})\right)=z_{i,j}.$
(5.4)

Written (5.4) in detail with discrete $\nabla^{+}$ and $\nabla^{-}$, taking
the interior points for example, we have

$\displaystyle\big{(}\gamma_{i,j}+(d_{i,j}^{1}+d_{i-1,j}^{1}+d_{i,j}^{2}+d_{i,j-1}^{2})\big{)}u_{i,j}-\big{(}d_{i,j}^{1}u_{i+1,j}+d_{i-1,j}^{1}u_{i-1,j}+d_{i,j}^{2}u_{i,j+1}$
$\displaystyle+d_{i,j-1}^{2}u_{i,j-1}\big{)}=z_{i,j}.$

For the boundary and corner points, we can also get the corresponding updates
and see Tables 1 and 2 for more details. We call this discretization scheme as
the normal formula of finite difference (shorten as NFFD).

There is another symmetric finite difference approximation of (5.1) [51] (see
section 3.4 of [51]) with the following additional
$\widetilde{\nabla}_{x}^{-},\widetilde{\nabla}_{x}^{+},\widetilde{\nabla}_{y}^{-},\widetilde{\nabla}_{y}^{+}$
defined by

$\widetilde{\nabla}_{x}^{-}u_{i,j}=\begin{cases}0\quad&i=1,\\\
u_{i,j}-u_{i-1,j}\quad&1<i\leq
m,\end{cases}\quad\widetilde{\nabla}_{y}^{-}u_{i,j}=\begin{cases}0\quad&j=1,\\\
u_{i,j}-u_{i,j-1}\quad&1<j\leq n,\end{cases}$ (5.5)

and

$\widetilde{\nabla}_{x}^{+}u_{i,j}=\begin{cases}u_{2,j}\quad&i=1,\\\
u_{i+1,j}-u_{i,j}\quad&1<i\leq m,\\\
-u_{m,j}\quad&i=m,\end{cases}\quad\widetilde{\nabla}_{y}^{+}u_{i,j}=\begin{cases}u_{i,2}\quad&j=1,\\\
u_{i,j+1}-u_{i,j}\quad&1<j\leq n,\\\ -u_{i,n}\quad&j=n.\end{cases}$ (5.6)

The symmetric finite difference approximation of (5.1) [51] reads as follows
(see section 3.4 of [51])

$\displaystyle\gamma_{i,j}u_{i,j}-\frac{1}{2}\big{(}\nabla_{x}^{-}(d^{1}_{i,j}\nabla_{x}^{+}u_{i,j})+\widetilde{\nabla}_{x}^{+}(d^{1}_{i,j}\widetilde{\nabla}_{x}^{-}u_{i,j})+\nabla_{y}^{-}(d^{2}_{i,j}\nabla_{y}^{+}u_{i,j})$
(5.7)
$\displaystyle+\widetilde{\nabla}_{y}^{+}(d^{2}_{i,j}\widetilde{\nabla}_{y}^{-}u_{i,j})\big{)}=z_{i,j},$

Written out in full and for the interior points, we obtain

$\displaystyle\big{(}\gamma_{i,j}+(\alpha_{i,j}^{1}+\alpha_{i-1,j}^{1}+\alpha_{i,j}^{2}+\alpha_{i,j-1}^{2})\big{)}u_{i,j}-\big{(}\alpha_{i,j}^{1}u_{i+1,j}+\alpha_{i-1,j}^{1}u_{i-1,j}$
(5.8)
$\displaystyle+\alpha_{i,j}^{2}u_{i,j+1}+\alpha_{i,j-1}^{2}u_{i,j-1}\big{)}=z_{i,j},$

where $\alpha_{i,j}^{1}=\frac{1}{2}(d_{i,j}^{1}+d_{i+1,j}^{1})$,
$\alpha_{i,j}^{2}=\frac{1}{2}(d_{i,j}^{2}+d_{i,j+1}^{2})$ and see the Table 3
for more details. We call this discretization scheme as the symmetric formula
of finite difference (shorten as SFFD). Here we mainly focus on the NFFD and
SFFD for these nonhomogeneous elliptic equations in divergence form.

With the stencils in Tables 2, 3 and the corresponding locations in Table 1,
one can do the classical symmetric Gauss-Seidel including the symmetric red-
black Gauss-Seidel iteration (shorten as SRBGS henceforth) conveniently [14,
47]. To do SRBGS iteratons, one should mark all points as red or black points
first. For example, let $\Omega_{red}:=\\{(i,j)\in\Omega|i+j\ \text{even}\\}$
and $\Omega_{black}:=\\{(i,j)\in\Omega|i+j\ \text{odd}\\}$. Then one can first
update all the red points with the five-point stencils in Table 2 or 3 using
all the black neighbor points followed by updating all the black points with
the five-point stencils in Table 2 or 3 with the updated red neighbor points.
Finally, update all the red points again. Then the one cycle Gauss-Seidel
iterations with updating order “red $\rightarrow$ black $\rightarrow$ red” is
finished and it is called one time SRBGS iteration. For more details of the
SRBGS, we refer to [15] (section 4.1.3). Henceforth, we will use the following
notation to denote the $n$ times SRBGS iterations for the equation (5.1) with
initial iteration value $u^{0}$

$\text{SRBGS}(\gamma(x),d^{1},d^{2},z,u^{0},n),$

which can be seen as one n-folds SRBGS preconditioned iteration according to
(3.11) for (5.1). Due to the positive definiteness requirement of $M_{k}$,
taking the $k$-th nonlinear preconditioned alternating minimization of the
anisotropic $L_{HL}$ model (2.13) for example, we actually do the following
preconditioned iteration with tiny positive constant $\eta$ and initial
iterative value $u^{k}$

$u^{k+1}:=\text{SRBGS}(1.0+\eta,b_{1}^{k},b_{2}^{k},z^{k}+\eta
u^{k},u^{k},n),$ (5.9)

according to (3.10) in Remark 1.

Actually, the discretization schemes (5.4) and (5.7) are both consistent with
the Neumann boundary condition as in (5.1) by the following Remark 2.
Furthermore, it can be readily checked that while $b_{1}$ and $b_{2}$ are both
constants, the schemes (5.4) and (5.7) are equivalent; see the following
Remark 3.

######  Remark 2.

We point out that (5.4) and (5.7) are both with the Neumann boundary condition
$\frac{\partial u}{\partial\nu}=0$, where $\nu$ is the outward pointing unit
normal vector of $\Omega$. For example, for the $\mathcal{S}_{S}$ points in
Table 1, we have $u_{m+1,j}=u_{m,j},1<j<n$. Moreover, with the notation
$\Lambda_{m.j}:=\alpha_{m-1,j}^{1}+\alpha_{m,j}^{2}+\alpha_{m,j-1}^{2}$ for
convenience, we have

$(\gamma_{i,j}+\Lambda_{m,j})u_{m,j}-(\alpha_{m-1,j}^{1}u_{m-1,j}+\alpha_{m,j}^{2}u_{m,j+1}+\alpha_{m,j-1}^{2}u_{m,j-1})=z_{m,j}$
(5.10)

in (5.8) and for the left item in (5.7), we have

$\displaystyle\gamma_{i,j}u_{m,j}-\frac{1}{2}[-d_{m-1,j}^{1}\nabla_{x}^{+}u_{m-1,j}-d^{1}_{m,j}\widetilde{\nabla}_{x}^{-}u_{m,j}+(d^{2}_{m,j}\nabla^{+}_{y}u_{m,j}-d^{2}_{m,j-1}\nabla^{+}_{y}u_{m,j-1})$
$\displaystyle+(d^{2}_{m,j+1}\widetilde{\nabla}^{-}_{y}u_{m,j+1}-d^{2}_{m,j}\widetilde{\nabla}^{-}_{y}u_{m,j})]$
$\displaystyle=\gamma_{i,j}u_{m,j}-\frac{1}{2}[-(d^{1}_{m-1,j}(u_{m,j}-u_{m-1,j}))+(-d^{1}_{m,j}(u_{m,j}-u_{m-1,j}))+(d^{2}_{m,j}(u_{m,j+1}$
$\displaystyle-
u_{m,j})-d^{2}_{m,j-1}(u_{m,j}-u_{m,j-1}))+(d^{2}_{m,j+1}(u_{m,j+1}-u_{m,j})+d^{2}_{m,j}(u_{m,j}-u_{m,j-1}))]$
$\displaystyle=(\gamma_{i,j}+\Lambda_{m,j})u_{m,j}-(\alpha_{m-1,j}^{1}u_{m-1,j}+\alpha_{m,j}^{2}u_{m,j+1}+\alpha_{m,j-1}^{2}u_{m,j-1}),$

which is equal to the left-hand side of (5.10).

######  Remark 3.

(5.4) and (5.7) are the same when $d^{1}$ and $d^{2}$ both are constants.

$i$$j$$\scriptstyle\mathcal{S}_{SW}$$\scriptstyle\mathcal{S}_{S}$$\scriptstyle\mathcal{S}_{SE}$$\scriptstyle\mathcal{S}_{W}$$\scriptstyle\mathcal{S}_{O}$$\scriptstyle\mathcal{S}_{E}$$\scriptstyle\mathcal{S}_{NW}$$\scriptstyle\mathcal{S}_{N}$$\scriptstyle\mathcal{S}_{NE}$

Table 1: Finite-difference stencils for (5.1) with homogeneous Neumann boundary conditions. The highlighted entry denotes the center element and $u$ is assumed to be extended by arbitrary values outside $\Omega$. $\underbrace{\begin{bmatrix}\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{1,1}$}&d_{1,1}^{2}\\\ -d_{1,1}^{1}&\end{bmatrix}}_{=\mathcal{S}_{NW}}$ | $\underbrace{\begin{bmatrix}-d_{1,j-1}^{2}&\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{1,j}$}&-d_{1,j}^{2}\\\ &-d_{1,j}^{1}&\end{bmatrix}}_{=\mathcal{S}_{N}}$ | $\underbrace{\begin{bmatrix}-d_{1,n-1}^{2}&\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{1,n}$}\\\ &-d_{1,n}^{1}\end{bmatrix}}_{=\mathcal{S}_{NE}}$
---|---|---
$\underbrace{\begin{bmatrix}-d_{i-1,1}^{1}&\\\ \hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\mu\Sigma_{i,1}$}&-d_{i,1}^{2}\\\ -d_{i,1}^{1}&\end{bmatrix}}_{=\mathcal{S}_{W}}$ | $\underbrace{\begin{bmatrix}&-d_{i-1,j}^{1}&\\\ -d_{i,j-1}^{2}&\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{i,j}$}&-d_{i,j}^{2}\\\ &-d_{i,j}^{1}&\end{bmatrix}}_{=\mathcal{S}_{O}}$ | $\underbrace{\begin{bmatrix}&-d_{i-1,n}^{1}&\\\ -d_{i,n-1}^{2}&\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{i,n}$}&\\\ &-d_{i,n}^{1}&\end{bmatrix}}_{=\mathcal{S}_{E}}$
$\underbrace{\begin{bmatrix}&-d_{m-1,1}^{1}&\\\ &\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{m,1}$}&-d_{m,1}^{2}\end{bmatrix}}_{=\mathcal{S}_{SW}}$ | $\underbrace{\begin{bmatrix}&-d_{m-1,j}^{1}&\\\ -d_{m,j-1}^{2}&\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{m,j}$}&-d_{m,j}^{2}\end{bmatrix}}_{=\mathcal{S}_{S}}$ | $\underbrace{\begin{bmatrix}&-d_{m-1,n}^{1}&\\\ -d_{m,n-1}^{2}&\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{m,n}$}&\end{bmatrix}}_{=\mathcal{S}_{SE}}$
Table 2: The five-point stencils for NFFD discretization scheme (5.4). $\Sigma_{i,j}:=d_{i,j-1}^{2}+d_{i,j}^{2}+d_{i-1,j}^{1}+d_{i,j}^{1}$ for the interior points. On the corner, $\Sigma_{1,1}:=d_{1,1}^{2}+d_{1,1}^{1}$ and other points on the corner are similar. For the interior boundary points, $\Sigma_{i,1}:=d_{i-1,1}^{1}+d_{i,1}^{1}+d_{i,1}^{2}$ and other points on the interior boundary are similar.  $\underbrace{\begin{bmatrix}\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{1,1}$}&-\alpha_{1,1}^{2}\\\ -\alpha_{1,1}^{1}&\end{bmatrix}}_{=\mathcal{S}_{NW}}$ | $\underbrace{\begin{bmatrix}-\alpha_{1,j-1}^{2}&\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{1,j}$}&-\alpha_{1,j}^{2}\\\ &-\alpha_{1,j}^{1}&\end{bmatrix}}_{=\mathcal{S}_{N}}$ | $\underbrace{\begin{bmatrix}-\alpha_{1,n-1}^{2}&\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{1,n}$}\\\ &-\alpha_{1,n}^{1}\end{bmatrix}}_{=\mathcal{S}_{NE}}$
---|---|---
$\underbrace{\begin{bmatrix}-\alpha_{i-1,1}^{1}&\\\ \hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{i,1}$}&-\alpha_{i,1}^{2}\\\ -\alpha_{i,1}^{1}&\end{bmatrix}}_{=\mathcal{S}_{W}}$ | $\underbrace{\begin{bmatrix}&-\alpha_{i-1,j}^{1}&\\\ -\alpha_{i,j-1}^{2}&\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{i,j}$}&-\alpha_{i,j}^{2}\\\ &-\alpha_{i,j}^{1}&\end{bmatrix}}_{=\mathcal{S}_{O}}$ | $\underbrace{\begin{bmatrix}&-\alpha_{i-1,n}^{1}&\\\ -\alpha_{i,n-1}^{2}&\hbox{\pagecolor{white!90!black}$\lambda+\Sigma_{i,n}$}&\\\ &-\alpha_{i,n}^{1}&\end{bmatrix}}_{=\mathcal{S}_{E}}$
$\underbrace{\begin{bmatrix}&-\alpha_{m-1,1}^{1}&\\\ &\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{m,1}$}&-\alpha_{m,1}^{2}\end{bmatrix}}_{=\mathcal{S}_{SW}}$ | $\underbrace{\begin{bmatrix}&-\alpha_{m-1,j}^{1}&\\\ -\alpha_{m,j-1}^{2}&\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{m,j}$}&-\alpha_{m,j}^{2}\end{bmatrix}}_{=\mathcal{S}_{S}}$ | $\underbrace{\begin{bmatrix}&-\alpha_{m-1,n}^{1}&\\\ -\alpha_{m,n-1}^{2}&\hbox{\pagecolor{white!90!black}$\gamma_{i,j}+\Sigma_{m,n}$}&\end{bmatrix}}_{=\mathcal{S}_{SE}}$
Table 3: The five-point stencils for SFFD scheme (5.7). Here
$\alpha_{i,j}^{1}=\frac{1}{2}(d_{i,j}^{1}+d_{i+1,j}^{1})$ with $1\leq
i<m,1\leq j\leq n$, $\alpha_{i,j}^{2}=\frac{1}{2}(d_{i,j}^{2}+d_{i,j+1}^{2})$
with $1\leq i\leq m,1\leq j<n$ and
$\Sigma_{i,j}:=\alpha_{i,j-1}^{2}+\alpha_{i,j}^{2}+\alpha_{i-1,j}^{1}+\alpha_{i,j}^{1}$
on the interior points. On the corner,
$\Sigma_{1,1}:=\alpha_{1,1}^{2}+\alpha_{1,1}^{1}$ and other points on the
corner are similar. On the interior boundary,
$\Sigma_{i,1}:=\alpha_{i-1,1}^{1}+\alpha_{i,1}^{1}+\alpha_{i,1}^{2}$ and other
points on the interior boundary are similar.  Preconditioned alternating
minimization for $L_{GY}$ model
---
Initialization: |  $(u^{0},l^{0})$, $\lambda>0$, $\mu>0$, $0<\bar{\eta}\leq\eta^{+}$, $N_{k}\equiv\kappa I$
| with $0<\kappa\leq\kappa^{+}$,
|  $n\geq 1$ inner iterations for SRBGS
Iteration: | $z^{k}$ |  $=u_{0}-\mu\operatorname{div}l^{k}$
| $u^{k+1}$ |  $=\text{SRBGS}(1.0+\eta,\mu,\mu,z^{k}+\eta u^{k},u^{k},n)$
| $l^{k+1}$ |  = (4.5) for the isotropic or anisotropic case
Table 4: Preconditioned alternating minimization for the isotropic or
anisotropic Geman-Yang model $L_{GY}$ (2.5) or $L_{GY}^{A}$ (2.7). The SRBGS
iterations are used as preconditioner for dealing with the equations of
constant coefficient for updating $u^{k+1}$.  Preconditioned alternating
minimization for $L_{GM}$, $L_{GR}$ or $L_{HL}$ model
---
Initialization: |  $(u^{0},b^{0})$ with $0<b^{0}\leq 1$, $\lambda>0$, $\mu>0$, $0<\eta\leq\eta^{+}$,
|  $N_{k}\equiv\lambda I/2$ for $L_{GR}$, $N_{k}\equiv\mu I/2$ for $L_{GM}$
and $L_{HL}$,
|  $n\geq 1$ inner iterations for SRBGS,
Iteration: | $u^{k+1}$ |  $=\begin{cases}\text{for $L_{GR}$ case:}\\\ \text{SRBGS}(1.0+\eta,\mu d^{1k},\mu d^{2k},u_{0}+\eta u^{k},u^{k},n)\quad\\\ \text{for $L_{GM}$ or $L_{HL}$ case:}\\\ \text{SRBGS}(1.0+\eta,\frac{\lambda}{\mu}d^{2k},\frac{\lambda}{\mu}d^{1k},u_{0}+\eta u^{k},u^{k},n)\quad\end{cases}$
| $b^{k+1}$ |  $=\begin{cases}\eqref{eq:proj:p}\ \text{for $L_{GR}$ case}\\\ \eqref{eq:update:b:gm}\ \text{for $L_{GM}$ case}\\\ \eqref{eq:update:b:hl}\ \text{for $L_{HL}$ case}\end{cases}$
Table 5: Preconditioned alternating minimization for the isotropic or
anisotropic Geman-McClure $L_{GM}$ models (2.9) or (2.10), the Geman-Reynolds
models $L_{GR}$ (2.4) and (2.6) and the Hebert-Leahy models $L_{HL}$ (2.12)
and (2.13). The SRBGS iterations are used as preconditioner for dealing with
equations of varying coefficients for updating $u^{k+1}$. Preconditioned
alternating minimization for $L_{MS}$ model
---
Initialization: |  $(u^{0},s^{0})$, $\lambda>0$, $\mu>0$, $\epsilon>0$, $0<\eta\leq\eta^{+}$, $0<\gamma\leq\gamma^{+}$,
|  $n\geq 1$ inner iterations for symmetric Gauss-Seidel
Iteration: | $u^{k+1}$ |  $=\text{SRBGS}(1.0+\eta,2\alpha(s^{k})^{2},2\alpha(s^{k})^{2},u_{0}+\eta u^{k},u^{k},n)$
| $s^{k+1}$ |  $=\text{SRBGS}(2\alpha|\nabla u^{k}|^{2}+\frac{\lambda}{2\epsilon}+\gamma,2\lambda\epsilon,2\lambda\epsilon,\frac{\lambda}{2\epsilon}+\gamma s^{k},s^{k},n)$
Table 6: Preconditioned alternating minimization the Ambrosio–Tortorelli
approximated Mumford-Shah model $L_{MS}$ (2.14). The preconditioned iterations
are used both for dealing with equations of varying coefficients for updating
$u^{k+1}$ and $s^{k+1}$.

###  5.2 Numerical tests: image denoising and segmentation

In this section, we present the detailed performance of these models. All
experiments are performed in Matlab 2019a on a 64-bit PC with an Inter(R)
Core(TM) i5-9300HQ CPU (2.40Hz) and 12 GB of RAM. It can be seen from Table 7
that the anisotropic $L_{GM}$, $L_{GR}$, $L_{HL}$, and $L_{GY}$ models are
competitive compared to the widely used TV (total variation) model.
Especially, the $L_{HL}$ model with the SFFD has some improvement compared
with TV. Furthermore, for most of the images tested, the SFFD performs better
than the NFFD. Generally, with SFFD, one can get better PSNR according to
Table 7.

Figure 1 shows the effectiveness for image denoising and segmentation of the
models discussed. It can be seen that (i.e., images (a)-(l) of Figure 1) all
these nonconvex regularizations including $L_{GM}$, $L_{GR}$, $L_{HL}$, and
$L_{GY}$ do not have stair-casing effect as TV regularization, no matter for
the isotropic or anisotropic case. The segmentation of the well-known $L_{MS}$
model is very appealing with our preconditioning technique; see images (m)-(p)
of Figure 1.

Moreover, Figure 2 shows our motivation for introducing the preconditioning.
Let’s take the $L_{HL}$ model for example. Other models including $L_{GM}$,
$L_{GR}$, $L_{MS}$, and $L_{GY}$ are similar to our observations and we omit
the comparisons for compactness. First, it can be seen that one can benefit by
getting more lower energy and better PSNR from the proximal terms in (3.1) by
solving with CG (conjugate gradient) directly compared to solving the original
systems (1.2) without any proximal terms. Furthermore, with the efficient
SRBGS preconditioners, the nonlinear iterations converge much faster while
obtaining the high PSNR and low energy more quickly compared to solving the
linear system with proximal terms by CG. Although CG solver with low accuracy
is also fast, however, there is no convergence guarantee for the whole
nonlinear iterations with CG solver. The proposed preconditioned framework can
obtain lower energy and better PSNR faster with a convergence guarantee, which
is very promising.

Furthermore, experimentally, we observed that all the energy functions of
$L_{GM}$, $L_{GR}$, $L_{HL}$, $L_{GY}$, and $L_{GM}$ are decreasing
monotonously, which are consistent with the proximal or preconditioned
alternating minimization framework (3.1). We select some representatives as in
image (b) of Figure 3 for compactness.

Moreover, for the local linear convergence rate of $L_{GY}^{A}$ whose KL
exponent can be proved to be $1/2$ in Lemma 7, the numerical test also shows
the asymptotic linear convergence rate as in image (a) of Figure 3, which is
also shown theoretically in Corollary 1.

Table 7: Comparison for the anisotropic image denoising models with PSNR (peak signal to noise ratio). The noisy images are as follows: Lena1, Lena2 with size 512 $\times$512 and Monarch1, Monarch2 with size 768 $\times$ 512\. The usual zero mean Gaussian white noise of variance $\sigma=0.1$ for Lena1 or Monarch1 and $\sigma=0.05$ for Lena2 and Monarch2. The parameters for the corresponding models are as follows. For $L_{GM}$ model, we choose $\mu=0.02$, $\lambda=0.05$ for $\sigma=0.1$ case and $\mu=0.007$, $\lambda=0.004$ for $\sigma=0.05$ case. For $L_{GR}$ and $L_{GY}$ models, we choose $\mu=3$, $\lambda=0.01$ for $\sigma=0.1$ case and $\mu=1.5$, $\lambda=0.05$ for $\sigma=0.05$ case. For $L_{HL}$ model, we choose $\mu=0.005$, $\lambda=0.001$ for $\sigma=0.1$ case and $\mu=0.002$, $\lambda=0.0005$ for $\sigma=0.05$ case. For the anisotropic TV model, the regularization parameter $\alpha$ is chosen as the variance of the noise, i.e., $\alpha=\sigma$. The models along with the finite difference schemes that obtained the best PSNR are highlighted in bold face.  | $L_{GM}$ | $L_{GR}$ | $L_{HL}$ | TV | $L_{GY}$
---|---|---|---|---|---
NFFD | SFFD | NFFD | SFFD | NFFD | SFFD | * | *
Lena1 | 29.17 | 29.57 | 28.41 | 29.23 | 29.23 | $\boldsymbol{30.02}$ | 29.22 | 29.31
Monarch1 | 28.93 | 29.75 | 28.48 | 29.49 | 29.31 | $\boldsymbol{30.29}$ | 29.14 | 29.62
Lena2 | 32.33 | $\boldsymbol{32.94}$ | 32.26 | 32.47 | 32.47 | 32.90 | 31.85 | 32.36
Monarch2 | 32.90 | 33.71 | 32.91 | 33.32 | 33.25 | $\boldsymbol{33.94}$ | 32.61 | 33.18

(a) Monarch-noise01

(b) Monarch-noise005

(c) Lena-noise01

(d) Lena-noise005

(e) Monarch-ani-GM-01

(f) Monarch-ani-GR-005

(g) Lena-ani-HL-01

(h) Lena-ani-GY-005

(i) Monarch-iso-GM-01

(j) Monarch-iso-GR-005

(k) Lena-iso-HL-01

(l) Lena-iso-GY-005

(m) Man

(n) Man-segmentation

(o) Tulips

(p) Tulips-segmentation

Figure 1: Image (a) shows the noisy Monarch image of size $512\times 512$
corrupted by Gaussian noise of variance $\sigma=0.1$. Image (e) and (i) show
the denoised images of (a) by the anisotropic $L_{GM}$ model with parameters
$\mu=0.02$, $\lambda=0.05$ and the isotropic $L_{GM}$ model with parameters
$\mu=0.02$, $\lambda=0.001$ correspondingly. Image (b) shows the Monarch noisy
image of size $512\times 512$ corrupted by Gaussian noise of variance
$\sigma=0.0.5$. Image (f) and (j) show the denoised images of (b) by the
anisotropic $L_{GR}$ model with parameters $\mu=1.5$, $\lambda=0.05$ and the
isotropic $L_{GR}$ model with parameters $\mu=1.5$, $\lambda=0.05$
correspondingly. Image (c) shows the corresponding noisy Lena image of size
$512\times 512$ corrupted by Gaussian noise of variance $\sigma=0.1$. Image
(g) and (k) show the denoised images of image (c) by the anisotropic $L_{HL}$
model with parameters $\mu=0.005$, $\lambda=0.001$ and the isotropic $L_{HL}$
model with parameters $\mu=0.005$, $\lambda=0.0005$ correspondingly. Image (d)
shows the noisy Lena image of size $512\times 512$ corrupted by Gaussian noise
of variance $\sigma=0.05$. Image (h) and (l) show the denoised images of (a)
by the anisotropic $L_{GY}$ model with parameters $\mu=1.5$, $\lambda=0.05$
and the isotropic $L_{GY}$ model with parameters $\mu=1.5$, $\lambda=0.005$
correspondingly. Image (m) shows the original 512$\times$512 Man image and
image (n) shows the segmented man image by the MS model $L_{MS}$ with
parameters $\alpha=5000$, $\lambda=0.1$ and $\epsilon=0.02$. Image (o) shows
the original 512$\times$512 tulips image. Image (p) shows the segmented tulips
image by the MS model $L_{MS}$ with parameters $\alpha=3000$, $\lambda=0.1$
and $\epsilon=0.02$.

(a)  Energy (b) PSNR

Figure 2: Figures (a) or (b) shows the comparison with computational time
between preconditioned alternating minimization with 10 times inner symmetric
Red-Black Gauss-Seidel (SRBGS) iterations and the conjugate gradient (CG)
solver. SFFD-pre denotes the preconditioned alternating minimization by 10
times inner SRBGS preconditioned iterations with SFFD scheme. CG-pro and CG-
pro2 denote solving (3.1a) and (3.1b) by CG iterations with $10^{-3}$ and
$10^{-6}$ error stopping criterion correspondingly, while CG-Npro denotes
solving the linear subproblems by CG iterations of the anisotropic $L_{HL}$
model of $10^{-3}$ error stopping criteron without proximal terms in (3.1a)
and (3.1b). The computation is based on the Lena image of size 512$\times$512
corrupted by Gaussian noise of variance $\sigma=0.1$ for the anistropic
$L_{HL}$ model with parameters$\mu=0.005,\lambda=0.05$.

(a) The local linear convergence rate for the anisotropic GY model
$L_{GY}^{A}$ (b) Energy descent of some models

Figure 3: Figure (a) shows the local linear convergence rate. The computation
is by the anistropic GY model $L_{GY}^{A}$ for the lena image of size 512
$\times$512 with parameters $\mu=1.5$, $\lambda=0.005$ for the Gaussian
variance $\sigma=0.05$ cases. Figure (b) shows the energy descent. The
computation is by different models for the Monarch image of size $768\times
512$. All models except MS model are used for denoising the noisy Monarch
image corrupted by Gaussian noise of variance $\sigma=0.1$ and the MS model
$L_{MS}$ is used to segmentation. The parameters of this models are as
follows: $\mu=0.02,\lambda=0.001$ for the GM-SFFD-iso case (the isotropic
$L_{GM}$ model with SFFD), $\mu=1.5,\lambda=0.05$ for the GR-NFFD-ani case
(the anisotropic $L_{GR}$ model with NFFD) and GY-ani cases,
$\mu=0.005,\lambda=0.001$ for the HL-NFFD-ani case (the anisotropic $L_{HL}$
model with SFFD) and $\alpha=5000,\lambda=0.1,\epsilon=0.02$ for the MS
$L_{MS}$ case.

##  6 Discussion and Conclusions

In this paper, we investigated the proposed framework of preconditioned
alternating minimization methods for some typical nonconvex and nonlinear
models. With specially designed proximal terms, we can reformulate solving the
linear subproblems as the classical preconditioned iterations. We thus can
avoid solving the linear subproblems exactly or high accurately. Meanwhile, we
can also get the global convergence guarantee. Moreover, we can also obtain
high-quality reconstructions including image denoising and segmentation more
efficiently as shown in numerics. For the future study, we will consider more
general cases including the constraints which are also very important.

## References

  * [1] M. Allain, J. Idier, and Y. Goussard.  On global and local convergence of half-quadratic algorithms.  IEEE T. Image Process., 15(5):1130–1142, 2006.
  * [2] Luigi Ambrosio and Vincenzo Maria Tortorelli.  Approximation of functional depending on jumps by elliptic functional via t-convergence.  Comm. Pure Appl. Math., 43(8):999–1036, 1990.
  * [3] Hedy Attouch and Jérôme Bolte.  On the convergence of the proximal algorithm for nonsmooth functions involving analytic features.  Math. Program., 116(1):5–16, Jan 2009.
  * [4] Hedy Attouch, Jérôme Bolte, and Benar Fux Svaiter.  Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward–backward splitting, and regularized gauss–seidel methods.  Math. Program., 137(1):91–129, Feb 2013.
  * [5] Hédy Attouch, Jérôme Bolte, Patrick Redont, and Antoine Soubeyran.  Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the kurdyka-Łojasiewicz inequality.  Math. Oper. Res., 35(2):438–457, 2010.
  * [6] Gilles Aubert and Luminita Vese.  A variational method in image recovery.  SIAM J. Numer. Anal., 34(5):1948–1979, 1997.
  * [7] A. Auslender.  Asymptotic properties of the fenchel dual functional and applications to decomposition problems.  J. Optim. Theory Appl., 73(3):427–449, Jun 1992.
  * [8] Leah Bar, Nahum Kiryati, and Nir Sochen.  Image deblurring in the presence of impulsive noise.  Int. J. Comput. Vis., 70(3):279–298, Dec 2006.
  * [9] Michael J. Black and Anand Rangarajan.  On the unification of line processes, outlier rejection, and robust statistics with applications in early vision.  Int. J. Comput. Vis., 19(1):57–91, Jul 1996.
  * [10] A. Blake and A. Zisserman.  Visual Reconstruction.  The MIT Press, 1987.
  * [11] Jérôme Bolte, Aris Daniilidis, and Adrian Lewis.  The łojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems.  SIAM J. Optim., 17(4):1205–1223, 2007.
  * [12] Jérôme Bolte, Aris Daniilidis, Adrian Lewis, and Masahiro Shiota.  Clarke subgradients of stratifiable functions.  SIAM J. Optim., 18(2):556–572, 2007.
  * [13] Jérôme Bolte, Shoham Sabach, and Marc Teboulle.  Proximal alternating linearized minimization for nonconvex and nonsmooth problems.  Math. Program., 146(1):459–494, 2014.
  * [14] Kristian Bredies and Hong Peng Sun.  Preconditioned douglas–rachford algorithms for tv- and tgv-regularized variational imaging problems.  J. Math. Imaging Vis., 52(3):317–344, Jul 2015.
  * [15] Kristian Bredies and Hongpeng Sun.  Preconditioned douglas–rachford splitting methods for convex-concave saddle-point problems.  SIAM J. Numer. Anal., 53(1):421–444, 2015.
  * [16] Kristian Bredies and Hongpeng Sun.  A proximal point analysis of the preconditioned alternating direction method of multipliers.  J. Optim. Theory Appl., 173(3):878–907, 2017.
  * [17] Antonin Chambolle and Thomas Pock.  A remark on accelerated block coordinate descent for computing the proximity operators of a sum of convex functions.  Journal of computational mathematics, 1:29–54, 2015.
  * [18] Raymond Chan, Alessandro Lanza, Serena Morigi, and Fiorella Sgallari.  Convex non-convex image segmentation.  Numerische Mathematik, 138(3):635–680, 2018.
  * [19] Tony F Chan and Pep Mulet.  On the convergence of the lagged diffusivity fixed point method in total variation image restoration.  SIAM J. Numer. Anal., 36(2):354–367, 1999.
  * [20] Pierre Charbonnier.  Reconstruction d’image: régularisation avec prise en compte des discontinuités. PhD thesis, Université de Nice-Sophia Antipolis, 1994.
  * [21] Pierre Charbonnier, Laure Blanc-Feraud, Gilles Aubert, and Michel Barlaud.  Two deterministic half-quadratic regularization algorithms for computed imaging.  In Proceedings of 1st ICIP, volume 2, pages 168–172. IEEE, 1994\.
  * [22] Pierre Charbonnier, Laure Blanc-Féraud, Gilles Aubert, and Michel Barlaud.  Deterministic edge-preserving regularization in computed imaging.  IEEE T. Image Process., 6(2):298–311, 1997.
  * [23] Emilie Chouzenoux, Jérôme Idier, and Saïd Moussaoui.  A majorize–minimize strategy for subspace optimization applied to image restoration.  IEEE T. Image Process., 20(6):1517–1528, 2010.
  * [24] Frank H Clarke.  Optimization and nonsmooth analysis.  SIAM, 1990.
  * [25] Shengxiang Deng and Hongpeng Sun.  A preconditioned difference of convex algorithm for truncated quadratic regularization with application to imaging.  J. Sci. Comput., 88(2):1–28, 2021.
  * [26] David C Dobson and Curtis R Vogel.  Convergence of an iterative method for total variation denoising.  SIAM J. Numer. Anal., 34(5):1779–1791, 1997.
  * [27] Stuart Ganan and D McClure.  Bayesian image analysis: An application to single photon emission tomography.  Amer. Statist. Assoc, pages 12–18, 1985.
  * [28] Donald Geman and George Reynolds.  Constrained restoration and the recovery of discontinuities.  IEEE Trans. Pattern Anal. Mach. Intell., 14(3):367–383, 1992.
  * [29] Donald Geman and Chengda Yang.  Nonlinear image recovery with half-quadratic regularization.  IEEE T. Image Process., 4(7):932–946, 1995.
  * [30] Frank R Hampel, Elvezio M Ronchetti, Peter J Rousseeuw, and Werner A Stahel.  Robust statistics: the approach based on influence functions, volume 196.  John Wiley & Sons, 2011.
  * [31] Tom Hebert and Richard Leahy.  A generalized em algorithm for 3-d bayesian reconstruction from poisson data using gibbs priors.  IEEE Trans. Med. Imaging, 8(2):194–202, 1989.
  * [32] Michael Hintermüller, Steven-Marian Stengl, and Thomas M Surowiec.  Uncertainty quantification in image segmentation using the ambrosio–tortorelli approximation of the mumford–shah energy.  J. Math. Imaging Vis., pages 1–23, 2021.
  * [33] Jérôme Idier.  Convex half-quadratic criteria and interacting auxiliary variables for image restoration.  IEEE T. Image Process., 10(7):1001–1009, 2001.
  * [34] Askold G Khovanskiĭ.  Fewnomials, volume 88.  American Mathematical Soc., 1991.
  * [35] Christian Labat and Jérôme Idier.  Convergence of conjugate gradient methods with a closed-form stepsize formula.  J. Optim. Theory Appl., 136(1):43–60, 2008.
  * [36] Ta Lê Loi.  Łojasiewicz inequalities for sets definable in the structure r exp.  In Annales de l’institut Fourier, volume 45, pages 951–971, 1995\.
  * [37] Hoai An Le Thi and Tao Pham Dinh.  Difference of convex functions algorithms (dca) for image restoration via a markov random field model.  Optim. Eng., 18(4):873–906, 2017.
  * [38] Guoyin Li and Ting Kei Pong.  Calculus of the exponent of kurdyka–łojasiewicz inequality and its applications to linear convergence of first-order methods.  Found. Comput. Math., 18(5):1199–1232, 2018.
  * [39] Hao Li and Xiangxiong Zhang.  Superconvergence of high order finite difference schemes based on variational formulation for elliptic equations.  J. Sci. Comput., 82(2):1–39, 2020.
  * [40] Jean-Marie Lion and Jean-Philippe Rolin.  Théoreme de préparation pour les fonctions logarithmico-exponentielles.  In Annales de l’institut Fourier, volume 47, pages 859–884, 1997\.
  * [41] Boris S Mordukhovich.  Variational analysis and generalized differentiation I: Basic theory, volume 330.  Springer Science & Business Media, 2006.
  * [42] David Mumford and Agnès Desolneux.  Pattern theory: the stochastic analysis of real-world signals.  CRC Press, 2010.
  * [43] Mila Nikolova.  Markovian reconstruction using a gnc approach.  IEEE T. Image Process., 8(9):1204–1220, 1999.
  * [44] Mila Nikolova and Michael K Ng.  Analysis of half-quadratic minimization methods for signal and image recovery.  SIAM J. Sci. Comput., 27(3):937–966, 2005.
  * [45] Marc C Robini, Feng Yang, and Yuemin Zhu.  Inexact half-quadratic optimization for linear inverse problems.  SIAM J. Imaging Sci., 11(2):1078–1133, 2018.
  * [46] R Tyrrell Rockafellar and Roger J-B Wets.  Variational analysis, volume 317.  Springer Science & Business Media, 2009.
  * [47] Yousef Saad.  Iterative methods for sparse linear systems.  SIAM, 2003.
  * [48] Evgeny Strekalovskiy and Daniel Cremers.  Real-time minimization of the piecewise smooth mumford-shah functional.  In ECCV, pages 127–141. Springer, 2014.
  * [49] Meng Tang, Dmitrii Marin, Ismail Ben Ayed, and Yuri Boykov.  Kernel cuts: Kernel and spectral clustering meet regularization.  Int. J. Comput. Vis., 127(5):477–511, 2019.
  * [50] L. P. D. van den Dries.  Tame topology and o-minimal structures, volume 248.  Cambridge university press, 1998.
  * [51] Pieter Wesseling.  Introduction to multigrid methods. John Wiley & Sons, 1992.
  * [52] Alex J Wilkie.  Model completeness results for expansions of the ordered field of real numbers by restricted pfaffian functions and the exponential function.  J. Am. Math. Soc., 9(4):1051–1094, 1996.
  * [53] Gerhard Winkler.  Image analysis, random fields and Markov chain Monte Carlo methods: a mathematical introduction, volume 27.  Springer Science & Business Media, 2012.