Minimizing $L_{1}$ over $L_{2}$ norms on the gradient

Regularization methods play an important role in inverse problems to refine the solution space by prior knowledge and/or special structures. For example, the celebrated total variation (TV) [1] prefers piecewise constant images, while total generalized variation (TGV) [2] and fractional-order TV [3, 4] tend to preserve piecewise smoothness of an image. TV can be defined either isotropically or anisotropically. The anisotropic TV [5] in the discrete setting is equivalent to applying the $L_{1}$ norm on the image gradient. As the $L_{1}$ norm is often used to enforce a signal being sparse, one can interpret the TV regularization as to promote the sparsity of gradient vectors.

To find the sparsest signal, it is straightforward to minimize the $L_{0}$ norm (counting the number of nonzero elements), which is unfortunately NP-hard [6]. A popular approach involves the convex relaxation of using the $L_{1}$ norm to replace the ill-posed $L_{0}$ norm, with the equivalence between $L_{1}$ and $L_{0}$ for sparse signal recovery given in terms of restricted isometry property (RIP) [7]. However, Fan and Li [8] pointed out that the $L_{1}$ approach is biased towards large coefficients, and proposed to minimize a nonconvex regularization, called smoothly clipped absolute deviation (SCAD). Subsequently, various nonconvex functionals emerged such as minimax concave penalty (MCP) [9], capped $L_{1}$ [10, 11, 12], and transformed $L_{1}$ [13, 14, 15]. Following the literature on sparse signal recovery, there is a trend to apply a nonconvex regularization on the gradient to deal with images. For instance, Chartrand [16] discussed both the $L_{p}$ norm with $0<p<1$ for sparse signals and $L_{p}$ on the gradient for magnetic resonance imaging (MRI), while MCP on the gradient was proposed in [17].

Recently, a scale-invariant functional $L_{1}/L_{2}$ was examined, which gives promising results in recovering sparse signals [18, 19, 20] and sparse gradients [21]. In this paper, we rely on a constrained formulation to characterize some scenarios, under which the quotient of the $L_{1}$ and $L_{2}$ norms on the gradient performs well. In particular, we borrow the analysis of a super-resolution problem, which refers to recovering a sparse signal from its low-frequency measurements. Candés and Fernandez-Granda [22] proved that if a signal has spikes (locations of nonzero elements) that are sufficiently separated, then the $L_{1}$ minimization yields an exact recovery for super-resolution. We innovatively design a certain type of piecewise constant signals that lead to well-separated spikes after taking the gradient. Using such signals, we empirically demonstrate that the TV minimization can find the desired solution under a similar separation condition as in [22]. We also illustrate that $L_{1}/L_{2}$ can deal with less separated spikes in gradient, and is better at preserving image contrast than $L_{1}$. These empirical evidences show $L_{1}/L_{2}$ holds great potentials in promoting sparse gradients and preserving image contrasts. To the best of our knowledge, it is the first time to relate the exact recovery of gradient-based methods to minimum separation and image contrast in a super-resolution problem.

Numerically, we consider the same splitting scheme used in an unconstrained formulation [21] to minimize the $L_{1}/L_{2}$ on the gradient, followed by the alternating direction method of multipliers (ADMM) [23]. We formulate the linear constraint using an indicator function, which is not strongly convex. As a result, the convergence analysis in the unconstrained model [21] is not directly applicable to this problem. We utilize the property of indicator function as well as the optimality conditions for constrained optimization problems to prove that the sequence generated by the proposed algorithm has a subsequence converging to a stationary point. Under a stronger assumption, we can establish the convergence of the entire sequence, referred to as global convergence.

