Nonmonotone variable metric Barzilai-Borwein method for composite minimization problem

1.   Introduction

In this paper, we consider the following composite optimization problem:

where $f(x)$ is a differentiable function with a Lipschitz constant $L > 0$ for the gradient and $h(x)$ is a convex function that is possibly nonsmooth. Many models in sparse optimization can be formulated as problems (1.1), such as the $l_1$-regularized minimization in compressive sensing ^[1,2,3] and total variation regularization in image processing ^[4,5,6].

A special case of model (1.1) is the constrained problem over a convex set, when $h(x)$ is the indicator function of a convex set $\Omega$, i.e.,

and

For this case, a simple and famous algorithm is the gradient projection method ^[7]. To accelerate this projection algorithm, a scheme with line-search along the feasible direction is proposed in ^[8]. Using the Barzilai-Borwein (BB) stepsize ^[9] and nonmonotone line search strategy, Birgin et al. ^[10] developed the nonmonotone spectral projected gradient method, which is particularly appealing for large-scale problems. Dai and Fletcher ^[11] induced a new projected gradient method that alternately used the two BB steplengths. In ^[12], the authors proposed a modified spectral conjugate gradient projection algorithm and applied it to the total variation image restoration.

The other special case of model (1.1), which has attracted much interest in signal and image processing and machine learning, is the well-known $l_2$-$l_1$ problem, i.e.,

One of the most popular methods is the iterative shrinkage thresholding algorithm ^[13]. Using Nesterov's technique, Beck and Teboulle ^[14] proposed the famous fast iterative shrinkage thresholding algorithm (FISTA) by forward-backward splitting. Hale et al. ^[2] developed the fixed-point continuation method for $l_1$-minimization. They also used BB stepsize and nonmonotone line search to enhance this algorithm in ^[15]. Another closely related method is the sparse reconstruction by separable approximation algorithm (SpaRSA) ^[16], which involves minimizing a non-smooth convex problem with separable structures. Hager et al. ^[17] showed that SpaRSA has a sublinear convergence rate for general convex functions or R-linear when the objective function is strongly convex. An improved version of SpaRSA based on a cyclic BB iteration and adaptive line search was also presented in ^[17]. Huang and Liu ^[18] proposed a BB-type method for minimizing composite function. A nonmonotone strategy was employed for determining the step length along the search direction, which was defined by the difference between the minimizer of the approximation objective function and the current iteration point. Recently, Cheng et al. ^[19] induced an algorithm framework for the more general type of (1.1). At each step, the search direction was obtained by solving an optimization problem. Xiao et al. ^[20] proposed a nonmonotone BB gradient algorithm for $l_1$-regularized nonsmooth minimization in compressive sensing. At each iteration, the search direction can be easily derived by minimizing a local approximal quadratic model. Tseng and Yun ^[21] gave a block coordinate gradient descent method. But numerical results showed the efficiency of this method only for small and medium-scale problems.

There are many variable metric forward-backward (FB) approaches for solving (1.1) in recent literatures. In this type of algorithm, the suitable symmetric positive definite scaling matrices changed at each iteration. A clever choice of them can lead to an improvement in convergence speed ^[22,23]. Bonettini et al. ^[24] developed a variable metric inexact-line-search-based method for nonsmooth optimization. An Armijo-like rule was used to determine the stepsize along the descent direction. For the general nonconvex case, they proved that all limit points of the iterate sequence are stationary. In the last two decades, more variants of FB approaches have been proposed. These FB variants introduce a kind of extrapolation, or inertial, step exploiting the two previous iterations and can be traced back to two main typologies: FISTA-like methods ^[25] and heavy-ball ones ^[26]. Other recent results can be found in ^[27,28], where the analysis applies also to the nonconvex case, assuming that the Kurdyka-Lojasiewicz property holds.

