Total variation regularization with variable Lebesgue prior

TL;DR

This paper introduces a convex total variation regularization method using a variable Lebesgue exponent, enabling local adaptation between edge preservation and smoothing, with efficient computation and promising results in denoising and tomography.

Contribution

It proposes a convex TV regularization with a fixed, spatially varying exponent function, improving flexibility and computational efficiency over previous variable exponent methods.

Findings

Favorable comparison to TGV and TV in denoising tasks

Efficient evaluation of proximal operators

Convex formulation with fixed variable exponent

Abstract

This work proposes the variable exponent Lebesgue modular as a replacement for the 1-norm in total variation (TV) regularization. It allows the exponent to vary with spatial location and thus enables users to locally select whether to preserve edges or smooth intensity variations. In contrast to earlier work using TV-like methods with variable exponents, the exponent function is here computed offline as a fixed parameter of the final optimization problem, resulting in a convex goal functional. The obtained formulas for the convex conjugate and the proximal operators are simple in structure and can be evaluated very efficiently, an important property for practical usability. Numerical results with variable $L^p$ TV prior in denoising and tomography problems on synthetic data compare favorably to total generalized variation (TGV) and TV.

Tables

Table 1. Table 1: Speed comparison of computing the proximal factor U τ subscript 𝑈 𝜏 U_{\tau} from 16 and the convex conjugate integrand R 𝑅 R as defined in 10 for 10 6 superscript 10 6 10^{6} points using 10 Newton iterations. The test machine has an 8-core Intel Core i7-6700K CPU and an NVidia GeForce GTX 1070 graphics card. The tested implementations are NumPy (NP), Cython (CY), Numba (NB) with CPU, parallel and CUDA targets, and libgpuarray. Numbers in parentheses are speed-ups compared to NumPy. Software versions: Python 3.5, NumPy 1.12, Cython 0.25.2, Numba 0.31.0, libgpuarray 0.6.0.

	NP	CY	NB CPU	NB par	NB CUDA	GPUArr
prox	1.33	0.95 (1.4)	1.06 (1.3)	0.29 (4.5)	0.39 (3.4)	0.028 (47)
cc	0.16	0.20 (0.8)	0.17 (0.9)	0.08 (1.8)	0.18 (0.86)	0.006 (26)

Equations

$$f_{\lambda }=\underset{f\in X}{\mathrm{argmin}}\left[\mathcal{D}(f)+\lambda \text{TV}(f)\right]$$

$${\text{TV}}(f)=\sup\Big{\{}\int_{\Omega}f\operatorname{div}\phi\,{\mathrm{d}}x\nonscript\>\Big{|}\nonscript\>\mathopen{}\allowbreak\phi\in L^{\infty}(\Omega,{\mathbb{R}}^{d}),\ \lVert\phi\rVert_{\infty}<1\Big{\}}$$

$$f^{\ast }=\mathrm{argmin}\left[\mathcal{D}(f)+\lambda {\int }_{\Omega }|\nabla f(x){|}^{p(\nabla f(x))}\mathrm{d}x\right],$$

$$p(\nabla f)=1+\frac{1}{1+k|G_{\sigma }\ast \nabla f{|}^2}.$$

$${\rho }_p(f):={\int }_{{\Omega }^{\ast }}|f(x){|}^{p(x)}\mathrm{d}x+\Vert f_{\infty }{\Vert }_{\infty },$$

$${\rho }_p(f)={\int }_{\Omega }|f(x){|}^{p(x)}\mathrm{d}x.$$

$$\Vert f{\Vert }_p:=\inf \{\lambda >0\text{nonscript}|\text{nonscript}{\rho }_p(f/\lambda )<1\}$$

$$\nabla {\rho }_p(f)=p|f{|}^{p-2}f.$$

$${\mathcal{S}}^{*}(f)\operatorname{{\coloneqq}}\sup_{g\in{X}}\big{[}\langle f,\,g\rangle_{X}-{\mathcal{S}}(g)\big{]},\quad f\in{X}.$$

$$M_{\tau }(\mathcal{S})(f):=\underset{g\in X}{\inf }\left[\mathcal{S}(g)+\frac{1}{2\tau }\Vert f-g{\Vert }_2^2\right].$$

$${\mathrm{prox}}_{\tau \mathcal{S}}(f):=\underset{g\in X}{\mathrm{arginf}}\left[\mathcal{S}(g)+\frac{1}{2\tau }\Vert f-g{\Vert }_2^2\right].$$

$${\mathcal{I}}_{F}\colon{X}\to({-\infty},{\infty}],\quad{\mathcal{I}}_{F}(u)\operatorname{{\coloneqq}}\int_{\Omega}F\big{(}x,u(x)\big{)}\,{\mathrm{d}}x$$

$$\underset{u\in X}{\inf }{\mathcal{I}}_F(u)={\int }_{\Omega }\left[\underset{z\in {\mathbb{R}}^m}{\inf }F(x,z)\right]\mathrm{d}x,$$

$$\bar{u}\in \underset{u\in X}{\mathrm{arginf}}{\mathcal{I}}_F(u)⟺\bar{u}(x)\in \underset{z\in X}{\mathrm{arginf}}F(x,z)\text{for almost all }x\in \Omega .$$

