# A Descent Conjugate Gradient Algorithm for Optimization Problems and Its Applications in Image Restoration and Compression Sensing

Authors: Junyue Cao; Jinzhao Wu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

It is well known that the nonlinear conjugate gradient algorithm is one of the effective algorithms for optimization problems since it has low storage and simple structure properties. This motivates us to make a further study to design a modified conjugate gradient formula for the optimization model, and this proposed conjugate gradient algorithm possesses several properties: (1) the search direction possesses not only the gradient value but also the function value; (2) the presented direction has both the sufficient descent property and the trust region feature; (3) the proposed algorithm has the global convergence for nonconvex functions; (4) the experiment is done for the image restoration problems and compression sensing to prove the performance of the new algorithm.


---

## Body

## 1. Introduction

Consider the following model defined by(1)minfxx∈Rn,where f:Rn⟶R is a continuous function. The above problem (1) has many practical applied fields, such as economics, biology, and engineering. It is well known that the nonlinear conjugate gradient (CG) method is one of the most effective methods for (1). The CG algorithm has the following iterative formula with(2)xk+1=xk+αkdk,k=0,1,2,…,where αk denotes the steplength, xk is the kth iterative point, and dk is the search direction designed by(3)dk=−gk+βkdk−1,ifk≥1,−gk,ifk=0,where gk=∇fxk is the gradient and βk is a scalar which determines different CG algorithms ([1–7], etc.), where the Polak–Ribière–Polak (PRP) formula [6, 7] is one of the well-known nonlinear CG formulas with(4)βkPRP=gkTgk−gk−1gk−12,where gk−1=∇fxk−1 and ⋅ is the Euclidean norm. The PRP method has been studied by many scholars, and many results are obtained (see [7–12], etc.) since the PRP algorithm has the superior numerical performance but for convergence. At present, under the weak Wolfe–Powell (WWP) inexact line search and for nonconvex functions, the global convergence of the PRP algorithm is still open, and it is one of the well-known open problems in optimization fields. Based on the PRP formula, many modified nonlinear CG formulas are done ([13–16], etc.) because many scholars want to use the perfect numerical attitude of it. Recently, Yuan et al. [17] open up a new way by modifying the WWP line search technique and partly proved the global convergence of the PRP algorithm. Further results are obtained (see [18–20], etc.) by this technique. It has been proved that the nonlinear CG algorithms can be used in nonlinear equations, nonsmooth optimization, and image restoration problems (see [21–24], etc.).We all know that the sufficient descent property designed by(5)dkTgk≤−cgk2,c>0,plays an important role for convergence analysis in CG methods (see [13, 14, 24], etc.), where c>0 is a constant. There is another crucial condition about the scalar βk≥0 that has been pointed out by Powell [10] and further emphasized in the global convergence [11, 12]. Thus, under the assumption of the sufficient descent condition and the WWP technique, a modified PRP formula βkPRP+=max0,βkPRP is presented by Gilbert and Nocedal [13], and its global convergence for nonconvex functions is established. All of these observations tell us that both property (5) and βk≥0 are very important in the CG algorithms. To get one of the conditions or both of them, many scholars made a further study and got many interesting results. Yu [25] presented a modified PRP nonlinear CG formula designed by(6)βkmPRP1=gk+1Tykgk−μykgk+1Tdkgk4,where μ>1/4 is a positive constant and yk=gk+1−gk, which has property (5) with c=1−1/4μ. Yuan [12] proposed a further formula defined by(7)βkmPRP2=gk+1Tykgk2−mingk+1Tykgk2,μykgk+1Tdkgk4,which possesses not only property (5) with c=1−1/4μ but also the scalar βkmPRP2≥0. To get a greater drop, a three-term FR CG formula is given by Zhang et al. [26]:(8)dk+1=−θkgk+1+βkFRdk,θk=dkTykgk2,where it has (5) with c=−1. Dai and Tian [27] gave another CG direction designed by(9)dk+1=−1+βkdkTgk+1gk2gk+1+βkdk,ifk≥0,−gk,ifk=0,which also possesses (5) with c=−1. The global convergence of the above CG method is proved by Dai and Tian [27] for βk=βkYuanPRP and βk=βkYuPRP. For nonconvex functions and the effective Armijo line search, they did not analyze them. One the main reasons lies in the trust region feature. To overcome it, we [28] proposed a CG formula designed by(10)dk+1=−gk,ifk=0,−gk+1+βkdk−βkdkTgk+1/gk+12gk+1γk,ifk≥0,where γk=βkdk/gk+1, which possesses not only (5) with c=−1 but also the trust region property. It has been proved that the CG formula will have better numerical performance if it possesses not only the gradient value information but also the function value [29]. This motivates us to present a CG formula based on (10) designed by(11)dk+1=−gk,ifk=0,−gk+1+βk∗dk−βk∗dkTgk+1/gk+12gk+1γk∗,ifk≥0,where γk∗=βk∗dk/gk+1, βk∗=βkmPRP2∗=gk+1Tyk∗/gk2−mingk+1Tyk∗/gk2,μyk∗gk+1Tdk/gk4, and yk∗=yk+ρk with ρk=maxϱk,0/sk2 and ϱk=2fxk−fxk+1+gk+1+gkTsk and sk=xk+1−xk. The new vector yk∗ [30] has been proved that it has some good properties in theory and experiment. Yuan et al. [29] use it in the CG formula and get some good results. These achievements inspire us to propose the new CG direction (11), and this paper possesses the following features:(i)
The sufficient property and the trust region feature are obtained(ii)
The new direction possesses not only the gradient value but also the function value(iii)
The given algorithm has the global convergence under the Armijo line search for nonconvex functions(iv)
The experiments for image restoration problems and compression sensing are done to test the performance of the new algorithmThe next section states the given algorithm. The convergence analysis is given in Section3, and experiments are done in Section 4, respectively. The last section proposes one conclusion.


