An efficient gradient-free projection algorithm for constrained nonlinear equations and image restoration

1.   Introduction

Recently, using Solodov and Svaiter's projection technique ^[1], several conjugate gradient methods for solving large-scale unconstrained optimization problems have been extended to solve nonlinear equations with convex constraints (see, ^{[2,3,4,5,6,7,8,9]} and the references therein). Due to its simplicity, low storage requirement, and applications, the method has been of interest to various research communities ^{[10,11,12,13,14]}. As known, the Fletcher-Reeves (FR) ^[15], Conjugate Descent (CD) ^[16] and Dai-Yuan (DY) ^[17] conjugate gradient methods have strong convergence properties, but due to jamming, they do not do well in practice. Having said that, the Hestenes-Stiefel (HS) ^[18], Polak-Ribiére-Polyak (PRP) ^[19,20], and Liu-Storey (LS) ^[21] conjugate gradient methods do not necessarily converge, but they often work better than FR, CD and DY. In ^[22], in order to combine the numerical efficiency of the LS method and the strong convergence of the FR method, Djordjević proposed a hybrid LS-FR conjugate gradient method for solving the unconstrained optimization problem. In her work, the conjugate gradient parameter was computed as a convex combination of the LS and FR conjugate gradient parameter. The hybridization parameter for the convex combination was obtained in such a way that the direction of the proposed method satisfies the condition of the Newton direction but also at the same time, it satisfies the famous Dai-Liao conjugacy condition.

In an attempt to extend the LS-FR method of Djordjević to solve monotone nonlinear equations with convex constraints, Ibrahim et al. ^[23] proposed a derivative-free hybrid LS-FR conjugate gradient method with a conjugate gradient parameter computed as a convex combination of derivative-free LS and FR conjugate gradient parameter. The hybridization parameter of the convex combination in their work was obtained to satisfy the famous conjugacy condition. Numerical results show that the method is efficient for solving nonlinear monotone equations with convex constraints. It is noteworthy to state that, several conditions were imposed on the hybridization parameter used in ^[23] in order for the hybridization parameter to take values within the interval $(0, 1).$

Our motivation is the following: Can we extend the LS-FR method proposed by Djordjević to construct an efficient hybrid gradient-free projection algorithm where the hybridization parameter has no condition imposed on it and the hybridization parameter will always take values in the interval $[0, 1])$? In this paper, we give a positive answer to this question. The remainder of the paper is organized as follows. In Section 2, we describe the algorithm and some properties. In Section 3, we analyze the global convergence of the method. Numerical example and application are presented in Section 4 and 5 respectively.

2.   Algorithms

Consider the following unconstrained optimization problem

where $g:\mathbb{R}^n \rightarrow \mathbb{R}$ is a continuously differentiable function whose gradient at $z_k$ is denoted by $f(z_k): = \nabla(z_k)$. Given any starting point $z_0 \in \mathbb{R}^n$, the algorithm in ^[22] is to generate a sequence of approximation $\{z_k\}$ to the minimum $z^*$ of $g$, in which

where $t_k > 0$ is the steplength which is computed by a certain line search and $j_k$ is the search direction defined by

with $\beta_k$ defined by

where $\theta_k$ is a hybridization parameter chosen to satisfy the Dai-Liao's condition, that is, {for $t > 0,$}

where ${s_{k-1} = z_{k+1} -z_k}$.

Motivated by (2.3) and (2.4), we propose a gradient free projection algorithm for solving the following nonlinear equation with convex constraints:

where $\Omega \subseteq \mathbb{R}^n$ is a nonempty closed convex set, and $\rho: \mathbb{R}^n \rightarrow \mathbb{R}^n$ is a continuous mapping. Our propose gradient-free projection iterative method first generates a trial point say $\{c_k\}$ using the relation:

the search direction $j_k$ is computed by

where $\beta_k$ is computed

and $\pi_k$ is obtained to satisfy the descent condition, that is, for $\alpha > 0$,

For $k = 0$, (2.8) obviously holds. For $k\in \mathbb{N},$ we have

To satisfy (2.8), we only need that

