Global optimality analysis and solution of the $\ell_0$ total variation signal denoising model

1.   Introduction

The total variation (TV) regularization term has been used broadly in many areas, such as statistics ^[1,2,3], signal and image processing ^[4,5,6,7], in order to get piecewise constant estimator. The TV-based signal denoising methods are very effective for recovering piecewise constant (PWC) signals compared with conventional linear time-invariant filtering approach because the TV regularizer can capture the jump-sparsity in data. The classical TV (L1-TV) denoising model ^[4] based on $\ell_1$ norm penalty function is formulated as a strongly convex optimization problem. For the observed signal $y\in R^n$, the L1-TV denoising model is defined as

where $\lambda > 0$ is a regularization parameter, and $D\in R^{(n-1)\times n}$ is the first-order difference operator

In this case, the optimization problem (1.1) has unique global optimal solution, and it can be solved in finite time by using very fast direct (noniterative) algorithms ^[8,9]. However, it has been shown in the literatures that nonconvex penalties can lead to more accurate estimation by comparison to $\ell_1$ penalty, see, e.g., ^{[10,11,12,13]}. To enhance the performance of L1-TV denoising method, numerous nonconvex penalties function are used to substitute the $\ell_1$ norm, see, e.g., ^{[5,14,15,16,17,18,19]}. In view of the sparsity of the derivative of the underlying signal, the number of discontinuous points is a natural and cogent regularization term ^[17,18,19], namely,

However, it is a nonconvex and discontinuous function. To escape the computational challenge arising from this regularizer, some continuous surrogates of $\|Dx\|_0$ are proposed in ^[5,14,15,16]. In ^[5], the authors adopt the logarithmic penalty and arctangent penalty to substitute $\ell_1$ norm in (1.1), and show that the corresponding objective functions are both convex if regularization parameter $\lambda$ is less than a threshold. In ^[14], the authors propose the MCP-TV denoising method by using minimax concave penalty ^[12] to induce the sparsity of $Dx$. Based on the generalized Moreau envelope, Selesnick et al. ^[16] define the generalized Moreau envelope TV (GME-TV) penalty with matrix parameter $B$ and use in denoising of signals. Note that the L1-TV, MCP-TV and Moreau enhanced TV ^[15] are the special cases of the GME-TV regularizer.

In this paper, we revisit the $\ell_0$ TV regularized signal denoising model (L0TV), which is also called L2-Potts ^[17,18]. The corresponding optimization model is

where $\|Dx\|_0$ is the $\ell_0$ pseudo-norm of $Dx$, for counting the number of non-zero elements of $Dx$. The objective function of (1.2) is discontinuous and nonconvex. In ^[17,18], the model (1.2) is solved by dynamic programming algorithm. To well resolve the challengeable problem, we first equivalently rewrite the L0TV as a partial $\ell_0$ regularized (PL0R) optimization problem by resorting the interesting structure of $D$. For the PL0R problem, we define a P-stationary point and show that it is a global optimal solution. These fascinating theoretical results allow us to find a global optimal solution of the model (1.2). Although many existing numerical methods can be used to solve the $\ell_0$ regularized (L0R) problem such as iterative hard-thresholding algorithm ^[20,21], penalty decomposition ^[22], active set Barzilar-Borwein ^[23], these methods cannot be used to solve the model (1.2) and may suffer from low accuracy and slow convergence due to only make use of the first-order information of the involved functions. Very recently, an efficient Newton-type algorithm was proposed for L0R problem, and its global and quadratic convergence was established ^[24]. Inspired by this work, we design an effective algorithm based on the Newton's method to pursuit the global optimal solution of the L0TV model (NL0TV). The algorithm has a low computational complexity since a small-scale linear equation system is solved to update the Newton direction in per iteration. Comparing with several state-of-the-art methods, our experimental results show that the NL0TV achieves an excellent performance. The main contributions of this paper can be summarized as follows.

● An equivalent reformulation of the L0TV denoising model (1.2) is proposed in Theorem 1. After that, its necessary and sufficient conditions of global optimal solution are established in Theorem 2.

● We systematically analyze the relationships among the L0TV model (1.2), PL0R model (2.1) and the optimization problem (2.3) in terms of global optimal solution.