## 2. Algorithm

Based on the discussions of the above section, the CG algorithm is listed in Algorithm1.Algorithm 1: Three-term conjugate gradient algorithm.
Initial step: given any initial pointx0∈Rn and positive constants ϵ∈0,1,σ0>0,δ∈0,1,μ>1/4, σ∈0,1, set d0=−g0=−∇fx0 and k:=0.Step 1: stop ifgk≤ε is true.Step 2: findαk=σ0σik such thatfxk+αkdk≤fxk+δαkgxkTdk,whereik=min0,1,2,… satisfying the equation in Step 2 mentioned in Algorithm 1.Step 3: setxk+1=xk+αkdk.Step 4: stop ifgk+1≤ε holds.Step 5: computedk by (11).Step 6: setk:=k+1 and go to Step 2.Theorem 1.
The directiondk is defined by (11); then there exists a positive β>0 satisfying(12)dkTgk=−gk2,∀k≥0,(13)dk≤βgk,∀k≥0.Proof.
By (11), we directly get (12) and (13) for k=0 with β=−1. If k≥0, using (11) again, we have(14)gk+1Tdk+1=gk+1T−gk+1+βkmPRP2∗dk−βkmPRP2∗dkTgk+1/gk+12gkγk∗=−gk+1Tgk+1+βkmPRP2∗gk+1Tdk−βkmPRP2∗dkTgk+1/gk+12gk+1Tgk+1βkmPRP2∗dk/gk+1=−gk+12+βkmPRP2∗gk+1Tdk−βkmPRP2∗dkTgk+1βkmPRP2∗dk/gk+1=−gk+12,then (12) is true. By (11) again, we can get(15)dk+1=−gk+1+βkmPRP2∗dk−βkmPRP2∗dkTgk+1gk+12gk+1/γk∗≤gk+1+βkmPRP2∗dk+βkmPRP2∗dkgk+1/gk+12gk+1βkmPRP2∗dk/gk+1=gk+1+2βkmPRP2∗dkβkmPRP2∗dk/gk+1=3gk+1,which implies that (13) holds by choosing β∈3,+∞. We complete the proof.Remark 1.
The relation (13) is the so-called trust region feature, and the above theorem tells us that direction (11) has not only the sufficient descent property but also the trust region feature. Both these relations (12) and (13) will make the proof of the global convergence of Algorithm 1 be easy to be established.