In this paper, following the same lines as in ^[24] that exploiting some variable metric information and nonmonotone line search in ^[18,19] can lead to fast image restorations, we propose a BB-type algorithm with adaptive nonmonotone line search and scaling technique. The nonmonotone strategy includes a convex combination of the maximum function value of some previous iterations and the current function value. At each step, the search direction of the algorithm is obtained by solving an optimization subproblem involving a quadratic term with a variable metric matrix and BB steplength plus $h(x)$. Then, we perform a nonmonotone line search along this direction. We prove that the method with the nomonotone line search techniques is globally convergent. Similar to ^[17,18,24], the new method does not need to know the value of the Lipschitz constant $L$. Then, we apply this method to total variation image restoration, such as parallel magnetic resonance imaging ^[29,30], Poisson and Cauchy noise removal ^[24]. Numerical experiments show that our approach is competitive with several other known methods.

The rest of the paper is organized as follows: In Section 2, we present the new method for solving the model (1.1). The convergence of the new algorithm is proved in Section 3. In Section 4, we show the results of the numerical experiments. Finally, the conclusion is given in Section 5.

Throughout this paper, we denote $\langle x, y\rangle = x^Ty$ as the inner product of two vectors $x, y\in R^n$. $\|\cdot\|$ denotes the Euclidean norm. The scaled Euclidean norm associated with a symmetric positive definite matrix $D$ is

Given

$M_{[\mu_{\min}, \mu_{\max}]}$ denotes the set of all symmetric positive definite matrix with all eigenvalues contained in the interval $[\mu_{\min}, \mu_{\max}]$. For any

we have the result

2.   Proposed method

In this section, we first introduce the search direction and then introduce a new nonmonotone line search to determine the step length along this direction. At last, the new method is presented.

Our iteration step for solving (1.1) can be described as follows:

where $\lambda_k\in (0, 1]$ is the step length and $d_k\in R^n$ is the search direction, which is given by

where $\alpha_k > 0$, and $D_k$ is a symmetric positive definite matrix. Let

from (2.2), we can get that $\tilde{x}_k$ is the minimizer of

In fact, the problem (2.3) can be rewritten as the proximity operator style associated with the convex function $h(x)$, which is defined as

where

So, (2.1) is equivalent to

The setting of $D_k$ and $\alpha_k$ will affect the performance of the method (2.5). A clever choice of them can lead to significant improvements in the practical convergence speed ^[31,32]. Firstly, about the metric selection, we assume that

Different problems can choose a different scaling matrix $D_k$. So, we will discuss how to choose $D_k$ in the section on numerical experiments. Next, we will focus on how to choose the parameter $\alpha_k$. Similar to ^{[17,18,19,20]}, we choose $\alpha_k$ as follows:

where

In fact, this step length was first proposed by Barzilai and Borwein for solving unconstrained problems ^[9]. Due to its efficiency, the BB approach has received considerable attention. As said in ^[9], a good property of $\alpha_k$ is that it satisfies the quasi-Newton condition

We take

where $\alpha_{\max} > \alpha_{\min} > 0$ are fixed constants.

When the search direction is determined, a suitable stepsize along the decent direction should be found to determine the next iterative point. In the following, we will show how to determine the step length $\lambda_k$. We know that the performance of the BB algorithm benefits from the nonmonotone line search (see ^[10,11]). The earliest nonmonotone line search technique, which differs from the traditional Armijo line search or the Wolfe-Powell line search, was given by Grippo et al. ^[33] for smooth, unconstrained optimization. In order to overcome some disadvantages of the nonmonotone line search, Amini et al. introduced a more relaxed nonmonotone strategy that included a convex combination of the largest objective function in some recent past iteration and the current function value ^[34]. However, they are not directly applied to the nonsmooth problems. In ^[18], the authors proposed a line search for solving the nonsmooth problem (1.1) as follows:

where $\gamma\in (0, 1)$ is a constant and

where $M$ is a fixed integer and $l(k)$ is an integer such that

Motivated by the idea from ^[18,34], we modified the line search (2.8) and used the following acceptance criterion to determine the step length $\lambda_k$:

where

and $0\leq\eta\leq1$. Since

we have

Given all the above derivations, we now describe the nonmonotone variable metric Barzilai-Borwein type algorithm as Algorithm 1 shows.