● Based on these theoretical results, an effective Newton algorithm is designed to pursuit the global optimal solution of the model (1.2).

● The excellent performance of our approach is illustrated by comparing with several state-of-the-art methods.

The remainder of this paper is organized as follows. Section 2 includes an equivalent reformulation of the L0TV denoising model, optimality conditions and Newton algorithm. Experiments and conclusions are presented in Sections 3 and 4, respectively.

2.   Our method

The model (1.2) is a natural and cogent model for denoising of signals. But it is generally difficult to find the global optimal solution due to the intrinsic combinatorial property and inseparability of the regularization term. Next, we will give an equivalent reformulation of (1.2), which is helpful for designing effective optimization algorithm.

For convenience, we first definite several symbols. Let $z = Dx\in R^{n-1}$ and

Based on these symbols, we give an equivalent reformulation of (1.2) in next theorem.

Theorem 1. The L0TV signal denoising model (1.2) is equivalent to the following optimization problem

where $G = [e, H]\in R^{n\times n}$ and $e$ is an $n$-dimensional column vector with each component is one.

Proof. It follows from $z = Dx\in R^{n-1}$ that $z_i = x_{i+1}-x_{i}$, for any $i = 1, 2, \cdots, n-1$, by means of the fascinating structure of $D$. Thus,

Together with $w = (x_1; z)$, we have

and $w = G^{-1}x$. Hence, the L0TV signal denoising model (1.2) can be equivalently written as (2.1) by substituting $z$ and $w$ into (1.2).

The model (2.1) is named as partial $\ell_0$ regularized (PL0R) optimization problem because the regularization term in (2.1) is with respect to $z$ rather than $w$. Note that its the first term $f(w) = \frac{1}{2}\|Gw - y\|_2^2$ is a smooth and strongly convex function. Specifically, the gradient and Hessian of $f(w)$ are

and $\nabla^2f(w) = G^\top G$, respectively. Moreover, $\nabla^2f(w)$ is a positive definite matrix due to the nonsingularity of $G$. Hence, there exist two positive constants $L_f$ and $l_f$ such that all eigenvalues of $\nabla^2f(w)$ are greater than $l_f$ and less than $L_f$.

2.1.   Optimality conditions

Some first-order optimality conditions of $\ell_0$ regularized problem have been established in Theorems 2.2 and 2.4 of ^[22]. But it is worth noting that they analyzed merely the relationship between Karush-Kuhn-Tucker (KKT) type stationary points and local minimizer. Inspired by the definition of the L-stationarity in [25,Definition 4.8], we introduce a P-stationary point of the PL0R optimization problem. For convenience, we first briefly review the definition of L-stationarity. Let $L > 0$. A vector $x^\ast$ is called an L-stationary point of the $\ell_0$ regularized optimization problem $\min_{x} g(x) + \lambda\|x\|_0$ if

Here $\text{Prox}_{\frac{\lambda}{L}\|\cdot\|_0}(\cdot)$ is the proximal operator of $\frac{\lambda}{L}\|\cdot\|_0$, and defined as

Next we present the definition of P-stationary point for the PL0R optimization problem (2.1).

Definition 1. Let $\alpha > 0$. We say that $w^\ast = (x_1^\ast; z^\ast)$ is a P-stationary point with $\alpha$ of (2.1) if

It is worth mentioning that the proximal operator of $\alpha\lambda\|\cdot\|_0$ is multi-valued due to the non-convexity of $\alpha\lambda\|\cdot\|_0$. Fortunately, the analytic formula of $\text{Prox}_{\alpha\lambda\|\cdot\|_0}(\cdot)$ can be characterized by using the intrinsic discreteness and separability. Specifically, for any $i = 1, 2, \cdots, n-1$,

where $\mathcal{J}: R\rightarrow \{0, 1\}$ is the function defined by $\mathcal{J}(0) = 0$ and $\mathcal{J}(v_i) = 1$ for $v_i\neq0$. It is easy to show that the minimum value of $\mathcal{K}(v_i)$ is the smaller one of $\alpha\lambda$ and $\frac{1}{2}u_i^2$. If $\alpha\lambda > \frac{1}{2}u_i^2$, then the minimizer of $\mathcal{K}(v_i)$ is $0$. If $\alpha\lambda = \frac{1}{2}u_i^2$, then the minimizer of $\mathcal{K}(v_i)$ is $0$ or $u_i$. If $\alpha\lambda < \frac{1}{2}u_i^2$, then the minimum value and minimizer of $\mathcal{K}(v_i)$ are $\alpha\lambda$ and $u_i$ respectively. Hence, the proximal operator of $\alpha\lambda\|\cdot\|_0$ can be characterized as follows:

Whereafter, the relationship between P-stationary point and global optimal solution of (2.1) is revealed.

Theorem 2. For the PL0R optimization problem (2.1), the following results hold.

(a) (Sufficiency) Let $(x_1^\ast; z^\ast)$ be a P-stationary point with $\alpha\geq1/l_f$, then it is a global optimal solution.

(b) (Necessity) If $(x_1^\ast; z^\ast)$ is a global optimal solution, then it is a P-stationary point with $\alpha\in(0, 1/L_f)$.

Proof. (a) Since $(x_1^\ast; z^\ast)$ is a P-stationary point. It follows from Definition 1 that $0 = \nabla_{x_1}f(w^\ast)$ and

According to the optimality of $z^\ast$, we have

After some simple manipulations, we can further obtain that

From the strong convexity of $f$, one can see, for any $w\in R^n$,

Combining the above two aspects, we have

Here, the last inequality of the above equation holds from $\alpha\geq1/l_f$. Hence, $(x_1^\ast; z^\ast)$ is a global solution of (2.1).

(b) Suppose $w^\ast = (x_1^\ast; z^\ast)$ is a global optimal solution of (2.1). Then, from Fermat's rule [26,Theorem 10.1], we have

which imply the first equation in (2.2) holds. On the other hand, from the strong smoothness of $f$, we obtain that

Let $x_1 = x_1^\ast$ and $z\in\Omega = \text{Prox}_{\alpha\lambda\|\cdot\|_0}(z^\ast - \alpha\nabla_{z}f(w^\ast))$. We know

where last inequality holds from the fact that $(x_1^\ast; z^\ast)$ is a global solution of (2.1). This together with $\alpha < 1/L_f$ leads to

yielding $z = z^\ast$. Therefore,

Namely, $(x_1^\ast; z^\ast)$ is a P-stationary point with $\alpha\in(0, 1/L_f)$.

To solve the PL0R problem well, we first consider an optimization problem that is closely relevant to it, namely,

Here, $h(z) = \frac{1}{2}(Hz-y)^\top M(Hz-y)$, $M = I_n-ee^\top/n$ and $I_n$ is an identity matrix. Thus, the gradient and Hessian matrix of $h(z)$ are

respectively. Moreover, $h(z)$ is a strongly convex function because its Hessian matrix is positive definite. Hence, there exists a constant $l_h > 0$ such that, for any $z^1, z^2$,

Similarly, we say that $z^\ast$ is a P-stationary point with $\alpha > 0$ of (2.3) if

The following theoretical results present systematically the relationship between (2.1) and (2.3).

Theorem 3. A point $z^\ast$ is a P-stationary point with $\alpha$ of (2.3) if and only if

is a P-stationary point with $\alpha$ of the PL0R model (2.1).

Proof. It follows that

and $\text{Prox}_{\alpha\lambda\|\cdot\|_0}(z^\ast - \alpha\nabla_{z}f(v^\ast)) = \text{Prox}_{\alpha\lambda\|\cdot\|_0}(z^\ast - \alpha\nabla h(z^\ast))$. Together with the equivalence of the first equation in (2.2) and $x_1^\ast = e^\top(y-Hz^\ast)/n$, the desired conclusions are proved.

Theorem 4. Let $z^\ast$ be a global optimal solution of (2.3). Then it is a P-stationary point of (2.3) with $\alpha > 1/l_h$, and $v^\ast = (e^\top(y-Hz^\ast)/n, z^\ast)$ is also a global optimal solution of (2.1).

Proof. Since the strong convexity of $h(z)$, we obtain that

Let $z^\ast$ be a global optimal solution of (2.3) and $z$ be a P-stationary point of (2.3). Then we have