We present some algorithmic insights on computational efficiency of our proposed algorithm for nonconvex optimization. In particular, we show an additional box constraint not only prevents the solution from being stuck at local minima, but also stabilizes the algorithm. Furthermore, we discuss algorithmic behaviors on two types of applications: MRI and computed tomography (CT). For the MRI reconstruction, a subproblem in ADMM has a closed-form solution, while an iterative solver is required for CT. As the accuracy of the subproblem varies between MRI and CT, we shall alter internal settings of our algorithm accordingly. In summary, this paper relies on a constrained formulation to discuss theoretical and computational aspects of a nonconvex regularization for imaging problems. The major contributions are three-fold:

1. (i)

   We reveal empirical evidences towards exact recovery of piecewise constant signals and demonstrate the superiority of $L_{1}/L_{2}$ on the gradient over TV.

2. (ii)

   We establish the subsequential convergence of the proposed algorithm and explore the global convergence under the certain assumptions.

3. (iii)

   We conduct extensive experiments to characterize computational efficiency of our algorithm and discuss how internal settings can be customized to cater to specific imaging applications, such as MRI and limited-angle CT reconstruction. Numerical results highlight the superior performance of our approach over other gradient-based regularizations.

The rest of the paper is organized as follows. Section 2 defines the notations that will be used through the paper, and gives a brief review on the related works. The empirical evidences for TV’s exact recovery and advantages of the proposed model are given in Section 3. The numerical scheme is detailed in Section 4, followed by convergence analysis in Section 5. Section 6 presents three types of imaging applications: super-resolution, MRI and CT reconstruction problems. Finally, conclusions and future works are given in Section 7.

We use a bold letter to denote a vector, a capital letter to denote a matrix or linear operator, and a calligraphic letter for a functional space. We use $\odot$ to denote the component-wise multiplication of two vectors. When a function (e.g., sign, max, min) applies to a vector, it returns a vector with corresponding component-wise operation.

We adopt a discrete setting to describe the related models. Suppose a two-dimensional (2D) image is defined on an $m\times n$ Cartesian grid. By using a standard linear index, we can represent a 2D image as a vector, i.e., the $((i-1)m+j)$-th component denotes the intensity value at pixel $(i,j).$ We define a discrete gradient operator,

where $D_{x},D_{y}$ are the finite forward difference operator with periodic boundary condition in the horizontal and vertical directions, respectively. We denote $N:=mn$ and the Euclidean spaces by $\mathcal{X}:=\mathds{R}^{N},\mathcal{Y}:=\mathds{R}^{2N}$, then $\mathbf{u}\in\mathcal{X}$ and $D\mathbf{u}\in\mathcal{Y}$. We can apply the standard norms, e.g., $L_{1},L_{2},$ on vectors $\mathbf{u}$ and $D\mathbf{u}.$ For example, the $L_{1}$ norm on the gradient, i.e., $\|D\mathbf{u}\|_{1}$, is the anisotropic TV regularization [5]. Throughout the paper, we use TV and “$L_{1}$ on the gradient” interchangeably. Note that the isotropic TV [1] is the $L_{2,1}$ norm, i.e., $\|(D_{x}\mathbf{u},D_{y}\mathbf{u})^{T}\|_{2,1},$ although Lou et al. [24] claimed to consider a weighted difference of anisotropic and isotropic TV based on the $L_{1}$-$L_{2}$ functional [25, 26, 27, 28] (isotropic TV is not the $L_{2}$ norm on the gradient.)

We examine the $L_{1}/L_{2}$ penalty on the gradient in a constrained formulation,

One way to solve for (2) involves the following equivalent form

with two auxiliary variables $\mathbf{d}$ and $\mathbf{h}$. For more details, please refer to [18] that presented a proof-of-concept example for MRI reconstruction. Since the splitting scheme (3) involves two block variables of $\mathbf{u}$ and $(\mathbf{d},\mathbf{h})$, the existing ADMM convergence results [29, 30, 31] are not applicable. An alternative approach was discussed in our preliminary work [21] for an unconstrained minimization problem,