Algorithm 1 is closely related to the method in ^[18]. However, compared (2.10) with (2.8), it can be easily found differences. Algorithm 1 use the adaptive nonmonotone line search and the variable metric matrix $D_k$ to get the step length $\lambda_k$. If $\eta = 1$ and $D_k = I$, Algorithm 1 is reduced to the algorithms in ^[18]. Similar to the idea in ^[24], the variable metric technique is adopted to compute the descent direction, while the paper ^[24] used the modified Armijo line search to determine the step length $\lambda_k$.

3.   Convergence analysis

In this section, we analyze the convergence of Algorithm 1.

First, we give some useful properties of $d_k$. In the nonsmooth case (1.1), $d_k$ is a descent direction for $\phi(x)$ at $x_k$ if $\phi'(x_k; d_k) < 0$, where $\phi'(x; d)$ is the direction derivative of $\phi(x)$ at $x$ in the direction $d\in R^n$ and is defined as ^[35]

A feasible point $x$ is said to be a stationary point of $\phi(x)$ if $\phi'(x; d)\geq0$ for all $d$. Using the definition of $\phi'(x; d)$ and the property of stationary point, we show that the direction defined by (2.2) is descent if $d_k\neq0$ and if $d_k = 0$, $x_k$ is a stationary point of $\phi(x)$ in the following two lemmas. The proof is similar to the proofs of Lemmas 2.2 and 2.3 in ^[19]. So, we omit it here.

Lemma 3.1. For any $x_k\in R^n$ and $d_k\in R^n$ determined by (2.2), we have

Lemma 3.2. $x_k\in R^n$ is a stationary point of (1.1) if and only if $d_k = 0$.

Next, we show that Algorithm 1 is well defined.

Lemma 3.3. Let $\{x_k\}$ be a sequence generated by Algorithm 1. If $x_k$ is not a stationary point, then the line search (2.10) is satisfied whenever

Moreover, if $\lambda_k$ satisfies the inequality (2.10), then

where $\rho\in (0, 1)$.

Proof. By the inequality (2.12), the Lipschitz continuity of $\nabla f(x)$ with constant $L$, the convexity of $h(x)$, and Lemma 3.1, we can get

The use of inequality (2.10) yields

Using the assumption (1.2) and Lemma 3.2, we can have

Now, we give the lower bound of $\lambda_k$. Assume that either $\lambda_k = 1$ or the line search (2.10) has failed at least once. Therefore, there is a constant a $\rho\in(0, 1)$ such that

Combining the inequality (3.2) yields

So,

This completes the proof. □

The following theorem implies that every accumulation point of $\{x_k\}$ is a stationary point of (1.1). Moreover, if $f(x)$ is a convex function, Algorithm 1 converges to a global solution of (1.1).

Theorem 3.1. Let the sequences $\{x_k\}$ and $\{d_k\}$ be generated by Algorithm 1, then

and

where $\bar{\phi}$ is a constant.

Proof. Using the definition $R_k$, the inequalities (2.10) and (2.12), we have

It indicates that

On the other hand, from (2.9), we get

Then the sequence $\{\phi(x_{l(k)})\}$ is nonincreasing. Using the same assumption as the one in ^[19] that $\phi(x)$ is bounded below, there is a constant $\bar{\phi}$ such that

By applying (3.4) with $k$ replaced by $l(k)-1$, we have

Using (3.6), we obtain

By the continuity of $\phi(x)$, we deduce

Let $k = l(k)-j$, by induction and using the similar proof of Theorem 3.1 in ^[19], we can show that the following limits are satisfied for $j = 1, 2, \cdots, M$,

For

we have

This means that

The facts with the continuity of $\phi(x)$ and (3.6) indicate that

Taking the limits of (3.4) and (3.7), we obtain

From Lemma 3.2, we can also know that any limit point of $\{x_k\}$ is a stationary point of $\phi(x)$. Therefore, the proof is complete. □

4.   Numerical experiments

In this section, we give some numerical experiments to illustrate the effectiveness of the proposed method. The first test problem of the form (1.1) is the total variation (TV) of magnetic resonance imaging (MRI). As a special type of (1.1), the parallel magnetic resonance imaging model ^[36,37] can be written as follows:

where $x\in C^n$ and $n$ is the total number of pixels in the image. $A\in C^{m\times n}$ is a sample matrix, and $b\in C^m$ is the vector of measured Fourier coefficients. $\|\cdot\|_{TV}$ is the total variation norm, which was first introduced in ^[5]. $\tau > 0$ is the regularization parameter. Obviously,

is Lipschitz continuous,

is nonsmooth, and $\phi(x)$ is a convex function.

Our new method is named NMBBA, NVMBB, and NVMBBA, respectively, in different parameter settings. We compare it with the famous algorithms SpaRSA ^[16] and NMBB ^[18]. In all the experiments, we set

In the first experiment, we use a Poisson mask with an acquisition rate 0.243. The image size is $256\times256$ (named Data1). The regularization parameter is set as $\tau = 0.0005$. In NMBBA, $\eta = 0.5$, $D_k = I$. In NVMBB,

In NVMBBA,

In the second experiment, we use radial mask with an acquisition rate 0.424. The image size is $512\times512$ (named Data2). The regularization parameter is set as $\tau = 0.00037$. In NMBBA, $\eta = 0.5$, $D_k = I$. In NVMBB,

In NVMBBA,

For the subproblem (2.4), we use the FISTA ^[14] method to solve the minimizing total variation proximal function. The inner iteration number is set to 5. The stopping criterion of the five algorithms is:

We use signal-to-noise ratio (SNR) as a means of judging performance, which is defined as follows:

where $x_0$ is the original image, $x_k$ is the recovery, and $\hat{x}$ is its mean value. All algorithms are implemented in the Matlab programming environment (Version R2011a). The experiments are performed on a laptop with 2 GB of RAM, a 2 GHz processor, and the Windows 7 operating system.

The reconstructed images are shown in Figures 1 and 2. It is seen that all algorithms adequately recovered the image from a small set of data samples. We can also see that NMBB appears to have higher accuracy than other methods for Data1 and NMBBA foe Data2.

This can also be seen in Table 1. From the table, we can see that the proposed method is very efficient in the sense that it reaches a similar SNR value in a shorter time (second) and with a lower iteration number (Iter) compared to the other methods.

In Figure 3, we plot the objective function value versus CPU time. From the figure, we can see that, among these five algorithms, the decay rates of SpaRSA are much slower than those of the others.

The reason may be that it has to solve the subproblem for each trial point on each line search, as said in ^[18]. NMBBA does not use the variable metric technique. This demonstrates that the convex combination of the maximum function value of some previous iteration and the current function value can improve the convergence speed. If we choose the suitable scaling matrix $D_k$ for NVMBB and NVMBBA, they may converge quickly; for example, see their performance Data1. Otherwise, NVMBB and NVMBBA may be a little slower than NMBBA but also faster than SpaRSA and NMBB. So, a suitable combination of the stepsize $\alpha_k$ and the variable metric $D_k$ is important for the practical problems discussed in ^[31].

The second test problem is TV Poisson noise image deconvolution. When the image is corrupted by Poisson- type noise, the Kullback-Leibler divergence is often used as the data-fitting term. Let $b$ be the observed blurry and noisy image, then the data discrepancy can be given as

where

is a blur matrix, $\textbf{1}$ is the vector of all ones, and $g$ is a constant background term. To deal with this ill-posed problem, based on the prior information of the image, TV regularization with a nonnegative constraint is added to control the noise and artifacts, i.e.,

where $\rho > 0$ is a regularization parameter. Hence, the Poisson noise image restoration problem is reformulated as a minimization problem

So, Algorithm 1 can be used to solve the convex composite optimization problem.