where last inequality holds from the fact that $z^\ast$ is a global optimal solution of (2.3). This together with $\alpha > 1/l_h$, we obtain that $z^\ast = z$. Therefore, $z^\ast$ is a P-stationary point of (2.3). From Theorem 3 and Theorem 2 (a), we obtain that the desired conclusion immediately.

2.2.   Newton algorithm

These optimality conditions establish a foundation for finding the global optimal solution of (1.2) effectively. Particularly, our algorithm can be divided into three steps, see Algorithm 1 for details.

After Step 1 and Step 2, we obtain $(x_1^*; z^*)$ which is a P-stationary point of (2.1) by Theorem 3. Further, it follows from Theorem 2 (a) that $(x_1^*; z^*)$ is a global optimal solution of the PL0R model (2.1). Therefore, we acquire the global optimal solution of the L0TV model (1.2) in Step 3. The main computational cost of NL0TV is to find a P-stationary point of (2.3). Inspired by ^[24], we adopt Newton-type algorithm to solve the model (2.3) and thereby obtaining its P-stationary point. To express the solution of (2.4) more explicitly, we introduce the following stationary equation

where $T = \left\{i: \vert z_i - \alpha\nabla_{i}h(z)\vert\geq\sqrt{2\lambda\alpha}\right\}$ and $\overline{T}$ is the complementary set of $T$. The relationship between stationary equation (2.5) and a P-stationary point of (2.3) is established in next theorem.

Theorem 5. Let $z^\ast$ be a solution of (2.5). Then it is a P-stationary point of (2.3).

Proof. Suppose that $z^\ast$ is a solution of (2.5). Then $T^\ast = \left\{i: \vert z_i^\ast -\alpha\nabla_{i}h(z^\ast)\vert\geq\sqrt{2\lambda\alpha}\right\}$ and

Moreover, $\vert z_i^\ast\vert\geq\sqrt{2\lambda\alpha}$ for any $i\in T^\ast$, and $\vert\nabla_{i}h(z^\ast)\vert < \sqrt{2\lambda/\alpha}$ for any $i\in\overline{T}^\ast$. Together with [20,Lemma 2], we obtain $z^\ast$ is a P-stationary point of (2.3).

To find a solution of (2.5), we first need to locate the index set $T$ that is unknown in general and then solve the equation. Therefore, we employ an adaptively updating rule as follows. Let $z^k$ be the $k$-th iteration point. We first calculate an approximation $T_k$, and then apply the Newton method to solve (2.5) with $T_k$. Namely, update $d^k$ and $z^{k+1}$ by, respectively,

and

Here $t_k$ is a step size generated by the Amijio line search as described in ^[24]. One can observe that the major computational cost of Newton algorithm arises from solving the equations (2.6). Its complexity is approximately $O(\lvert T_k\rvert^3)$ in the $k$-th iteration. Note that, a sparse solution $z^*$ is admitted, namely, $\|z^\ast\|_0\ll n-1$, then $\lvert T_k\rvert$ can be quite small. Hence, NL0TV has a low computational complexity. As shown in ^[24], the generated sequence of Newton algorithm converges to a P-stationary point of (2.3) globally and quadratically. Thus, the global optimal solution of (2.1) can be obtained according to Theorem 4.

3.   Experiments

This section conducts numerical experiments to illustrate the effectiveness of our proposed NL0TV method. All experiments are implemented by MATLAB (R2020a) on a personal laptop with 8 GB of RAM and Inter Core i7 2.3 GHz CPU.