In this paper, we choose $\pi_k$ as

It is important to note that, $\theta_k$ has the following property:

Thus,

$\theta_k = \frac{\|y_{k-1}\|^2}{y_{k-1}^T w_{k-1}^*} \in (0, 1), \quad \forall k.$

The definition of $w_{k-1}^*$ is from the ideas of Li and Fukushima ^[24,25]. The definition of $\theta_k$ was originally proposed by Birgin and Martinez ^[26] and similar idea can be found in ^[27,28] and other optimization literature. The proposed algorithm is described immediately after recalling the definition of the projection operator.

Definition 2.1. Let $\Omega \subseteq \mathbb{R}^n$ be a nonempty closed convex set. Then for any $x \in \mathbb{R}^n$, its projection onto $\Omega$, denoted by $P_\Omega[x]$, is defined by

The projection operator $P_\Omega$ has a well-known property, that is, for any $x, y \in \mathbb{R}^n$ the following nonexpansive property hold

3.   Convergence analysis

In what follows, we assume that $\rho$ satisfies the following assumptions.

Assumption 1. The solution set $\Omega^*$ is nonempty.

Assumption 2. The mapping $\rho$ is Lipschitz continuous on $\mathbb{R}^n$. That is,

Assumption 3. For any $y \in \Omega^*$ and $x \in \mathbb{R}^n,$ it holds that

Lemma 3.1. Suppose that Assumption 1 holds. Then there exists a step-size $t_k$ satisfying the line search (2.13) for $k\geq 0.$

Proof. Assume there exist $k_0 \geq 0$ such that (2.13) fails to hold for any $i\geq 0$, that is

Applying the continuity property of $\rho$ and letting $i \rightarrow \infty$ yields

which negates (2.8). Hence proved.

Lemma 3.2. Suppose Assumption 1-3 is satisfied and the sequences $\{z_k, c_k, t_k, j_k\}$ are generated by Algorithm 1. Then

Proof. Note that from (2.13), if $t_k \neq 1,$ then ${\bar t_k} = \tau^{-1}t_k$ does not satisfy $(2.13),$ that is,

Combining the above inequality with the descent condition (2.8), we have

Since $\rho$ satisfies Assumption 2 then, (3.3) holds. Thus, from $(3.3),$

This proves Lemma $3.2.$

Lemma 3.3. Suppose that Assumptions 1-3 hold and let $\{z_k\}$ and $\{c_k\}$ be the sequences generated by Algorithm 1. Then, $\rho(c_k)$ is an ascent direction of the function $\|z-z^*\|^2$ at the point $z_k,$ where $z_* \in \Omega^*$.

Proof. At $z_k$, the function $\frac{1}{2}\|x - z_*\|^2$ has a gradient of $z_k -z_*.$ By the weakly monotonicity property (3.1), it can be seen that

The inequality above, i.e., (3.5) points out that $-\rho(c_k)$ is a descent direction of the function $\|z -z_*\|$ at the iteration point $z_k.$

Lemma 3.4. Let Assumption 1-3 hold and the sequence $\{z_k\}$ be generated by Algorithm 1. Suppose that $z_*$ is a solution of problem (2.5) with $\rho(z_*) = 0.$ Then there exists a positive $\delta > 0$ such that

Proof. Remember, by using the well-known property of $P_\Omega,$ we can deduce that for any $z_* \in \Omega^*,$

From inequality (3.8) we see that $\lbrace \|z_k -z_*\| \rbrace$ is a decreasing sequence and hence $\{z_k\}$ is bounded. That is,

Furthermore, we obtain

Using the Lipchitz continuity of $\rho$, we have

Setting $\delta = L\|z_0 -z_*\|$ proves Lemma 3.4.

Lemma 3.5. Suppose Assumption 1-3 hold and the sequence $\{z_k\}$ and $\{c_k\}$ are generated by Algorithm 1. Then,

(a) $\{c_k\}$ is bounded

(b) $\lim_{k \rightarrow \infty} \|z_k -c_k \| = 0$