where $\lambda>0$ is a weighting parameter. By only introducing one variable $\mathbf{h},$ the new splitting scheme (4) can guarantee the ADMM framework with subsequential convergence.

In this paper, we incorporate the splitting scheme (4) to solve the constrained problem (2), which is crucial to reveal theoretical properties of the gradient-based regularizations for image reconstruction, as elaborated in Section 3. Another contribution of this work lies in the convergence analysis, specifically for different optimality conditions of the constrained problem, as opposed to unconstrained formulation presented in [21]. It is true that the constrained formulation limits our experimental design in a noise-free fashion, but it helps us to draw conclusions solely on the model, ruling out the influence from other nuisances such as noises and tuning parameters. Our model (2) is parameter-free, while there is a parameter $\lambda$ in the unconstrained problem (4).

We aim to demonstrate the superiority of $L_{1}/L_{2}$ on the gradient over TV for a super-resolution problem [32], in which a sparse vector can be exactly recovered via the $L_{1}$ minimization. A mathematical model for super-resolution is expressed as

where $i$ is the imaginary unit, $\mathbf{u}\in\mathbb{R}^{N}$ is a vector to be recovered, and $\mathbf{b}\in\mathbb{C}^{n}$ consists of the given low frequency measurements with $n=2f_{c}+1<N$. Recovering $\mathbf{u}$ from $\mathbf{b}$ is referred to as super-resolution in the sense that the underlying signal $\mathbf{u}$ is defined on a fine grid with spacing $1/N$, while a direct inversion of $n$ frequency data yields a signal defined on a coarser grid with spacing $1/n$. For simplicity, we use matrix notation to rewrite (5) as $\mathbf{b}=S_{n}F\mathbf{u}$, where $S_{n}$ is a sampling matrix that collects the required low frequencies and $F$ is the Fourier transform matrix. A sparse signal can be represented by $\mathbf{u}=\sum_{j\in T}c_{j}\mathbf{e}_{j},$ where $\mathbf{e}_{j}$ is the $j$-th canonical basis in $\mathbb{R}^{N}$, $T$ is the support set of $\mathbf{u}$, and $\{c_{j}\}$ are coefficients. Following the work of [32], the sparse spikes are required to be sufficiently separated to guarantee the exact recovery of the $L_{1}$ minimization. To make the paper self-contained, we provide the definition of minimum separation in Definition 1 and an exact recovery condition in Theorem 1.

We empirically extend the analysis from sparse signals to sparse gradients. For this purpose, we construct a one-bar step function of length 100 with the first and the last $s$ elements taking value 0, and the remaining elements equal to 1, as illustrated in Figure 1. The gradient of such signal is 2-sparse with MS to be $\min(2s,100-2s)$ due to wrap-around distance. By setting $f_{c}=2,$ we only take $n=5$ low frequency measurements, and reconstruct the signal by minimizing either $L_{1}$ or $L_{1}/L_{2}$ on the gradient. For simplicity, we adopt the CVX MATLAB toolbox [33] for solving the TV model,

where we use $A=S_{n}F$ to be consistent with our setting (2). Note that the TV model (9) is parameter free, while we need to tune an algorithmic parameter for $L_{1}/L_{2}$. Please refer to Section 4 for more details on the $L_{1}/L_{2}$ minimization, in which one subproblem can be solved by CVX, and Section 6.1 for sensitivity analysis on this parameter.

By varying the value of $s$ that changes MS of the spikes in gradient, we compute the relative errors between the reconstructed solutions and the ground-truth signals. If we define an exact recovery for its relative error smaller than $10^{-6}$, we observe in Figure 1 that the exact recovery by $L_{1}$ occurs at $s\in[13,37],$ which implies that MS is larger than or equal to $26$. This phenomenon suggests that Theorem 1 might hold for sparse gradients by replacing the $L_{1}$ norm with the total variation. Figure 1 also shows the exact recovery by $L_{1}/L_{2}$ at $s\in[12,38],$ meaning that $L_{1}/L_{2}$ can deal with less separated spikes than $L_{1}$. Moreover, we further study the reconstruction results at $s=39,$ where both models fail to find the true sparse signal. The restored solutions by these two models as well as the different plots between restored and ground truth are displayed in Figure 2, showing that our ratio model has smaller relative errors than $L_{1}$.