To demonstrate the excellent performance of our approach, we compare NL0TV with five state-of-the-art TV-based denoising methods, including L1-TV ^[8], Atan-TV ^[5], MCP-TV ^[14], GME-TV ^[16] and L2-Potts ^[18]. For a fair comparison, the regularization parameters of all methods are traversed in $\{0.1, 0.2, \cdots, 5\}$ to output the best experimental results. The root-mean-square error (RMSE), mean-absolute-error (MAE) and signal-to-noise ratio (SNR) are used to quantify the performance of NL0TV and other compared methods. For comparison, we use the same test signal as described in the relevant literatures ^[5,15,16]. The true signal $x^\ast$ is named 'blocks', which is generated by the function "$\textsf{MakeSignal}$" in the Wavelab software library (see https://statweb.stanford.edu/ wavelab/ for details). The noisy signal $y$ is obtained by adding Gaussian noise to the true signal, i.e., $y = x^*+\sigma N$, where $N$ is the standard normal distribution, $\sigma$ is the noise factor. Figure 1 illustrates an example of the 'blocks' signal with noise ($\sigma = 0.7$). Experimental results of six methods on the noisy signals with different $\sigma$ are summarized in Figures 2-3 and Table 1.

Figure 2 reports the experimental results of six methods on the noisy signal with $\sigma = 0.7$. Some comments on Figure 2 can be made. (i) The RMSE of L1-TV, Atan-TV and MCP-TV methods are more than 0.3, while our approach achieves the minimum values 0.1460. (ii) Regarding MAE, NL0TV method obtains the minimum values compared with other five methods. Moreover, our approach also achieves the biggest SNR in comparison with other methods, which shows that the proposed method can remove more contaminants in the noisy signal. (iii) As presented in Figure 2, there are some small jumps in the signals estimated by L1-TV, Atan-TV, MCP-TV and GME-TV methods that are not present in the true signal, but this is rarely the case with NL0TV and L2-Potts methods. The reason may be that the $\ell_0$ TV regularization term is more suitable for piecewise constant data than its surrogates. (iv) It can be seen from Figure 2 (e) and (f) that the NL0TV possesses superior denoising performance comparing to the L2-Potts method, which indicates that our Newton algorithm can obtain more accurate solution than dynamic programming algorithm developed in ^[17,18,27]. Overall, the performance of NL0TV is better than other compared methods in terms of the RMSE, MAE and SNR.

To further evaluate the stability of the denoising performance of our approach, we repeat 100 times noise realizations for each $\sigma$ value and report the experimental results of six methods in Table 1 and Figure 3.

Table 1 presents the mean and standard deviation (SD) of RMSE, MAE and SNR for each method on the noisy signal with $\sigma = 0.6$. The values in bold represent the best performances. It is not difficult to see the following facts from this table. (i) The smallest values of average RMSE and average MAE are obtained by our approach. These experimental results imply that the proposed method possesses the best denoising performance. In particular, NL0TV method reduces the average RMSE by $48.49\%$ and $18.44\%$ compared with the L1-TV method and the L2-Potts method, respectively. (ii) The largest value of average SNR is achieved by the proposed approach. Specifically, NL0TV method increases the average SNR by $29.73\%$ compared with the L1-TV method. Moreover, our approach increases the average SNR by $6.51\%$ compared with the L2-Potts method, which indicates that the proposed Newton algorithm is effective. (iii) NL0TV method achieves the minimal standard deviation of RMSE and MAE. This phenomenon indicates that the proposed method is more stable compared with other several methods.

Figure 3 shows the experimental results of 100 replicates of each method on the noisy signal with $\sigma = 0.8$. This boxplot conveys several interesting phenomena. (i) In Figure 3 (a), the median (red line) of the RMSE of the NL0TV method is lower than that of the other five methods. (ii) In Figure 3 (a) and (b), the boxplot of our method does not have a red plus sign ({$+$}), which illustrates that our method does not produce abnormal experimental results. (iii) In Figure 3 (c), the median (red line) of the SNR of our method is higher than that of the other five methods. These phenomena show that the L0TV method is superior to other methods.

4.   Conclusions

In this paper, we first give an equivalent reformulation of the L0TV model and analyze the relationship between its P-stationary point and global optimal solution. Based on these theoretical results, an efficient approach based on Newton's method is designed to find the global optimal solution of the L0TV model. The algorithm has a considerably low computational complexity and enjoys an excellent performance when against five state-of-the-art methods. Furthermore, the techniques developed in this paper can be also used in conjunction with other TV-based methods in order to process more general signal and image processing problems.

Acknowledgment

We would like to thank the editors and anonymous referees for their insightful remarks and constructive comments that have helped to improve the quality of the paper. This work was supported in part by the Natural Science Foundation of Guangxi (2021GXNSFBA075012); in part by the Project funded by China Postdoctoral Science Foundation (2022M723327); in part by the Middle-aged and Young Teachers' Basic Ability Promotion Project of Guangxi(2020KY08005); and in part by the Doctoral Foundation of Guangxi University of Science and Technology (19Z40).

Conflict of interest

The authors declare there is no conflict of interest.