## 3. Global Convergence

For the nonconvex functions, the global convergence of Algorithm1 is established under the following assumptions.Assumption 1.
Assume that the functionfx has at least a stationary point x∗, namely, gx∗=0 is true. Suppose that the level set defined by L0=xfx≤fx0 is bounded.Assumption 2.
The functionfx is twice continuously differentiable and bounded below, and its gx is Lipschitz continuous. We also assume that there exists a positive constant L>0 such that(16)gx−gy≤Lx−y,x,y∈Rn.
Now, we prove the global convergence of Algorithm1.Theorem 2.
Let Assumption1 be true. Then, we get(17)limk⟶∞gk=0.Proof.
Using (12) and the Step 2 of Algorithm 1, we obtain(18)fxk+αkdk≤fxk+δαkgxkTdk<fk−δαkgxk2,which means that the sequence fxk is descent and the following relation(19)δαkgxk2≤fxk−fxk+αkdk,is true. For k=0 to ∞, by summing the above inequalities and Assumption 1, we deduce that(20)∑k=0∞δαkgxk2≤fx0−f∞<+∞,holds. Thus, we have(21)limk⟶∞αkgxk2=0.
This implies that(22)limk⟶∞gxk=0,or(23)limk⟶∞αk=0.
Suppose that (22) holds, the proof of this theorem is complete. Assuming that (23) is true, we aim to get (17). Let the stepsize αk satisfy the equation in Step 2 in Algorithm 1; for αk∗=αk/σ, we have(24)fxk+αk∗dk>fxk+δαk∗dkTgxk.
By (12) and (13) and the well-known mean value theorem, we obtain(25)fxk+αk∗dk−fxk=αk∗dkTgxk+αk∗2Odk2=−αkσgxk2+αk2σ2Odk2>δαk∗dkTgxk=−δαkσgxk2,which implies that(26)αk>σ1−δgxk2Odk2≥σ1−δO1β,is true. This is a contradiction to (23). Then, only relation (22) holds. We complete the proof.Remark 2.
We can see that the proof process of the global convergence is very simple since the defined direction (11) has not only the good sufficient descent property (12) but also the perfect trust region feature (13).


## 4. Numerical Results

The numerical experiments for image restoration problems and compression sensing will be done by Algorithm1 and the normal PRP algorithm, respectively. All codes are run on a PC with an Intel (R) Core (TM) i7-7700T CPU @ 2.9 GHz, 16.00 GB of RAM, and the Windows 10 operating system and written by MATLAB r2014a. The parameters are chosen as σ=0.5, σ0=0.1, δ=0.9, and μ=300.


### 4.1. Image Restoration Problems