Figure 2 illustrates that the TV solution can not reach the top of the bar in the ground-truth, which is referred to as loss-of-contrast. Motivated by this well-known drawback of TV, we postulate that the signal contrast may affect the performance of $L_{1}$ and $L_{1}/L_{2}$. To verify, we examine a two-bar step function, in which the contrast varies by the relative heights of the two bars. Following MATLAB’s notation, we set $\mathbf{u}(s+1:2s)=2,\mathbf{u}(\text{end}-2s+1:\text{end})=1,$ and the value of remaining elements uniformly as $t$; see Figure 3 for a general setting. We fix $s=12$, and vary the value of $t\in(1,2)$ to generate signals with different intensity contrasts. Considering four spikes in the gradient, we set $f_{c}=4$ or equivalently 9 low-frequency measurements to reconstruct the signal. The reconstruction errors are plotted in Figure 3, which shows that $L_{1}$ fails in all the cases, and $L_{1}/L_{2}$ can find the signals except for $t\in[1.5,1.65]$. We further examine a particular case of $t=1.65$ in Figure 4, where both models fail to get an exact recovery, but $L_{1}/L_{2}$ yields smaller oscillations than $L_{1}$ near the edges. Figures 3 and 4 demonstrate that $L_{1}/L_{2}$ is better at preserving image contrast than $L_{1}$.

We verify that all the solutions of $L_{1}$ and $L_{1}/L_{2}$ satisfy the linear constraint $A\mathbf{u}=\mathbf{b}$ with high accuracy thanks for CVX. We further investigate when the $L_{1}$ approach fails, and discover that it yields a solution that has a smaller $L_{1}$ norm compared to the $L_{1}$ norm of the ground-truth, which implies that $L_{1}$ is not sufficient to enforce gradient sparse. On the other hand, $L_{1}/L_{2}$ solutions often have higher objective value than the ground-truth, which calls for a better algorithm that can potentially find the ground-truth. We also want to point out that $L_{1}/L_{2}$ solutions depend on initial conditions. In Figure 1 and Figure 3, we present the smallest relative errors among 10 random vectors for initial values.

The minimum separation distance in 1D (Definition 1) can be naturally extended to 2D. In fact, there are two types of minimum separation definitions in 2D: one uses the $L_{\infty}$ norm to measure the distance [32], while another definition is called Rayleigh regularity [34, 35]. The exact recovery for 2D sparse vectors was characterized in [35] with additional restriction of positive signals. Both distance definitions were empirically examined in [12] for point-source super-resolution. When extending to 2D sparse gradient, one can compute the gradient norm at each pixel, and separation distance can be defined as the distance between any two locations with non-zero gradient norm. To the best of our knowledge, there is no analysis on the exact recovery of sparse gradients, no matter whether it is in 1D or 2D. In Section 3, we devote some empirical evidences, showing a similar relationship between sparse gradient recovery and minimum separation as Theorem 1, which calls for a theoretical justification in the future. Once the extension from 1D sparse vectors to 1D sparse gradients is established, it is expected that the analysis can be applied to sparse gradients in 2D to facilitate theoretical analysis in imaging applications.

Starting from (2), we incorporate an additional box constraint in the model, i.e.,

The notation $\mathbf{u}\in[p,q]^{N}$ means that every element of $\mathbf{u}$ is bounded by $[p,q].$ The box constraint is reasonable for image processing applications [36, 37], since pixel values are usually bounded by $[0,1]$ or $[0,255]$. On the other hand, the box constraint is particularly helpful for the $L_{1}/L_{2}$ model to prevent its divergence [19].