$$\rho_{p}^{*}(f)=\int_{\Omega}R\big{(}f(x),p(x)\big{)}\,{\mathrm{d}}x,$$

$$R(z,p)\operatorname{{\coloneqq}}\begin{cases}\iota_{\lvert\>\cdot\>\rvert\leq 1}(z),&\text{if }p=1,\\[4.30554pt] \lvert z\rvert^{2}/4,&\text{if }p=2,\\[4.30554pt] \lvert z\rvert^{\frac{p}{p-1}}\big{[}p^{-1/(p-1)}-p^{-p/(p-1)}\big{]},&\text{else},\end{cases}$$

$$⟨f,g{⟩}_X-{\rho }_p(g)={\int }_{\Omega }\left[⟨f(x),g(x){⟩}_{{\mathbb{R}}^m}-|g(x){|}^{p(x)}\right]\mathrm{d}x,$$

$$\rho_{p}^{*}(f)=\sup_{g\in{X}}\big{[}\langle f,\,g\rangle_{X}-\rho_{p}(g)\big{]}=\int_{\Omega}\sup_{y\in{\mathbb{R}}^{m}}\left[\langle f(x),\,y\rangle_{{\mathbb{R}}^{m}}-\lvert y\rvert^{p(x)}\right]{\mathrm{d}}x.$$

$$r^{*}(z)=\sup_{y\in{\mathbb{R}}^{m}}\big{[}\langle z,\,y\rangle_{{\mathbb{R}}^{m}}-\lvert y\rvert^{p}\big{]},$$

$$r^{*}(z)=\sup_{\alpha\geq 0}\sup_{\theta\in S^{m-1}}\big{[}\alpha\langle z,\,\theta\rangle_{{\mathbb{R}}^{m}}-\alpha^{p}\big{]}=\sup_{\alpha\geq 0}\big{[}\alpha(\lvert z\rvert-\alpha^{p-1})\big{]}.$$

$$v^{\prime }(\alpha )=|z|-p{\alpha }^{p-1},v^{\prime \prime }(\alpha )=-p(p-1){\alpha }^{p-2},$$

$$v(\bar{\alpha })=|z|{\left(\frac{|z|}{p}\right)}^{\frac{1}{p-1}}-{\left(\frac{|z|}{p}\right)}^{\frac{p}{p-1}}=|z{|}^{\frac{p}{p-1}}\left[p^{-1/(p-1)}-p^{-p/(p-1)}\right],$$

$$M_{\tau}(\rho_{p})(f)=\int_{\Omega}T_{\tau}\big{(}f(x),p(x)\big{)}\,{\mathrm{d}}x,$$