Settingx be the true image which has M×N pixels and i,j∈A=1,2,…,M×1,2,3,…,N. At a pixel location i,j, xi,j denotes the gray level of x. Then, defining a set N by(27)N≔i,j∈Aζ¯i,j≠ζi,j,ζi,j=sminorsmax,which is the index set of the noise candidates. Suppose that ζ is the observed noisy image of x corrupted by salt-and-pepper noise, we let ϕi,j=i,j−1,i,j+1,i−1,j,i+1,j be the neighborhood of i,j. By applying an adaptive median filter to the noisy image y, ζ¯ is defined by the image obtained. smax denotes the maximum of a noisy pixel, and smin denotes the minimum of a noisy pixel. The following conclusions can be obtained: (i) if i,j∈N, then ζi,j must be restored. A pixel i,j is identified as uncorrupted, and its original value is kept, which means that wi,j∗=ζi,j with the element wi,j∗ of the denoised image w by the two-phase method. (ii) If i,j∉N holds, wm,n∗=ζm,n is stetted and ζm,n is restored if m,n∈ϕi,j∩N. Chan et al. [31] presented the new function fα and minimized it for the restored images without a nonsmooth term, which has the following form:(28)fαw=∑i,j∈N∑m,n∈ϕi,j∖Nψαwi,j−ζm,n+12∑m,n∈ϕi,j∩Nψαwi,j−ζm,n,where α is a constant and ψα is an even edge-preserving potential function. The numerical performance of fα is noteworthy [32, 33].We choose Barbara512×512, man 256×256, Baboon 512×512, and Lena 256×256 as the tested images. The well-known PRP CG algorithm (PRP algorithm) is also done to compare with Algorithm 1. The detailed performances are listed in Figures 1 and 2.Figure 1
Restoration of Barbara, man, Baboon, and Lena by Algorithm1 and the PRP algorithm. From left to right: a noisy image with 20% salt-and-pepper noise and restorations obtained by minimizing z with the PRP algorithm and Algorithm 1.Figure 2
Restoration of Barbara, man, Baboon, and Lena by Algorithm1 and the PRP algorithm. From left to right: a noisy image with 40% salt-and-pepper noise and restorations obtained by minimizing z with the PRP algorithm and Algorithm 1.Figures1 and 2 tell us that these two algorithms (Algorithm 1 and the PRP algorithm) are successful to solve these image restoration problems, and the results are good. To directly compare their performances, the restoration performance is assessed by applying the peak signal-to-noise ratio (PSNR) defined in [34–36] which is computed and listed in Table 1. From the value of Table 1, we can see that Algorithm 1 is competitive to the PRP algorithm since its PSNR value is less than that of the PRP algorithm.Table 1
PSNR Algorithm1 and the PRP algorithm.
20% noiseBarbaraManBaboonLenaAverageAlgorithm131.11538.035529.439341.067434.9143PRP algorithm31.111837.958329.453441.35634.96940% noiseBarbaraManBaboonLenaAverageAlgorithm127.541534.006325.894736.649631.0230PRP algorithm27.615334.537525.857136.70131.1777


### 4.2. Compressive Sensing

In this section, the following compressive sensing images are tested: Phantom256×256, Fruits 256×256, and Boat 256×256. These three images are treated as 256 vectors, and the size of the observation matrix is 100×256. The so-called Fourier transform technology is used, and the measurement is the Fourier domain.Figures3–5 turn out that these two algorithms work well for these figures, and they can successfully solve them.Figure 3
Phantom: (a) the general images, (b) the recovered images by Algorithm1, and (c) the recovered images by the PRP algorithm.
(a)(b)(c)Figure 4
Fruits: (a) the general images, (b) the recovered images by Algorithm1, and the (c) recovered images by the PRP algorithm.
(a)(b)(c)Figure 5
Boat: (a) the general images, (b) the recovered images by Algorithm1, and (c) the recovered images by the PRP algorithm.
(a)(b)(c)


## 4.1. Image Restoration Problems