We use the indicator function to rewrite (10) into the following equivalent form

where $\Pi_{S}(\mathbf{t})$ denotes the indicator function that forces $\mathbf{t}$ to belong to a feasible set $S$, i.e.,

The augmented Lagrangian function corresponding to (11) can be expressed as,

where $\mathbf{g}$ is a dual variable and $\rho$ is a positive parameter. Then ADMM iterates as follows,

The update for $\mathbf{h}$ is the same as in [18], which has a closed-form solution of

where $\mathbf{r}^{(k)}$ is a random vector with the $L_{2}$ norm being $\sqrt[3]{\frac{\|D\mathbf{u}^{(k+1)}\|_{1}}{\rho}}$ and $\tau^{(k)}=\frac{1}{3}+\frac{1}{3}(\xi^{(k)}+\frac{1}{\xi^{(k)}})$ for

We elaborate on the $\mathbf{u}$-subproblem in (14), which can be expressed by the box constraint, i.e.,

To solve for (16), we introduce two variables, $\mathbf{v}$ for the box constraint and $\mathbf{d}$ for the gradient, thus getting

The augmented Lagrangian function corresponding to (17) becomes

where $\mathbf{w},\mathbf{y},\mathbf{z}$ are dual variables and $\beta,\gamma,\lambda$ are positive parameters. Here we have $k$ in the superscript of $\mathcal{L}$ to indicate that it is the Lagrangian for the $\mathbf{u}$-subproblem in (14) at the $k$-th iteration. The ADMM framework to minimize (17) leads to

where the subscript $j$ represents the inner loop index, as opposed to the superscript $k$ for outer iterations in (14). By taking derivative of $\mathcal{L}^{(k)}$ with respect to $\mathbf{u}$, we obtain a closed-form solution given by

where $I$ stands for the identity matrix. When the matrix $A$ involves frequency measurements, e.g. in super-resolution and MRI reconstruction, the update in (20) can be implemented efficiently by the fast Fourier transform (FFT) for periodic boundary conditions when defining the derivative operator $D$ in (1). For a general system matrix $A$, we adopt the conjugate gradient descent iterations [38] to solve for (20).

The $\mathbf{d}$-subproblem in (19) also has a closed-form solution, i.e.,

where $\mathbf{shrink}(\mathbf{x},\mu)=\mathrm{sign}(\mathbf{x})\odot\max\left\{|\mathbf{x}|-\mu,0\right\}$ is called soft shrinkage and $\odot$ denotes element-wise multiplication. We update $\mathbf{v}$ by a projection onto the $[p,q]$-box constraint, which is given by

In summary, we present an ADMM-based algorithm to minimize the $L_{1}/L_{2}$ on the gradient subject to a linear system with the box constraint in Algorithm 1. If we only run one iteration of the $\mathbf{u}$-subproblem (19), the overall ADMM iteration (14) is equivalent to the previous approach [18].

We intend to establish the convergence of Algorithm 1 with the box constraint, which is extensively tested in the experiments. Since our ADMM framework (14) share the same structure with the unconstrained formulation, we adapt some analysis in [21] to prove the subsequential convergence for the proposed model (10). For example, we make the same assumptions as in [21],

Unlike the unconstrained model (4), the strong convexity of $\mathcal{L}(\mathbf{u},\mathbf{h},\mathbf{g})$ with respect to $\mathbf{u}$ does not hold due to the indicator function $\Pi_{A\mathbf{u}=\mathbf{b}}(\mathbf{u})$. Besides, we have additional dual variable $\mathbf{w}$ which is not in the unconstrained model. To avoid redundancies to [21], we focus on the different strategies to the unconstrained case, such as optimality conditions and subgradient of the indicator function, when proving convergence for the constrained problem, e.g., in Lemmas 1-2 and Theorem 2.