(c) $\lim_{k \rightarrow \infty} \|z_k -z_{k+1} \| = 0$.

Proof. (a) From $(3.10)$, we know that the sequence $\lbrace z_k \rbrace$ is bounded. So by $(3.5)$, we have

By $(3.1)$ and $(3.6)$ we have

Combined with $(3.12)$, it is easy to deduce that

Then, we obtain,

Thus $\lbrace c_k\rbrace$ is bounded due to $\lbrace z_k \rbrace$ boundedness.

(b) From inequality $(3.7)$, we get

which means

Since the mapping $\rho$ is continuous, and the $\{c_k\}$ is bounded, we know that $\{ \|\rho(c_k)\|\}$ is bounded. Therefore a positive $\delta_1 > 0$ exists, such that $\|\rho(c_k)\| \leq \delta_1$ and moreover

Hence,

Using the property of the projection operator, i.e., (2.12), we have

The global convergence result for Algorithm 1 is established via the following theorem.

Theorem 3.6. Suppose Assumption 1-3 is satisfied and the sequences $\{z_k\}$ are generated by the Algorithm 1. Then we

Proof. Suppose $(3.14)$ does not hold, meaning there exist a constant $\varepsilon_0 > 0$ such that

By $(2.8)$, we know

which implies

By $(2.3)$, we have

for all $k \in \mathbb{N}.$ Since $(3.13)$ holds, it follows that for every $\varepsilon_1 > 0$ there exist $k_0$ such that $t_{k-1} \| j_{k-1} \| < \varepsilon_1$ for every $k > k_0.$ Choosing $\varepsilon_1 = \varepsilon_0$ and $\ell_0 = \max \{ \|j_0\|, \|j_1\|, \cdots, \|j_{k_0}\|, \ell_{0_1} \}$ where $\ell_{0_1} = \delta(c + 2 +2\delta/\varepsilon_0)$, it holds that

for every $k\in \mathbb{N}.$ Integrating with $(3.4), (3.15), (3.16)$ and $(3.17)$, we know that for any $k$ sufficiently large

The last inequality yields a contradiction with $(b)$ in Lemma 3.5. Consequently, (3.14) holds. The proof is completed.

4.   Numerical section

The Dolan and Moré performance profile ^[29] is used in this section to evaluate the efficiency of the proposed algorithm on a set of test problems with varying dimensions and initial points. Comparison is made with algorithm of the same class proposed in ^[30]. All codes were written in MATLAB environment and compiled on a HP laptop (CPU Corei3-2.5 GHz, RAM 8 GB) with Windows 10 operating system.

● Algo.1: The new method (Algorithm 1).

● Algo.2: MFRM method proposed in ^[30].

The parameters for Algo.1 are chosen as: $\tau = 0.9, \; \kappa = 10^{-4}, \; \eta = 1.2$. While parameters for Algo.2 are set as reported in ^[30]. All iterative procedure are terminated whenever $\|\rho(z_k)\| < 10^{-6}.$ The experiment is carried out on nine different problems with dimensions ranging from $n = 1000, 5000, 10, 000, 50, 000, 100, 000$ using seven different initial points: $z_1 = (0.1, \cdots, 0.1)^T, \; z_2 = (0.2, \cdots, 0.2)^T, \; z_3 = (0.5, \cdots, 0.5)^T, \; z_4 = (1.2, \cdots, 1.2)^T, z_5 = (1.5, \cdots, 1.5)^T, \; z_6 = (2, \cdots, 2)^T$ and $z_7 = rand(n, 1).$ The test problems considered are listed the below where the mapping $\rho(z) = (\rho_1(z), \rho_2(z), \cdots, \rho_n(z))^T$

Problem 1 ^[31] Exponential Function.

Problem 2 ^[31] Modified Logarithmic Function.

Problem 3 ^[32]

Problem 4 ^[31] Strictly Convex Function I.

Problem 5 ^[31] Strictly Convex Function II.

Problem 6 ^[33] Tridiagonal Exponential Function.

Problem 7 ^[34] Nonsmooth Function.

Problem 8 ^[31] The Trig exp function