Settingx be the true image which has M×N pixels and i,j∈A=1,2,…,M×1,2,3,…,N. At a pixel location i,j, xi,j denotes the gray level of x. Then, defining a set N by(27)N≔i,j∈Aζ¯i,j≠ζi,j,ζi,j=sminorsmax,which is the index set of the noise candidates. Suppose that ζ is the observed noisy image of x corrupted by salt-and-pepper noise, we let ϕi,j=i,j−1,i,j+1,i−1,j,i+1,j be the neighborhood of i,j. By applying an adaptive median filter to the noisy image y, ζ¯ is defined by the image obtained. smax denotes the maximum of a noisy pixel, and smin denotes the minimum of a noisy pixel. The following conclusions can be obtained: (i) if i,j∈N, then ζi,j must be restored. A pixel i,j is identified as uncorrupted, and its original value is kept, which means that wi,j∗=ζi,j with the element wi,j∗ of the denoised image w by the two-phase method. (ii) If i,j∉N holds, wm,n∗=ζm,n is stetted and ζm,n is restored if m,n∈ϕi,j∩N. Chan et al. [31] presented the new function fα and minimized it for the restored images without a nonsmooth term, which has the following form:(28)fαw=∑i,j∈N∑m,n∈ϕi,j∖Nψαwi,j−ζm,n+12∑m,n∈ϕi,j∩Nψαwi,j−ζm,n,where α is a constant and ψα is an even edge-preserving potential function. The numerical performance of fα is noteworthy [32, 33].We choose Barbara512×512, man 256×256, Baboon 512×512, and Lena 256×256 as the tested images. The well-known PRP CG algorithm (PRP algorithm) is also done to compare with Algorithm 1. The detailed performances are listed in Figures 1 and 2.Figure 1
Restoration of Barbara, man, Baboon, and Lena by Algorithm1 and the PRP algorithm. From left to right: a noisy image with 20% salt-and-pepper noise and restorations obtained by minimizing z with the PRP algorithm and Algorithm 1.Figure 2
Restoration of Barbara, man, Baboon, and Lena by Algorithm1 and the PRP algorithm. From left to right: a noisy image with 40% salt-and-pepper noise and restorations obtained by minimizing z with the PRP algorithm and Algorithm 1.Figures1 and 2 tell us that these two algorithms (Algorithm 1 and the PRP algorithm) are successful to solve these image restoration problems, and the results are good. To directly compare their performances, the restoration performance is assessed by applying the peak signal-to-noise ratio (PSNR) defined in [34–36] which is computed and listed in Table 1. From the value of Table 1, we can see that Algorithm 1 is competitive to the PRP algorithm since its PSNR value is less than that of the PRP algorithm.Table 1
PSNR Algorithm1 and the PRP algorithm.
20% noiseBarbaraManBaboonLenaAverageAlgorithm131.11538.035529.439341.067434.9143PRP algorithm31.111837.958329.453441.35634.96940% noiseBarbaraManBaboonLenaAverageAlgorithm127.541534.006325.894736.649631.0230PRP algorithm27.615334.537525.857136.70131.1777


## 4.2. Compressive Sensing

In this section, the following compressive sensing images are tested: Phantom256×256, Fruits 256×256, and Boat 256×256. These three images are treated as 256 vectors, and the size of the observation matrix is 100×256. The so-called Fourier transform technology is used, and the measurement is the Fourier domain.Figures3–5 turn out that these two algorithms work well for these figures, and they can successfully solve them.Figure 3
Phantom: (a) the general images, (b) the recovered images by Algorithm1, and (c) the recovered images by the PRP algorithm.
(a)(b)(c)Figure 4
Fruits: (a) the general images, (b) the recovered images by Algorithm1, and the (c) recovered images by the PRP algorithm.
(a)(b)(c)Figure 5
Boat: (a) the general images, (b) the recovered images by Algorithm1, and (c) the recovered images by the PRP algorithm.
(a)(b)(c)


## 5. Conclusion

This paper, by designing a CG algorithm, studies the unconstrained optimization problems. The given method possesses not only the sufficient descent property but also the trust region feature. The global convergence is proved by a simple way. The image restoration problems and compressive sensing problems are tested to show that the proposed algorithm is better than the normal algorithm. In the future, we will focus on the following aspects to be paid attention: (i) we believe there are many perfect CG algorithms which can be successfully used for image restoration problems and compressive sensing; (ii) more experiments will be done to test the performance of the new algorithm.



---

*Source: 6157294-2020-09-29.xml*