In this test, we compare it with the method named VMILA in ^[24]. In Algorithm 1, the variable metric technique is the same as VMILA, with $\eta = 1$. This is to say, we use the nonmonotone line search instead of the modified Armijo line search in VMILA to determine the step length $\lambda_k$. So, algorithm 1 is named NONVMILA. The two test images are Micro (128$\times$128) and Phantom (256$\times$256). For the micro image, $\rho = 0.09$ and $g = 0.5$. For the Phantom image, $\rho = 0.004$ and $g = 10$. We set $M = 3$ for Micro and $M = 9$ for Phantom, and the other parameters are the same as VMILA. We also used the FISTA ^[14] method to minimize the total variation of the proximal function. The inner iteration number is set to 1500. We refer the reader to ^[24] for more details. Both of the compared methods are stopped when the maximum number of iterations is more than 400. The peak-signal noise ratio (PSNR) is used to evaluate the quality of the restored image.

The numerical performance of image restoration is recorded in Table 2. We present the CPU times in seconds and the PSNR values of the two algorithms for various blurry observations.

From the table, one can see that the two approaches can produce similar behaviors for the two images. Figure 4 displays the visual results of image deblurring. In terms of visual quality, the difference is not significant. This can also be seen from the PSNR values. Figure 5 reports the relative decrease of the objective function with respect to the optimal value $\phi^*$ (when $x = x_0$, $x_0$ is the original image) as a function of the computational time. We observe that NONVMILA requires less time.

The last problem is Cauchy noise image deconvolution. When the image is corrupted by the Cauchy-type noise, the log function is used as the data fitting term. Let $b$ be the observed blurred and noisy image, then the data discrepancy is

where

is a blur matrix with the 9$\times$9 testing kernel and the variance 1, while has been set equal to 0.02. Based on the prior information of the image, TV regularization with a nonnegative constraint is added to control the noise and artifacts, i.e.,

where $\rho > 0$ is a regularization parameter. Hence, the Cauchy noise image restoration problem is reformulated as

So, Algorithm 1 can be used to solve this composite optimization problem.

In this test, we also compare it with the VMILA type in ^[38,39]. In Algorithm 1, the variable metric technique is the same as in VMILA, with $\eta = 0.9$, and is also named NONVMILA. We set $M = 5$ and the other parameters are the same as VMILA. The three test images are the cameraman (256$\times$256), the house (256$\times$256) and the bird (256$\times$256). All methods are stopped when the maximum number of iterations is more than 5000 or

The numerical performance of image restoration is reported in Table 3.

We present the running time in seconds and PSNR values of the two algorithms. From the table, we observe that the two approaches can get similar behaviors for three images in terms of PSNR. We also find that NONVMILA requires less time than VMILA when the stopping criterion is satisfied. Figure 6 displays the visual results of the deblurred image. From Table 3 and Figure 6, it seems that our approach is very efficient for this tested problem.

5.   Conclusions

In this paper, we introduce a variable metric algorithm frame for composite optimization problems. With the nonmonotone line search techniques, the convergence of the method is discussed. With the growing use of nonconvex objective functions in applied fields like image processing and machine learning, the need for numerical methods in a fully nonconvex setting has increased significantly. So, future work will focus on a general nonconvex non-Lipschitz function, including theoretical and numerical analysis. Moreover, there is good potential to improve the convergence results using the conditions about the regularity of the cost function and the boundedness of the Hessian matrix ^[40]. Also, a multi-step inertial forward-backward splitting algorithm will be considered in the future ^[41,42]. At last, it is interesting to develop a new scaling matrix and BB stepsize for solving the composite optimization problem.

Author contributions

Xiao Guo: writing-original draft, writing-review and editing, software. Chuanpei Xu: resources, writing-review and editing, supervision. Zhibin Zhu: methodology, writing-review and editing, supervision. Benxin Zhang: resources, writing-review and editing.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

The authors wish to thank the editor and anonymous referees for their constructive suggestions and comments, which greatly improved the presentation of this paper. This work was supported in part by the National Natural Science Foundation of China under Grants 11901137 and 62171147, in part by the Guangxi Natural Science Foundation under Grant 2024GXNSFAA010512, in part by the Guangxi Key Laboratory of Automatic Detecting Technology and Instruments under Grant YQ22108, and in part by the Innovation Project of GUET Graduate Education under Grant 2022YCXS161. We would like to thank the authors of ^[24,36,37] for making their code and data free for academic use.

Conflict of interest

The authors declare no conflicts of interest.