Problem 9 ^[35]

Figures 1-3 presents the results of the comparisons of the mentioned methods. Figure 1 shows the graph of the two methods where the performance measure is the total number of iterations. In the figure, we see that the Algo.1 obtain the most wins with the probability around 78 % and the Algo.2 method is in the second place. Figure 2 shows the performance of the considered methods relative to the total number of function evaluation. Graph of this measure shows that Algo.1 has better performance in comparison with Algo.2. In Figure 3 the performance measure is the CPU running time. The CPU running time figure also indicates that Algo.1 outperforms Algo.2. From the presented figures, it is clear that Algo.1 is the most efficient in solving the considered test problems. A detailed result of the numerical experiment for the test problems is reported in Table 2-10 in the appendix section.

5.   Application to image restoration

The restoration of images is a process in which a distorted or damaged image is restored to its original form. Having an algorithm that can perform such function with high restoration efficiency is of importance. We consider the signal-to-noise ratio (SNR), peak signal-to-noise ratio (PSNR) and the structural similarity index (SSIM) as a metric for measuring the restoration efficiency. SNR, PSNR and SSIM's larger values reflect better quality of the restored images and indicate that the restored images are closer to the original. Consider the following disturbed or incomplete observation

where $z\in \mathbb{R}^n, \; b\in \mathbb{R}^k$ is the observation data, $\rho \in \mathbb{R}^{k\times n}(k < < n)$ is a linear operator and $\omega \in \mathbb{R}^k$ is an error term. Our goal in this section is to recover the unknown vector $z.$ A well-known approach for obtaining $z$ is by solving the following $\ell_1$-regularization problem

where the regularization term $\sigma$ is positive, $\|\cdot \|_1,$ and $\|\cdot\|_2$ are the $\ell_1$-norm and $\ell_2$-norm respectively. See (Refs. ^{[36,37,38,39,40]}) for various algorithms for solving (5.2). For a comprehensive procedure on how to use our proposed algorithm to solve (5.2), see ^[41,42].

To assess the efficiency of Algo.1 in restoring the images degraded using a Gaussian blur kernel of standard deviation 0.1, we compare its performance with the modified Fletcher-Reeves conjugate Gradient method proposed in ^[30]. The algorithm is referred to as Algo.2. Four test images with different sizes are considered in this experiment. The images are labelled as A, B, C and D. The algorithms are implemented based on the following

● All codes were written and implemented in Matlab environment.

● Same starting point and stopping condition (with $Tol = 10^{-5}$) for all the algorithms.

● Parameters for Algo.1, are chosen as $\eta = 1, \; \tau = 0.55, \; \kappa = 10^{-4}.$ Parameters for Algo.2 are chosen as reported in the application section of ^[30].

● The linear operator $\rho$ in the experiment is choosen as the Gaussian matrix generated by the command $rand(k, n)$ in MATLAB.

● The signal-to-noise ratio (SNR) is defined as

where $\tilde z$ is recovered vector. The definition of the peak-to-signal and the structural similarity index (SSIM) ratio (PSNR) can be found in ^[43] and ^[44], respectively.

Figure 4 has four columns labelled ORI, BNI, RA1 and RA2. Images on the column labelled ORI are the original images, images on the column labelled BNI are the blurred and noisy images. RA1 are the images restored by Algo.1 and RA2 are images restored by Algo.2. Table 1 provides the SNR, PSNR and SSIM values for Algo.1 and Algo.2. It can be seen that Algo.1 has the highest SNR, PSNR and SSIM in all the images used for the experiment. This indicates that Algo.1 is more effective than Algo.2 in restoring blurred and noisy images.

Acknowledgments

"The authors acknowledge the support provided by the Theoretical and Computational Science (TaCS) Center under Computational and Applied Science for Smart research Innovation Cluster (CLASSIC), Faculty of Science, KMUTT. The first author was supported by the Petchra Pra Jom Klao Doctoral Scholarship, Academic for Ph.D. Program at KMUTT (Grant No.16/2561)."

Conflict of interest

The authors declare that they have no conflict of interest.

Appendix