A NEW RELAXED CQ ALGORITHM FOR SOLVING SPLIT FEASIBILITY PROBLEMS IN HILBERT SPACES AND ITS APPLICATIONS

. Inspired by the works of L´opez et al. [21] and the recent paper of Dang et al. [15], we devise a new inertial relaxation of the CQ algorithm for solving Split Feasibility Problems (SFP) in real Hilbert spaces. Under mild and standard conditions we establish weak convergence of the proposed algorithm. We also propose a Mann-type variant which converges strongly. The performances and comparisons with some existing methods are presented through numerical examples in Compressed Sensing and Sparse Binary Tomography by solving the LASSO problem.

1. Introduction. In this paper we are concerned with the Split Feasibility Problem (SFP) in real Hilbert spaces. Let C and Q be nonempty, closed and convex subsets of the real Hilbert spaces H 1 and H 2 , respectively. In addition, there is a bounded and linear operator A : H 1 → H 2 . The SFP is formulated as finding a point x * such that x * ∈ C and Ax * ∈ Q.
(1) We denote by Ω the solution set of (1), i.e., Ω := {x * ∈ C : Ax * ∈ Q}. ( The SFP was introduced by Censor and Elfving in [10], and since then it has been studied extensively and generalized in various ways for finite and infinite dimensional spaces. Censor et. al. [11,12] used the SFP as a modelling tool in the field of Intensity-Modulated Radiation Therapy (IMRT) thereanent planning. Due to its nature, the SFP has been applied successfully in many real-world problems, such as for signal processing, image reconstruction and many more, for theory and application, the reader is referred to [2,7,8,10,9,14,21,33,35,38] and the references therein.
Censor and Elfving [10] original algorithm for solving the SFP involves the computation of the inverse of A in (1) per each iteration, a fact that makes the algorithm non applicable in practice. Byrne in [7,8] proposed the CQ Algorithm that uses the adjoint (transpose) of A, that is A * (A T ) instead and hence make the algorithm very attractive. The iterative step of the CQ Algorithm is as follows.
where γ ∈ 0, 2 and P C and P Q are the orthogonal projections onto C and Q, respectively. In case when the sets C and Q are not "simple" to project onto, there is the need to solve a quadratic minimization problem twice per each iteration in order to evaluate P C and P Q , a fact that might effect the applicability of the method. In order to overcome this drawback, Yang [38] considered split feasibility problems in which the involved sets C and Q are given as sub-level sets of convex functions, i.e., where c : H 1 → R and q : H 2 → R are convex functions. Yang [38] introduced a relaxed CQ algorithm in a finite-dimensional Hilbert space as follows: where C k and Q k are given as and Q k = {y ∈ H 2 : q(Ax k ) ≤ ζ k , Ax k − y }, where ξ k ∈ ∂c(x k ) and ζ k ∈ ∂q(Ax k ) are subgradients of c and q at x k and Ax k , respectively; this will be explained further in the sequel. For more details on the SFP and in particular the CQ algorithm, the reader might referred to [2,7,8,9,14,16,21,33,34,35,38] and the many references therein.
Observe that the sets C k and Q k are either halfspaces or the all spaces, and then the projection is very simple and has a close formula, a fact which makes this projection method (5) very attractive in practice. A drawback of the above mention algorithms, both (3) and (5), is the need to calculate the stepsize γ which depends on the operator (matrix) norm A . This means that the implementation of algorithms (3) and (5), one needs to know (or at least approximate) a priori the norm of A, which is not an "easy" task in general and hence might affect the convergence of these algorithms. In order to overcome this difficulty, López et al. [21] introduced a new relaxation algorithm whose iterative step is as follows: where the stepsize τ k is calculated as where with 0 < ρ k < 4 and inf ρ k (4 − ρ k ) > 0. López et al. [21] show that any sequence {x k } ∞ k=0 generated by (8)-(10) converges weakly to a solution x * of the SFP (1).
Very recently, Dang et al. [15] applied the inertial technique of Alvarez and Attouch [1] to Yang's algorithm ((5)) and propose an inertial relaxed CQ algorithm for solving the SFP. Dang et al. iterative step is based on the following: choose arbitrary starting points x 0 , x 1 ∈ H, Although it is clear that this algorithm has a better theoretical properties, such as, but not only, weak convergence to a solution of the SFP (1) in real Hilbert spaces, there is still the need to calculate the stepsize γ which depends on the operator norm.
Hence, it is only natural to ask the following question.
Question 1. Can algorithm (11) with the inertial technique be combined with the self-adaptive rule used in (9)-(10)?
Motivated and inspired by the works of Alvarez and Attouch [1], Dang et al. [15] and Lopéz et al. [21], as well as of Shehu and Iyiola [30,29], we wish to provide an affirmative answer to this question. Our contribution is a new iterative method which is called self-adaptive relaxed CQ algorithm for solving the split feasibility problem 1 in real Hilbert spaces.
The outline of the paper is as follows. In Section 2, we collect definitions and results which are needed for our further analysis. In Section 3, our new self-adaptive CQ algorithm is introduced and analyzed. In Section 4 we propose a Mann-type variant of the self-adaptive CQ algorithm and establish its strong convergence. Finally, in Section 5 we illustrate the performances of the proposed algorithm and report some preliminary numerical results for solving the LASSO problem.
2. Preliminaries. Throughout this paper we denote by H 1 and H 2 real Hilbert spaces with the inner product ., . , the induced norm . and the identity operator I. Given a bounded and linear operator A : A useful and simple norm equality is the following for all x, y, z ∈ H and for all α, β, γ ∈ [0, 1] with α + β + γ = 1. Let C be a closed convex subset of H. For every element x ∈ H, there exists a unique nearest point in C, denoted by P C x such that

AVIV GIBALI, DANG THI MAI AND NGUYEN THE VINH
The operator P C is called the metric projection of H onto C and some of its properties are summarized in the next lemma, see e.g., [18].
Lemma 2.1. Let C ⊆ H be a closed convex set, P C fulfils the following: Definition 2.3. Let h : H → R be a convex function. The subdifferential set (or just subdifferential) of h is defined as Each element of ∂h(x) is called subgradient. In case that the function h is continuously differentiable then ∂h(x) = {∇h(x)}, this is the gradient of h.
Next we wish to present some properties which link directly to the SFP (1), for this we assume that the sets C ⊆ H 1 and Q ⊆ H 2 are nonempty, closed convex sets and A : H 1 → H 2 is a bounded and linear operator.
(1) The point x * solves the SFP (1); (2) The point x * solves the fixed point equation (3) The point x * solves the variational inequality problem with respect to the gradient of f , that is: find a point x ∈ C such that The next following lemmas are essential and be used in our sequel analysis of the proposed algorithms.
The next well-known result is due to Opial [28].
Then the following results hold: n=0 be a sequence of real numbers that does not decrease at infinity, in the sense that there exists a subsequence Then {τ (n)} n≥n0 is a nondecreasing sequence verifying lim n→∞ τ (n) = ∞ and, for all n ≥ n 0 , 3. The self-adaptive relaxed CQ algorithm.
3.1. The algorithm. In this section, we introduce a new self-adaptive relaxation CQ algorithm for solving SFP (1) in which we assume that the sets C ⊆ H 1 and Q ⊆ H 2 are given by where c : H 1 → R and q : H 2 → R are convex and lower semicontinuous functions. We assume that ∂c and ∂q are bounded operators (i.e., bounded on bounded sets). This assumption is a quite standard assumption and is hold for example when c and q are bounded on bounded sets, see [5, Proposition 1.1.11] and also [6,Proposition 7.8]. Set where ξ k ∈ ∂c(w k ), and where ζ k ∈ ∂q(Aw k ). Obviously, C k and Q k are half-spaces and moreover C ⊆ C k and Q ⊆ Q k for every k ≥ 0. Now we define the functions and can easily see that Algorithm 1. (Self-adaptive relaxed CQ algorithm) Select arbitrary starting points x 0 , x 1 ∈ C, θ ∈ [0, 1) and set k = 0. Iterative Step: Given the iterates x k−1 and Compute and calculate the next iterate as Stop Criterion: If x k+1 = w k then stop. Otherwise, set k := k + 1 and return to Iterative Step.
Remark 1. Evidently, we have from (30) that 3.2. Convergence analysis. We first wish to validate the stoping criterion of Algorithm 1.
Proof. Since x k+1 = w k , we get by (33) we have Since C ⊆ C k , we get w k ∈ C k . Combined with the fixed point relation Lemma 2.5(2), we also get that Aw k ∈ Q k . Following the representations of the sets C k and Q k ( (24) and (25)) we obtain that c(w k ) ≤ 0 and q(Aw k ) ≤ 0, and this implies that w k ∈ C and Aw k ∈ Q, which completes the proof. Now, we start analyzing the convergence of Algorithm 1.

Lemma 3.2.
Consider the SFP (1) with its solution set Ω. Let {x k } ∞ k=0 be any sequence generated by Algorithm 1. Then, for each z ∈ Ω, the following inequality holds: Proof. Denoting by Lemma 2.1 (2) and (33) we have So, and By (38), (39) and (40) we get Estimate by the Cauchy-Schwarz inequality the third term on the right-hand side of (41) yields: Consequently, we obtain And the proof is complete.
The main weak convergence theorem of Algorithm 1 is presented next.
Theorem 3.3. Consider the SFP (1) with its solution set Ω and assume that lim inf k→∞ ρ k (4 − ρ k ) > 0. Then any sequence {x k } ∞ k=0 generated by Algorithm 1 converges weakly to an element of Ω.
Proof. Let z ∈ Ω. From Lemma 3.2, we have On the other hand by applying (13) we get Therefore, Thus we get and Applying Lemma 2.6 in (50) with the data where ∞ k=1 δ k < ∞ (by (34)), we have that lim k→∞ x k − z exists. This leads to the boundedness of {x k } ∞ k=0 and therefore {w k } ∞ k=0 is also bounded. Passing to the limit in (51) and taking into account that {φ k } ∞ k=0 converges and δ k go to zero as k tends to ∞, we obtain The assumption lim inf Since ∇f k is Lipschitz, we have Together with (54), this implies that It remains to show that ω w (x k ) ⊂ Ω. Letx ∈ ω w (x k ) be an arbitrary element. Since {x k } ∞ k=0 is bounded, there exists a subsequence {x k l } ∞ l=0 of {x k } ∞ k=0 which converges weakly tox. On the other hand, from (31) we have It follows from (57) that w k l x ∈ C. Since P Q k l Aw k l ∈ Q k l , we have where ζ k l ∈ ∂q(Aw k l ). From the boundedness assumption of ζ k l and (56), we have From the weak lower-semicontinousness of the convex function q(x) and since w k l x, it follows from (90) that which means that Ax ∈ Q. Next, we show thatx ∈ C. First, we prove that Indeed, it is easy to see that Combining (48) and (62), we have It follows from x k+1 = P C k (w k − λ k ∇f k (w k )) that Note that Then Combining (63) and (66) we get Next, using the fact that x k l +1 ∈ C k l and by the definition of C k l , we get where ξ k l ∈ ∂c(w k l ). Following the assumption on the boundedness of ξ k l and (67), we have as l → ∞. Similarly, we obtain that c(x) ≤ 0, i.e.,x ∈ C. We arrive atx ∈ Ω. The choice ofx in ω w (x k ) was arbitrary, and so we conclude that ω w (x k ) ⊂ Ω. Hence, it follows from Lemma 2.7 that the result holds and the proof is complete. 4. Strong convergence theorem. In this section, motivated by the works of Li and Yao [20] and Suantai et al. [31] we propose a Mann-type variant of Algorithm 1 with strong convergence property. For this purpose, we define the following convex and differentiable functions and their gradients are Algorithm 2. (Self-adaptive relaxed CQ algorithm) where k = o(γ k ) means that the sequence { k } ∞ k=0 is an infinitesimal of higher order than {γ k } ∞ k=0 . Select arbitrary starting points x 0 , x 1 ∈ C, θ ∈ [0, 1) and set k = 0. Iterative Step: Given the iterates x k−1 and x k (k ≥ 1), choose α k such that 0 ≤ α k ≤ᾱ k , wherē Compute If ∇f k (w k ) = ∇g k (w k ) = 0 then stop (w k is a solution of the SFP (1)). Otherwise, compute and Let k := k + 1 and return to Iterative Step.
Remark 2. (i) In Algorithm 2, the choice of α k depends on the behaviour of the two previous iterates (x k , x k−1 ), while in [31, Algorithm 3.1], θ k ∈ [0, 1) is fixed and moreover, the next iterate in our case is not calculate as a convex combination of two terms.
(ii) It follows from (73) and (74) that For the convergence analysis of Algorithm 2, we need the following result.
Lemma 4.1. Let {x k } ∞ k=0 be any sequence generated by Algorithm 2, Then, for each z ∈ Ω, the following inequality holds: where Using similar arguments as in the proof of Lemma 3.2 we obtain Hence, for all k ∈ N, We note that Using the inequality x + y 2 ≤ x 2 + 2 y, x + y we get and the desired result is obtained.
For the algorithm's strong convergence theorem, which we present next, we recall the minimum-norm element of Ω which is a solution to the following problem. argmin{ x : x solves the SFP (1)}. (85) generated by Algorithm 2 converges strongly to the minimum-norm element of Ω.
Proof. We first show that {x k } ∞ k=0 is bounded. Let z ∈ Ω. From (77) we obtain where
Using Lemma 4.1 we have Besides, we obtain where Combining (97) and (98) we get It follows from (73) and (74) that From (96) we also have To apply Lemma 2.8, it remains to show that lim sup k→∞ x k+1 − z, −z ≤ 0. Indeed, since z = P Ω (0), by the characterization of the metric projection (Lemma 2.1 (1)), we get that By applying Lemma 2.8 to (34) with the data: we immediately deduce that the sequence {x k } ∞ k=0 converges strongly to z = P Ω (0). Furthermore, it follows again from Lemma 2.1 (1) that Hence z 2 ≤ p, z ≤ z p ∀p ∈ Ω, from which we infer that z is the minimum-norm solution of the SFP (1).

Numerical experiments.
In the next subsections we present two numerical examples which demonstrate the performances of our proposed scheme, Algorithm 1. The first example is in the field of compressed sensing and it is the recovery of a sparse and noisy signal from a limited number of sampling. In the second example we consider an image reconstruction from few tomographic projections, in particular a vascular system containing larger and smaller vessels.

Compressed sensing.
Here we consider the problem of recovering a noisy sparse signal from a limited number of sampling. The problem description as is follows. Let x 0 ∈ R n be a K-sparse signal, K << n. The sampling matrix A ∈ R m×n (m < n) is stimulated by standard Gaussian distribution and vector b = Ax 0 + e, where e is additive noise. When e = 0, it means that there is no noise to the observed data. Our task is to recover the signal x 0 from the data b. For further explanations the reader can consult [27].
For solving this problem we make use of the model of Tibshirani [32] (strongly related to the Basis Pursuit denosing problem [13]), which is known as the LASSO problem. Close and related approaches can be found in the recent results of [26,16].
where t > 0 is a given constant. So, in relation with the SFP (1) we consider C = {x : x 1 ≤ t} and Q = {b}. We also define the convex function c(x) = x 1 −t and consider the level set C k as in (24). Observe that the orthogonal projection onto C k can be calculated by the following, Following Definition 2.3, that is of the subdifferential set ∂c(x k ), we choose subgradient ξ k ∈ ∂c(x k ) as In the implementation of Algorithm 1 we choose the following set of parameters: the matrix A ∈ R m×n is generated randomly with m = 2 10 and n = 2 12 , x 0 ∈ R n with k spikes with amplitude ±1 distributed in the whole domain randomly. Moreover, ε k = 1/k 2 , θ = 0.9, ρ k ≡ 2 and α k = max{0, α k − 0.1}.
We compare the performances of Algorithm 1 with two relaxed algorithms with self adaptive step size, Lopez et al. algorithm [21] ((8)-(10)) and Zhou and Wang algorithm [39,Algorithm 3.1] in which the step size τ k of (9) is calculated as where σ k ≡ 0.5, f k and ∇f k are as in (10). In addition we choose = 10 −6 and the stopping criterion for our algorithm is x k+1 − w k ≤ , for Lopez et al.
∇f (x k ) ≤ and for Zhou and Wang algorithm x k+1 − x k / x k ≤ . Starting point for all algorithms is 0 ∈ R n and we take t = k. All computations were performed using MATLAB R2015a on an Intel Core i5-4200U 2.3GHz running 64bit Windows.
Next, in Table 1 we present our experiments with m = 2 10 , n = 2 12 for K = 20, 40 and also m = 2 12 , n = 2 13 for K = 50. Later, in Figures 1, 2 and 3 we illustrate and compared the recovered signal x 0 and the discrepancy of f k for each example and each algorithm. Iter f k (10) CPU time (sec.) Each pixel is associated with some unknown binary variable u i ∈ {0, 1}. Each entry in q ∈ R m , called tomographic measurement or single projection, corresponds to the integrated gray values of u along the single ray, see Figure 4. Hence the integral can be split into the sum of products M ij u j , where each matrix entry M ij ≥ 0 corresponds to the length of the intersection of the i-th ray with the j-th pixel. If ray i and pixel j do not intersect then M ij = 0. Stacking all equations for all the rays together leads to the linear equations in (119) and provides the following representation of the reconstruction problem Since problem (120) is NP-complete for more than two projections, alternative approaches have been suggested, for example the DC (Difference of Convex) programming and l 0 -superiorization, see Gibali and Petra [17] and l 1 -l 2 minimization, see [22,37] and [26].       So, here we apply the following l 1 approach min 1 2 M u − q 2 subject to , u ∈ [0, 1] n and u 1 ≤ t where [0, 1] n is the [0, 1] cube in R n and is a continues relaxation of the discrete set {0, 1} n . For the illustration we consider the reconstruction of the binary test image from [4] representing a vascular system containing larger and smaller vessels, see Figure 5 (left). and each block matrix M θi corresponds to a different projecting angle, see Figure 5 (right) for a parallel ray geometry that corresponds to one such angle. We use the MATLAB routine paralleltomo.m from the AIR Tools package [19] that implements such a tomographic matrix for a given vector of angles. We choose a square (N = 64) 64 × 64 size image and use the default value of p, i.e. the number of parallel rays for each angle p = round(sqrt(2)*N). So, the parameters in use are the same as in the previous subsection and m = 1365, n = 4096 and k = 916, 100 iterations and 15 projections. The results for all three algorithms are presented in Figure 6 as well as in Table 2.
K and m, n Methods = 10 −6 Iter f k (10) CPU time (sec.) Conclusions. In this research article we propose a new relaxed CQ algorithm for solving split feasibility problems in real Hilbert spaces. We combine the Dang et al. [15] inertial type algorithm (11) with the self-adaptive rule used by Lopéz et al. [21] (8)- (10). The algorithm advantages is that outer approximations of the sets C, Q are used and hence make the orthogonal projections simple for computation, moreover, there is no need to calculate or estimate the operator norm a-priori. The proposed algorithm converges weakly to a solution. For strong convergence, we also propose a Mann-type variant of our algorithm. We present two numerical examples in the Rescaled Original Lopez et al. Figure 6. Reconstructing a binary test image u ∈ R 64×64 from [4] representing a vascular system containing larger and smaller vessels. The results are for the recover u from a 15 (limited number of) tomographic projections.

Relaxed CQ Zhou and Wang
fields of Compressed Sensing and Sparse Binary Tomography which demonstrate the applicability and the potential advantages of our proposed schemes compared with some existing works.