$$T_{\tau }(z,p):=\{\begin{array}{ll}\frac{|z{|}^2}{2\tau },& \text{if }p=1\text{ and }|z|\le \tau ,\\ |z|-\frac{\tau }{2},& \text{if }p=1\text{ and }|z|\ge \tau ,\\ \frac{|z{|}^2}{1+2\tau },& \text{if }p=2,\\ (|z|-\bar{\alpha })\frac{2\bar{\alpha }+p(|z|-\bar{\alpha })}{2\tau p},& \text{ else}.\end{array}$$

$$\alpha +\tau p{\alpha }^{p-1}=|z|.$$

$$M_{\tau }({\rho }_p)(f)$$

$$={\int }_{\Omega }\underset{y\in {\mathbb{R}}^m}{\inf }\left[|y{|}^{p(x)}+\frac{1}{2\tau }|f(x)-y{|}^2\right]\mathrm{d}x.$$

$$\displaystyle M_{\tau}\big{(}\lvert\>\cdot\>\rvert^{p}\big{)}(z)$$

$$\displaystyle=\inf_{\alpha\geq 0}\inf_{\theta\in S^{m-1}}\left[\alpha^{p}+\frac{1}{2\tau}\big{(}\lvert z\rvert^{2}-2\alpha\langle z,\,\theta\rangle_{{\mathbb{R}}^{m}}+\alpha^{2}\big{)}\right]$$

$$\displaystyle=\inf_{\alpha\geq 0}\left[\alpha^{p}+\frac{1}{2\tau}\big{(}\lvert z\rvert^{2}-2\alpha\lvert z\rvert+\alpha^{2}\big{)}\right]$$

$$\displaystyle=\inf_{\alpha\geq 0}\underbrace{\left[\alpha^{p}+\frac{1}{2\tau}\big{(}\lvert z\rvert-\alpha\big{)}^{2}\right]}_{\operatorname{{\eqqcolon}}v(\alpha)}.$$

$$v^{\prime }(\alpha )$$

Full text

1 Introduction

Total variation (TV) regularization is nowadays an established technique in imaging. It consists in finding

[TABLE]

for a data discrepancy measure ${\mathcal{D}}$, a regularization parameter $\lambda>0$ and the TV functional

[TABLE]

defined on the space ${\text{BV}}(\Omega)$ of functions of bounded variation on $\Omega\subset{\mathbb{R}}^{d}$. This model has proved to possess favorable properties when edge preservation and effective noise suppression are primary goals in image reconstruction. There exist very efficient implementations of numerical solvers for TV-regularized inverse problems, which also makes the technique accessible for daily practical use in engineering and medical imaging. Despite its success, the classical TV model has also exhibited its limitations over time, most notably the staircasing effect [26], i.e. the preference of the regularizer for piecewise constant functions.

Several extensions of total variation have been proposed to tackle the problem of staircasing, with the latest contribution being total generalized variation (TGV) [9], a higher-order TV variant. It has shown to perform significantly better than TV on images with both edges and gradual intensity changes, and it will serve as a benchmark for the method proposed in this paper. On the downside, TGV requires handling of higher order tensors leading to higher numerical complexity, and it introduces an additional regularization parameter, which makes their selection more complex.

This paper takes a different approach in varying classical TV regularization based on the observation that functionals of the form $\lVert\nabla\cdot\rVert_{p}^{p}$ lead to very different behavior of minimizers of $f\mapsto{\mathcal{D}}(f)+\lambda\lVert\nabla f\rVert_{p}^{p}$ for $p=1$ and $p>1$, with a sharp transition between the two cases: While the case $p=1$ corresponds to TV regularization, any exponent $p>1$ leads to a smooth solution [10]. Therefore, a natural approach would be to select the exponent $p$ dependent on spatial location, taking the value $p=1$ where edges should be enhanced, and $p>1$ otherwise.

The idea of using a spatially varying exponent in a TV-like regularization method has been proposed as early as 1997 [6] and put into practice in 2006 [12]. Both papers as well as subsequent articles try to tackle variants of the problem

[TABLE]

where the exponent depends directly on the image, e.g.,

[TABLE]

For the resulting non-convex problem, only partial solutions have been proposed, e.g., for a smoothed version of the integrand with a weak notion of solution [30].

A parallel track that has recently received growing interest are variable exponent sequence spaces and their application in imaging problems. The sequence usually represents coefficients in wavelet-type decompositions, and a sequence of exponents serves to define Besov-type spaces in this case [20]. A very recent paper studies convergence and stability properties of a variable exponent prior in this context [22]. That approach has no direct connection to the variable exponent functionals considered here, although the variable Lebesgue spaces can be characterized in terms of wavelet coefficients [28].

The method proposed in this paper computes the exponent $p$ in an offline step and keeps it as a fixed parameter in the then convex minimization problem. There are two natural types of imaging problems where this approach can be applied:

(1) single-channel imaging where first the exponent is computed from the given data and then applied as prior in the subsequent minimization problem (“bootstrapping”) and

(2) dual-channel imaging where the secondary channel provides the exponent map that is used for regularization of the primary channel.

We start by introducing the variable Lebesgue modular, the central functional in this work. The following section studies the properties of the modular in the light of convex optimization and derives expressions for the building blocks needed to solve the TV minimization problem introduced thereafter. A short section presents the method for computing the exponent map, before numerical results are shown.

2 The variable Lebesgue modular

Let $\Omega\subset{\mathbb{R}}^{d}$ be an open and connected domain. Let further $p\colon\Omega\to[{1},{\infty}]$ be a measurable exponent function. The variable Lebesgue modular [14] or $p$-modular of a function $f\colon\Omega\to{\mathbb{R}}$ is defined as

[TABLE]

where $\Omega^{*}=\{x\in\Omega\nonscript\>|\nonscript\>\mathopen{}\allowbreak p(x)<\infty\}$ and $f_{\infty}$ is the restriction to $\Omega_{\infty}=\Omega\setminus\Omega^{*}$. For vector-valued functions $F\colon\Omega\to{\mathbb{R}}^{m}$, we use the same notation, in which case $\lvert F(x)\rvert$ stands for the Euclidean norm in ${\mathbb{R}}^{m}$. Occasionally we denote this Euclidean norm with $\lvert\>\cdot\>\rvert_{2}$ if we deem it necessary for notational clarity.

We are mainly interested in exponents between 1 and 2, and therefore we assume that $p\colon\Omega\to[{1},{2}]$ and thus $\Omega=\Omega^{*}$. As a direct consequence, the modular simplifies to

[TABLE]

Remark**.**

Although the $p$-modular is a direct generalization of the $p$-th power of the norm $\lVert\,\cdot\,\rVert_{p}$ for constant exponent, it is not homogeneous and thus cannot be easily used to define a norm. To achieve this, the construction

[TABLE]

can be used. The corresponding space $L^{p}(\Omega)$ of functions where this norm is finite, as well as related Sobolev spaces, have been studied thoroughly [15, 14]. Recently also spaces of functions with bounded variation have attracted attention in the context of PDEs [18] and imaging [19].

3 Convex optimization with variable Lebesgue modular

We now consider the $p$-modular as a functional on $L^{2}(\Omega)$, where we assign the value $\infty$ to functions in $L^{2}(\Omega)\setminus L^{p}(\Omega)$. This identification is justified by the fact that for bounded exponent functions, the set of finite modular and finite norm are identical [14, Proposition 2.12].

In accordance with optimization literature, we denote by $\Gamma_{0}({X})$ the set of proper, convex and lower semicontinuous functionals on ${X}$. These minimum requirements are satisfied by the $p$-modular as a consequence of the integrand possessing these properties.

Proposition 3.1**.**

The $p$-modular is proper, convex and lower semicontinuous as a functional $\rho_{p}\colon{X}\to[{0},{\infty}]$.

It is similarly straightforward to see why this functional suffers from the same non-differentiability issue at 0 as the 1-norm, unless the exponent is equal to 2 everywhere.

Proposition 3.2**.**

The gradient of the $p$-modular is given by

[TABLE]

It is well-defined for all functions $f$ that are nonzero at almost all points $x$ where $p(x)\neq 2$.

As a consequence, smooth optimization cannot be applied unless the functional itself is replaced by a smoothed version, e.g., $\rho_{p}(\lvert f\rvert_{\epsilon})$, using the smoothed norm $\lvert f\rvert_{\epsilon}=(\lvert f\rvert^{2}+\epsilon^{2})^{1/2}$.

We do not follow this route and instead resort to non-smooth convex optimization. In the following we define the basice notions of convex conjugate, Moreau envelope and proximal operator. Note that all definitions can be formulated for general Banach spaces, but we focus on the special case of the Hilbert space $L^{2}$.

Definition 3.3**.**

Let ${X}=L^{2}(\Omega,{\mathbb{R}}^{m})$ with $m=1$ for the scalar case and $m>1$ for the vectorial case. Let further ${\mathcal{S}}\in\Gamma_{0}({X})$. The convex conjugate ${\mathcal{S}}^{*}\colon{X}\to({-\infty},{\infty}]$ of ${\mathcal{S}}$ is the functional

[TABLE]

The Moreau envelope or Moreau-Yosida regularization $M_{\tau}({\mathcal{S}})\colon{X}\to({-\infty},{\infty}]$ of ${\mathcal{S}}$ is defined as

[TABLE]

The proximal operator $\operatorname{prox}_{\tau{\mathcal{S}}}\colon{X}\to{X}$ of ${\mathcal{S}}$ is the mapping

[TABLE]

We derive these three fundamental mappings for the functional $\rho_{p}$ in the following. As an important tool, we state a simplified version of a theorem that allows to rephrase the optimization of an integrating functional as the pointwise optimization of its integrand.

Theorem 3.4**.**

[27, Theorem 14.60]** Let ${X}=L^{2}(\Omega,{\mathbb{R}}^{m})$ and $F\colon\Omega\times{\mathbb{R}}^{m}$ be a normal integrand, i.e., a function whose epigraphical mapping $E(x)=\mathrm{epi}F(x,\cdot)$ is closed-valued and measurable. Then, the integral functional

[TABLE]

has the property that minimization and integration can be exchanged, i.e.,

[TABLE]

as long as ${\mathcal{I}}_{F}\not\equiv+\infty$ on ${X}$. Likewise, unless this common value is $-\infty$, the minimizers of the functional satisfy

[TABLE]

The condition on the integrands to be normal is rather straightforward to check in the cases considered below. According to [27, Example 14.29], all Carathéordory integrands, i.e., functions that are measurable in the first variable and continuous in the second, are normal integrands. This characterization applies to all integrands in this paper.

Remark**.**

Some of the subsequent results can be pieced together from the literature, e.g., the convex conjugate of the norm, the Huber function [4] or expressions for the proximal operator [13] in certain cases. However, for the sake of completeness, and due to the coherent kind of argumentation, we provide complete proofs of all cases.

Lemma 3.5**.**

The convex conjugate of $\rho_{p}$ is given by

[TABLE]

where

[TABLE]

for $z\in{\mathbb{R}}^{m}$ and $p\in[{1},{2}]$. Here, $\iota_{\lvert\>\cdot\>\rvert\leq 1}$ denotes the indicator function of the unit ball in ${\mathbb{R}}^{m}$, i.e., the function that takes the value 0 in the set and $\infty$ outside.

Proof**.**

Since $\sup A=-\inf(-A)$ for any set $A\subset{\mathbb{R}}$, 8 holds also for $\sup$ instead of $\inf$. For the $p$-modular, the functional to be maximized can be written as

[TABLE]

which means that 7 can be applied, i.e., integration and optimization can be exchanged:

[TABLE]

Hence, the optimization of the integrand can be performed for each point $x\in\Omega$ separately. In what follows, we consider $p=p(x)$ to be a constant and set $z=f(x)\in{\mathbb{R}}^{m}$. The new problem is to find

[TABLE]

which is nothing but the convex conjugate of $r=\lvert\>\cdot\>\rvert^{p}$. Apparently, $r^{*}(0)=0$. For other points, using polar coordinates, this problem can be rephrased as

[TABLE]

If $p=1$, the expression in brackets is only bounded from above if $\lvert z\rvert\leq 1$, in which case the supremum is [math]. Thus, the convex conjugate is $r^{*}=\iota_{\lvert\>\cdot\>\rvert\leq 1}$, the indicator function of the unit ball.

Otherwise, we can differentiate $v(\alpha)=\alpha\lvert z\rvert-\alpha^{p}$ with respect to $\alpha$ and get

[TABLE]

which shows that $\bar{\alpha}=(\lvert z\rvert/p)^{1/(p-1)}$ is a critical value of $v$, and further that $v$ is strictly concave and the critical value is the only local maximum. The maximum value is

[TABLE]

which simplifies to $v(\bar{\alpha})=\lvert z\rvert^{2}/4$ for $p=2$. Since this value is larger than $v(0)=0$, it is also the global maximum.

For the Moreau envelope, we can largely use the same technique to find pointwise expressions.

Lemma 3.6**.**

The Moreau envelope of $\rho_{p}$ can be written as

[TABLE]

where

[TABLE]

Here, $\bar{\alpha}=\bar{\alpha}(\lvert z\rvert,p,\tau)$ denotes the unique solution $0\leq\bar{\alpha}<\lvert z\rvert$ to the equation

[TABLE]

Proof**.**

We have, again due to 7, that

[TABLE]

Using the same notation as in the proof of Lemma 3.5, we can rewrite the inner optimization problem as the the determination of

[TABLE]

Differentiation of $v$ yields

[TABLE]

Thus the only critical point of $v$ is given as the solution to the implicit equation

[TABLE]

which is the same as 13. In the special case $p=2$, this equation is easily solved for $\alpha$, resulting in $\bar{\alpha}=\lvert z\rvert/(1+2\tau)$. The positivity of $v^{\prime\prime}$ lets us conclude that this point is the only local minimum of $v$. For the minimum value, we get

[TABLE]

This value is clearly smaller than $v(0)=\lvert z\rvert^{2}/(2\tau)$, therefore $\bar{\alpha}$ is also the global minimum.

For $p=1$, on the other hand, the critical point lies at $\bar{\alpha}=\lvert z\rvert-\tau$. Here we must distinguish two cases:

1. $\lvert z\rvert<\tau$:

The critical point is negative and thus not eligible as minimum. To find the global minimum, we observe that

[TABLE]

from which immediately follows that the global minimum is at $\alpha=0$, with minimum value $v(\bar{\alpha})=\lvert z\rvert^{2}/(2\tau)$. 2. $\lvert z\rvert\geq\tau$:

Now the critical point $\bar{\alpha}=\lvert z\rvert-\tau$ is valid, and the positivity of the second derivative $v^{\prime\prime}(\alpha)=1/\tau$ classifies this point as a local minimum. The corresponding function value is

[TABLE]

Finally, for $1<p<2$, there is no general closed form solution to 13, except for $z=0$ with solution $\alpha=0$. However, we can use that $v^{\prime\prime}$ is positive, which implies that any critical point is a local minimum. Further, if $z\neq 0$, then it can be readily deduced that the left-hand side of 13 is smaller than $\lvert z\rvert$ for $\alpha=0$ and larger than $\lvert z\rvert$ for $\alpha=\lvert z\rvert$. Therefore, by continuity it follows that there must be a solution $0<\bar{\alpha}<\lvert z\rvert$, which must be the global minimum due to the global convexity of $v$. To simplify the expression for the minimum value, we multiply 13 with $\alpha$ and solve for $\alpha^{p}$:

[TABLE]

Inserting this identity into $v$ with $\alpha=\bar{\alpha}$ finally yields

[TABLE]

The proximal operator can now be easily deduced from the above proof.

Corollary 3.7**.**

The proximal operator $\operatorname{prox}_{\tau\rho_{p}}\colon{X}\to{X}$ is given as

[TABLE]

with the interpretation $\operatorname{prox}_{\tau\rho_{p}}(f)(x)=0$ for $f(x)=0$, where

[TABLE]

and $\bar{\alpha}$ is defined as in 13.

Proof**.**

The values of the function $U_{\tau}$ correspond precisely to the minima as calculated in the proof of 13. For the directional part $f(x)/\lvert f(x)\rvert$, it suffices to observe that in 14,

[TABLE]

is achieved in the unique vector $\bar{\theta}=z/\lvert z\rvert$.

Often optimization algorithms require the proximal operator of the convex conjugate of a given functional. Thanks to the Moreau decomposition, this operator can be deduced from the proximal of the functional itself.

Corollary 3.8**.**

The proximal operator $\operatorname{prox}_{\tau\rho_{p}^{*}}$ of the convex conjugate is given by

[TABLE]

with

[TABLE]

As before, $\bar{\alpha}=\bar{\alpha}(\lvert z\rvert,p,\tau)$ is implicitly defined as solution to 13.

Proof**.**

Since ${\mathcal{S}}={\mathcal{S}}^{**}$ for all ${\mathcal{S}}\in\Gamma_{0}$, we can write the Moreau decomposition [4, Theorem 14.3] for $\rho_{p}^{*}$ as

[TABLE]

This implies that

[TABLE]

hence by the substitutions $z=f(x)$ and $p=p(x)$, the term in brackets can for $p=1$ be written as

[TABLE]

Likewise, for $p=2$, we have

[TABLE]

for the bracketed expression. In the case $1<p<2$, we apply the same substitution to 13 and get

[TABLE]

using the new variable $\beta=\tau\alpha$. This equation is again of the form 13 with $\tau^{1-p}$ in lieu of $\tau$.

Solving the implicit equation

Since the Moreau envelope and the proximal operators involve solving a non-linear equation for each point where $p$ is between 1 and 2, we formulate and analyze in the following a numerical scheme for its solution. Newton’s method seems well-suited for this task since the function whose zero is to be determined is smooth and strictly convex. The largest practical issue is the determination of a suitable start point that guarantees valid iterates and provides fast convergence to the solution.

Lemma 3.9**.**

The Newton iteration for solving 13 takes the form

[TABLE]

If the start value $\alpha_{0}$ satisfies

[TABLE]

the sequence $(\alpha_{k})_{k\geq 1}$ is positive and strictly increasing, and converges quadratically to the unique root $\bar{\alpha}$.

Proof**.**

We aim at finding a root of the function $s(\alpha)=\alpha+\tau p\alpha^{p-1}-\lvert z\rvert$ whose derivative is given by $s^{\prime}(\alpha)=1+\tau p(p-1)\alpha^{p-2}$. As noted already in the proof of 13, $s$ is strictly increasing and has exactly one root between $\alpha=0$ and $\alpha=\lvert z\rvert$, provided that $z\neq 0$.

For the Newton iterates, we have

[TABLE]

The first iterate $\alpha_{1}$ is positive if the numerator $\lvert z\rvert+\tau p(p-2)\alpha_{0}^{p-1}$ is positive. This is equivalent to 20.

We now show that the sequence $\alpha_{k}$ is strictly increasing from $k=1$ on, which ensures positivity in case $\alpha_{1}$ is positive. The function $s$ is twice differentiable, strictly concave and strictly increasing. Using Taylor expansion around the unique root $\bar{\alpha}$, we have

[TABLE]

for some $\eta_{k}$ between $\bar{\alpha}$ and $\alpha_{k}$. Dividing by $s^{\prime}(\alpha_{k})>0$ and rearranging terms yields

[TABLE]

i.e. the approximation error $e_{k}=\bar{\alpha}-\alpha_{k}$ satisfies

[TABLE]

Besides the quadratic convergence of Newton’s method, this also shows that $e_{k+1}>0$ for all $k\geq 0$ since $s^{\prime\prime}<0$ and $s^{\prime}>0$ everywhere. Finally, from the iteration we see that

[TABLE]

since $\alpha_{k}$ is smaller than the root $\alpha$, and $s^{\prime}>0$. This concludes the proof.

Lemma 3.9 formulates a necessary condition that a start value has to meet in order to produce a valid sequence of iterates. However, the bound does not imply any practical rule how to choose the value and may, in fact, yield impractically large or small numbers. This situation occurs when $p$ is close to 1 since the resulting power $1/(p-1)$ applied to the right-hand side of 20 is very large and yields either huge or tiny numbers, none of which are useful.

To find a better strategy for the selection of a start value, we repeat the observation from the proof of 11 that $\bar{\alpha}<\lvert z\rvert$, which will act as a second bound on the start value. Finally, since the solutions $\alpha^{(1)}$ and $\alpha^{(2)}$ for the extreme cases $p=1$ and $p=2$, respectively, can be computed explicitly, and since $\alpha^{(2)}$ is positive for $z\neq 0$, we choose the convex combination

[TABLE]

as an initial guess for a start value and set

[TABLE]

as the effective start value.

4 TV regularization with variable exponent

Let ${\mathcal{T}}\colon{X}\to{Y}$ be a bounded linear operator, where ${X}=L^{2}(\Omega)$ for $\Omega\subset{\mathbb{R}}^{d}$ and ${Y}$ is another Hilbert space. We consider the inverse problem of finding $f\in{X}$ such that ${\mathcal{T}}(f)=g$ for a given $g\in{Y}$. TV regularization in the context of $L^{2}$ spaces can be formulated as the problem of finding

[TABLE]

with a given Fréchet-differentiable data discrepancy functional ${\mathcal{D}}(\cdot;g)$ and the total variation

[TABLE]

We now replace the 1-norm with the $p$-modular for a given exponent function $p\colon\Omega\to[{1},{2}]$, leading to the optimization problem

[TABLE]

where the $p$-TV functional is defined as ${\text{TV}}^{p}=\rho_{p}\circ\nabla$. The functional ${\text{TV}}^{p}$ imposes two kinds of prior knowledge onto a given inverse problem to determine $f$ from $g={\mathcal{T}}(f)$. Firstly, it locally enforces either a sparse gradient ($p=1$) or a smooth solution ($p=2$), depending on the exponent value at the given location. Secondly, the exponent itself contains prior knowledge as to which of the two constraints is appropriate to enforce at a specific location. Hence it is important to find a robust and reliable way to construct an exponent function from given data. This aspect is investigated further in Section 5.

5 Construction of the exponent function

In the formulation of the problem 25, the exponent $p\colon\Omega\to[{1},{2}]$ is assumed to be known. For a concrete application this means that the function needs to be determined beforehand by a method providing a reasonable distinction between edge features corresponding to $p=1$ and non-edges that are mapped to other exponent values. Due to numerical instability, it is not feasible to allow exponents arbitrarily close to 1, as seen in 20.

In general, given a function $f\in{X}$, we construct the exponent function in the following way:

Compute the smoothed Laplacian $l=\varDelta(G_{\sigma_{1}}\ast f)$, where $G_{\sigma}$ is a Gaussian of width $\sigma$. The smoothing kernel should suppress noise and be chosen rather narrow. 2. 2.

Take the absolute value and smooth again to make the region of detected edges larger: $a=G_{\sigma_{2}}\ast\lvert l\rvert\colon\Omega\to[{0},{\infty})$. Typically, this convolution kernel should be significantly wider, i.e., $\sigma_{2}>\sigma_{1}$. 3. 3.

Multiply the result with a constant $c$ and threshold at 1, thus effectively clipping at the value $c$: $t=\min\{c\cdot a,\ 1\}\colon\Omega\to[{0},{1}]$. 4. 4.

The final exponent is given by $p=2-c\colon\Omega\to[{1},{2}]$.

We consider two scenarios of generating the exponent prior from given data: “bootstrapping” and bimodal imaging.

Bootstrapping

This technique is applied when only a single dataset is available. In this case, the first step consists in computing an image from the data that can be used in the procedure for extracting an exponent as shown in Figure 1. For denoising problems, the data is already in image space, and no action is required. In case of tomography, a simple method like filtered back-projection can be used to gain the initial image. After that, the procedure as described above is applied to the image, yielding an exponent function that encodes likely occurrences of edges in the prior image. This step can be enhanced by directly reconstructing edge images by Lambda tomography [17, 25] or Approximate Inverse for feature reconstruction [23, 24].

The exponent acquired in this way can then be used to define the variable $L^{p}$ modular that serves as a building block for the ${\text{TV}}^{p}$ regularizer in the minimization problem 25.

Bimodal imaging

A regular situation in bimodal imaging, i.e., when the same object is imaged with two different modalities, is that one type of data is much more reliable than the other in terms of signal-to-noise ratio, resolution or similar meaures of quality.

Typical examples are PET-MRI [16], PET-CT [21] and SPECT-CT [3]. In all these applications, the actual quantity of interest, a radioactive tracer distribution, cannot be reconstructed without a secondary source of information about the attenuation of the object. This secondary parameter can be obtained at high resolution, while the primary channel is usually characterized by strong noise and poor spatial resolution. Nevertheless, edges from the secondary channel may still be valuable prior information for the reconstruction of the primary quantity of interest, but should not be strictly enforced.

The variable exponent in the ${\text{TV}}^{p}$ prior seems to be well-suited for this scenario – it suggests edges without enforcing them. Here, a reconstruction of the secondary quantity is used as input $f$ to the exponent calculation, and the resulting function $p$ is then used in the regularized problem 25 for the primary channel alone.

Remarks**.**

While the exponent prior is not very prone to creating new edges in the primary channel, it may very well lead to edges that are only present in the primary channel being missed by the proposed method. Therefore, it should likely be complemented by another reconstruction method that is not biased in this sense. 2. 2.

The recipe for generating the exponent $p$ contains 3 parameters $\sigma_{1}$, $\sigma_{2}$ and $c$, which may seem like an impractically large parameter space. However, there are two aspects which greatly simplify parameter selection, namely

(1) that the range of sensible choices is very limited and can be made solely based on the size of the image domain and the level of noise, and

(2) that different settings can be tested with immediate and visual feedback, due to the simplicity of the involved operations.

This is confirmed by the fact that the examples in Section 6 use nearly identical settings.

6 Numerical results

In this section we study the numerical characteristics of the proposed regularization functional. First we demonstrate that structure of the variable $L^{p}$ modular itself and its derived properties allow a very efficient implementation, thus allowing the functional to be applied in large-scale inverse problems. We then test the effect of the ${\text{TV}}^{p}$ functional in isolation by using it in several single-channel image denoising problems. Finally, the method of applying the exponent prior from a secondary channel to regularize a primary channel is evaluated in a simulated bimodal tomography example. The code for all examples is available at GitHub111https://github.com/kohr-h/odl/tree/variable_lp.

6.1 Implementation

The code for the evaluation of the involved functionals and operators is implemented in the ODL (Operator Discretization Library) framework [1] for inverse problems. It is written in the Python programming language and features classes and functions for handling vector spaces, operators and functionals in a transparent and efficient way, offering a syntax close to mathematical notation. The library also contains a number of optimization methods that are adequate for solving convex problems with non-differentiable functionals and compositions with linear operators, among others the Chambolle-Pock method [11] and two methods based on Forward-Backward splitting [7] and Douglas-Rachford splitting [8], respectively. Hence, the implementation of the functionals $\rho_{p}$ and $\rho_{p}^{*}$ along with their proximal operators is sufficient to enable numerical tests with regularization using the ${\text{TV}}^{p}$ functional.

For running time comparison, the proximal operators are implemented in four variants:

(1) with NumPy222http://www.numpy.org/ [29], a package for fast vectorized array computations,

(2) in the Cython333http://cython.org/ language [5], a superset of Python with types that is compiled as native C code,

(3) using Numba444http://numba.pydata.org/, a just-in-time compiler to accelerate pure Python functions, and

(4) as GPU kernels using libgpuarray555https://github.com/Theano/libgpuarray, a subproject of the deep learning framework Theano [2] providing NumPy-like arrays on the GPU.

Table 1 shows a comparison of running times of two algorithms that are relevant for regularization with the variable $L^{p}$ modular. The first one is the computation of the proximal factor $U_{\tau}$ as defined in 16, which involves a Newton iteration as detailed in Lemma 3.9. In the second test, the integrand $R$ of the convex conjugate as given in 10 is computed for the same number of points. Clearly, the GPU implementation gives a significant speed-up in both cases ($47\times$/$26\times$), while the other acceleration techniques only provide moderate speed gains in the first scenario ($1.4\times$ to $4.5\times$) or even result in slower execution in the second case. This behavior can mainly be attributed to the fact that the function $R$ is too simple in structure and requires too little work for those techniques to pay off, while the evaluation of $U_{\tau}$ involves a Newton iteration that clearly benefits from acceleration.

6.2 Denoising examples

We start out by testing ${\text{TV}}^{p}$ regularization in a pure denoising context, to keep concerns separate from issues with forward operators. The first test case is an image whose grey values given as

[TABLE]

on the domain $\Omega=[-10,10]^{2}$. It serves as a simple prototype of an image with sharp edges and gradual intensity variations. The noisy data is generated by adding standard white noise scaled by 0.1 times the dynamic range of the image. We will henceforth refer to this as “10 % noise”.

For comparision, we denoise the image using two other priors, namely classical TV and $\operatorname{TGV}^{2}$ [9]. The latter is implemented by taking the first component $f^{*}$ of the solution to the product space problem

[TABLE]

with the component-wise gradient $\mathcal{E}(v)=(\nabla v_{1},\dots,\nabla v_{d})\colon\Omega\to{X}^{d\times d}$.

The results in this test scenario are show in Figure 2. Clearly, TGV excels at the given problem since the true image falls into the class of piecewise linear functions that this method tries to promote. On the other hand, TV exhibits the typical staircasing effect, which makes it a poor fit for this particular phantom.

Looking at the results with ${\text{TV}}^{p}$ prior, one can make some interesting observations. First, the staircasing problem in the interior of the square does not show here since the exponent map takes values between 1.5 and 2.0, resulting in the smoothing behavior that is also typical for the Tikhonov-type prior $\lVert\nabla f\rVert_{2}^{2}$. Second, moderate variations of the exponent due to noise have little influence on the values of the reconstruction (see also Figure 4). Thirdly, the transition from $p=1$ to $p>1$ does not introduce artificial edges. Any artificial boundaries stem from the TV prior itself that is effective in regions where $p=1$.

On the downside, the variable $L^{p}$ regularization blurs the top and bottom edges in the middle where the contrast to the background is low, which is to be expected given the exponent map not capturing these edges. A more fine-tuned exponent calculation would help alleviate this issue. Further, the left and right edges appear slightly jagged, compared to the other two alternatives. This issue can likely be tackled by choosing a slightly larger regularization parameter.

The second test scenario consists in denoising of a natural image with more complex edge structure and intensity variation. The test image is taken to be the transformation to greyscale of the image generated by the SciPy666https://scipy.org/ command scipy.misc.face(), and its noisy variant contains 15 % white noise.

Results for this test case are shown in Figure 3. It is remarkable that the difference between TV and TGV are only marginal. While the exact reason for this behavior remains unclear, it can be said that the parameter $\lambda_{2}$ in 26 has practically no influence even if varied by 10 orders of magnitude. Hence, any lowering of the parameter value for $\lambda_{1}$ – to reduce staircasing – below the value used in the TV method results in insufficient noise suppression.

By contrast, the ${\text{TV}}^{p}$ approach works very well in this scenario and produces a visually pleasing and naturally looking result with smooth intensity variations in the regions where they ought to be expected, while preserving most of the sharp edges. The exact behavior of the method can easily be interpreted from the generated exponent map.

6.3 Bimodal tomography examples

This application scenario considers tomography with two datasets generated from the same phantom (taken from [9]). It is intended to indicate the potential of variable $L^{p}$ TV regularization in bimodal imaging, not to model realistic imaging conditions.

The datasets are the 2D divergent beam transform of a simple digital phantom corrupted by low noise (1 %) in the secondary and high noise (15 %) in the primary imaging channel, respectively. In Figure 5, both TV and TGV reconstructions make use of the more noisy primary channel only, while ${\text{TV}}^{p}$ takes the secondary, less noisy data into account for computing the exponent map. For a fair comparison, Figure 6 shows the result of ${\text{TV}}^{p}$ regularization with the “bootstrapping” approach of using only the primary channel. In both cases, the exponent function $p$ is computed from a filtered back-projection reconstruction with the approach described in Section 5.

As expected, the TV reconstruction recovers the edges well at the cost of strong staircasing within the gradually varying regions. TGV, on the other hand, removes the staircasing, at the cost of turning some edges into linear slopes. Clearly, ${\text{TV}}^{p}$ yields the visually best result when the secondary channel information is taken into account in form of an accurate exponent map. Only at the boundaries some parts of the edges are blurred due to the regions of $p=1$ not extending all the way to the boundaries, a consequence of the convolution used for smoothing. When taking only the primary channel into account for both exponent and reconstruction, the result contains a larger amount of the blocky structures of the TV reconstruction since the regions with $p=1$ are overestimated. More sophisticated methods for finding edges would potentially help alleviate this issue.

7 Future work

While the ${\text{TV}}^{p}$ prior using the modular $\rho_{p}$ leads to very encouraging results in the test scenarios considered in this paper, it has the potential drawback of not being scale-invariant. In other words, the absolute magnitude of the unknown image may not only have an impact on the choice of regularization parameter, but also on the relative importance of regions with different exponents to the functional value. This issue can be solved by replacing the variable $L^{p}$ modular with its corresponding norm 3. Furthermore, the $p$-norm seems to be better suited for the derivation of guarantees for convergence and stability. This topic, along with efficient ways of computing the norm, its convex conjugate and its proximal operator, will be the subject of future work.

Acknowledgments

Part of this work was financially supported by the Netherlands Organization for Scientific Research (NWO), project 639.073.506.