On Two Projection Algorithms for the Multiple-Sets Split Feasibility Problem

Abstract

We present a projection algorithm which modifies the method proposed by Censor and Elfving (1994) and also introduce a self-adaptive algorithm for the multiple-sets split feasibility problem (MSFP). The global rates of convergence are firstly investigated and the sequences generated by two algorithms are proved to converge to a solution of the MSFP. The efficiency of the proposed algorithms is illustrated by some numerical tests.

1. Introduction

The multiple-sets split feasibility problem (MSFP) is to find 𝑥^∗ satisfying 𝑥^∗∈𝐶∶=^𝑡_𝑖=1𝐶_𝑖suchthat𝐴𝑥^∗∈𝑄∶=^𝑟_𝑗=1𝑄_𝑗, (1) where 𝐴 is an 𝑀×𝑁 real matrix, 𝐶_𝑖⊆ℝ^𝑁, 𝑖=1,…,𝑡, and 𝑄_𝑗⊆ℝ^𝑀, 𝑗=1,…,𝑟, are the nonempty closed convex sets. This problem was firstly proposed by Censor et al. in [1] and can be a model for many inverse problems where constraints are imposed on the solutions in the domain of a linear operator as well as in the operator's range. Many researchers studied the MSFP and introduced various algorithms to solve it (see [1–7] and the references therein). If 𝑡=𝑟=1, then this problem reduces to the feasible case of the split feasibility problem (see, e.g., [8–11]), which is to find 𝑥^∗∈𝐶 with 𝐴𝑥^∗∈𝑄.

Assume that the MSFP (1) is consistent; that is, its solution set, denoted by Γ, is nonempty. For convenience reasons, Censor et al. [1] considered the following constrained MSFP: find𝑥^∗∈Ωsuchthat𝑥^∗solvestheMSFP, (2) where Ω⊆ℋ₁ is an auxiliary simple nonempty closed convex set containing at least one solution of the MSFP. For solving the constrained MSFP, Censor et al. [1] defined a proximity function 𝑝(𝑥) to measure the distance of a point to all sets 𝑝(𝑥)=12^𝑡_𝑖=1𝛼_𝑖𝑥−𝑃_𝐶𝑖(𝑥)²+12^𝑟_𝑗=1𝛽_𝑗𝐴𝑥−𝑃_𝑄𝑗𝐴(𝑥)², (3) where 𝛼_𝑖>0 and 𝛽_𝑗>0 for all 𝑖 and 𝑗, respectively, and ∑^𝑡_𝑖=1𝛼_𝑖+∑^𝑟_𝑗=1𝛽_𝑗=1. We see that ∇𝑝(𝑥)=^𝑡_𝑖=1𝛼_𝑖𝑥−𝑃_𝐶𝑖(𝑥+_𝑟𝑗=1)𝛽_𝑗𝐴^𝑇𝐼−𝑃_𝑄𝑗𝐴𝑥. (4) Censor et al. [1] proposed a projection algorithm as follows: 𝑥_𝑛+1=𝑃_Ω𝑥_𝑛−𝑠∇𝑝𝑥_𝑛, (5) where 𝑠 is a positive number such that 0<𝑠_𝐿≤𝑠≤𝑠_𝑈<2/𝐿(𝑝) and 𝐿(𝑝) is the Lipschitz constant of ∇𝑝.

Observe that in the algorithm (5) the determination of the stepsize 𝑠 depends on the operator (matrix) norm ‖𝐴‖ (or the largest eigenvalue of 𝐴∗𝐴). This means that, in order to implement the algorithm (5), one has first to compute (or, at least, estimate) operator norm of 𝐴, which is in general not an easy work in practice. To overcome this difficulty, Zhang et al. [4] and Zhao and Yang [5, 7] proposed self-adaptive methods where the stepsize has no connection with matrix norms. Their methods actually compute the stepsize by adopting different self-adaptive strategies.

Note that the algorithms proposed by Censor et al. [1], Zhang et al. [4], and Zhao and Yang [5, 7] involve the projection to an auxiliary set Ω. In fact, the set Ω is introduced just for the convenience of the proof of the convergence and it may be difficult to determine Ω in some cases. Considering this, Zhao and Yang [6] presented simple projection algorithms which does not need projection to an auxiliary set Ω.

In this paper, we introduce two projection algorithms for solving the MSFP, inspired by Beck and Teboulle's iterative shrinkage-thresholding algorithm for linear inverse problem [12]. The first algorithm modifies Censor et al.'s method which does not need projection to an auxiliary set Ω. The second algorithm is self-adaptive and adopts the backtracking rule to determine the stepsize. We firstly study the global rate of convergence of two algorithms and prove that the sequences generated by the proposed algorithms converge to a solution of the MSFP. Some numerical results are presented, which illustrate the efficiency of the proposed algorithms.

2. Preliminaries

In this section, we review some definitions and lemmas which will be used in the main results.

The following lemma is not hard to prove (see [1, 13]).

Lemma 1. Let 𝑝 be given as in (3). Then (i)𝑝is convex and continuously differential,(ii)∇𝑝(𝑥) is Lipschitz continuous with 𝐿(𝑝)=∑^𝑡_𝑖=1𝛼_𝑖+𝜌(𝐴^𝑇𝐴)∑^𝑟_𝑗=1𝛽_𝑗 as the Lipschitz constant, where 𝜌(𝐴^𝑇𝐴) is the spectral radius of the matrix 𝐴^𝑇𝐴.

For any 𝜏>0, consider the following quadratic approximation of 𝑝(𝑥) at a given point 𝑦: 𝑅_𝜏𝑦+𝑦+𝜏𝑥,𝑦∶=𝑝𝑥−𝑦,∇𝑝2𝑥−𝑦², (6) which admits a unique minimizer 𝐹_𝜏𝑦∶=argmin𝑅_𝜏𝑥,𝑦∶𝑥∈ℝ^𝑁. (7) Simple algebra shows that (ignoring constant terms in 𝑦) 𝐹_𝜏𝑦=argmin_𝑥𝜏21𝑥−𝑦−𝜏∇𝑝𝑦²1=𝑦−𝜏∇𝑝𝑦. (8)

The following lemma is well known and a fundamental property for a smooth function in the class 𝐶^1,1; for example, see [14, 15].

Lemma 2. Let 𝑓∶ℝ^𝑛→ℝ be a continuously differentiable function with Lipschitz continuous gradient and Lipschitz constant 𝐿(𝑓). Then, for any 𝐿>𝐿(𝑓), 𝑓(𝑥)𝑦+𝑦+𝐿≤𝑓𝑥−𝑦,∇𝑓2𝑥−𝑦²,𝑓𝑜𝑟𝑒𝑣𝑒𝑟𝑦𝑥,𝑦∈ℝ^𝑛. (9)

We are now ready to state and prove the promised key result.

Lemma 3 (see [12]). Let 𝑦∈ℝ^𝑛 and 𝜏>0 be such that 𝑝𝐹_𝜏𝑦≤𝑅_𝜏𝐹_𝜏𝑦,𝑦. (10)
Then for any 𝑥∈ℝ^𝑛, 𝑝(𝑥)−𝑝𝐹_𝜏𝑦≥𝜏2𝐹_𝜏𝑦−𝑦²+𝜏𝑦−𝑥,𝐹_𝜏𝑦−𝑦. (11)

Proof. From (10), we have 𝑝(𝑥)−𝑝𝐹_𝜏𝑦≥𝑝(𝑥)−𝑅_𝜏𝐹_𝜏𝑦,𝑦. (12) Now, from the fact that 𝑝 are convex, it follows that 𝑝(𝑥)≥𝑝𝑦+𝑥−𝑦,∇𝑝𝑦. (13) On the other hand, by the definition of 𝑅_𝜏(𝑥,𝑦), one has 𝑅_𝜏𝐹_𝜏𝑦,𝑦=𝑝𝑦+𝐹_𝜏𝑦𝑦+𝜏−𝑦,∇𝑝2𝐹_𝜏𝑦−𝑦². (14) Therefore, using (12)–(14), it follows that 𝑝(𝑥)−𝑝𝐹_𝜏𝑦𝜏≥−2𝐹_𝜏𝑦−𝑦²+𝑥−𝐹_𝜏𝑦𝑦𝜏,∇𝑝=−2𝐹_𝜏𝑦−𝑦²+𝜏𝑥−𝐹_𝜏𝑦,𝑦−𝐹_𝜏𝑦=𝜏2𝐹_𝜏𝑦−𝑦²+𝜏𝑦−𝑥,𝐹_𝜏𝑦−𝑦, (15) where in the first equality above we used (8).

Remark 4. Note that, from Lemmas 1 and 2, it follows that if 𝜏≥𝐿(𝑝), then the condition (10) is always satisfied for 𝐹_𝜏(𝑦).

3. Two Projection Algorithms

In this section, we propose two projection algorithms which do not need an auxiliary set Ω; one modifies the algorithm introduced by Censor et al. [1] and the other is a self-adaptive algorithm which solves the MSFP without prior knowledge of spectral radius of the matrix 𝐴^𝑇𝐴.

Algorithm 5. Let 𝐿₁≥𝐿(𝑝) be a fixed constant and take 𝜏_𝑛∈(𝐿(𝑝),𝐿₁). Let 𝑥₀ be arbitrary. For 𝑛=0,1,2,…, compute 𝑥_𝑛+1=𝑥_𝑛−1𝜏_𝑛∇𝑝𝑥_𝑛. (16)

Remark 6. Algorithm 5 is different from Censor et al.'s algorithm (5) in [1] and it does not need the projection to an auxiliary simple nonempty closed convex set Ω. In Algorithm 5, we take 𝜏_𝑛>𝐿(𝑝) instead of 𝜏_𝑛>𝐿(𝑝)/2 (as in Censor et al.'s algorithm) which is restricted to a smaller range.

Algorithm 7. Given 𝛾>0 and 𝜂>1, let 𝑥₀ be arbitrary. For 𝑛=0,1,2,…, find the smallest nonnegative integer 𝑚_𝑛 such that 𝜏_𝑛=𝛾𝜂^𝑚𝑛 and 𝑥_𝑛+1=𝑥_𝑛−1𝜏_𝑛∇𝑝𝑥_𝑛, (17) which satisfies 𝑝𝑥_𝑛+1−𝑝𝑥_𝑛+∇𝑝𝑥_𝑛,𝑥_𝑛−𝑥_𝑛+1≤𝜏_𝑛2𝑥_𝑛−𝑥_𝑛+1². (18)

Remark 8. Note that the sequence of function values {𝑝(𝑥_𝑛)} produced by Algorithms 5 and 7 is nonincreasing. Indeed, for every 𝑛≥1, 𝑝𝑥_𝑛+1≤𝑅_𝜏𝑛𝑥_𝑛+1,𝑥_𝑛≤𝑅_𝜏𝑛𝑥_𝑛,𝑥_𝑛=𝑝𝑥_𝑛, (19) where the first inequality comes from Lemma 2 for Algorithm 5 and from (18) for Algorithm 7, and the second inequality follows from (7). 𝜏_𝑛 in (19) is either chosen by the backtracking rule (18) or 𝜏_𝑛∈(𝐿(𝑝),𝐿₁), where 𝐿(𝑝) is a given Lipschitz constant of ∇𝑝.

Lemma 9. There holds𝛽𝐿𝑝≤𝜏_𝑛≤𝛼𝐿𝑝, (20) where 𝛼=𝐿₁/𝐿(𝑝), 𝛽=1 in Algorithm 5 and 𝛼=𝜂, 𝛽=𝛾/𝐿(𝑝) in Algorithm 7.

Proof. It is easy to verify (20) for Algorithm 5. By 𝜂>1 and the choice of 𝜏_𝑛, we get 𝜏_𝑛≥𝛾. From Lemma 2, it follows that inequality (18) is satisfied for 𝜏_𝑛≥𝐿(𝑝), where 𝐿(𝑝) is the Lipschitz constant of ∇𝑝. So, for Algorithm 7 one has 𝜏_𝑛≤𝜂𝐿(𝑝) for every 𝑛≥1.

Remark 10. From Lemma 9, it follows that backtracking rule (18) is well defined.

Remark 11. In algorithm “ISTA with backtracking” proposed by Beck and Teboulle [12], they took 𝜏_𝑛=𝜏_𝑛−1𝜂^𝑚𝑛, with 𝜏₀>0 and 𝜂>1. It is obvious that 𝜏_𝑛 increases with 𝑛. It is verified that small 𝜏_𝑛 is more efficient than a larger one in numerical experiments (see Table 1). So, in Algorithm 7, we take 𝜏_𝑛=𝛾𝜂^𝑚𝑛 for backtracking rule which is smaller than the one in the algorithm of Beck and Teboulle.

Table 1 

Computational results for Example 13 with different algorithms.

Theorem 12. Let {𝑥_𝑛} be a sequence generated by Algorithm 5 or Algorithm 7. Then {𝑥_𝑛} converges to a solution of the MSFP (1), and furthermore for any 𝑛≥1 it holds that 𝑝𝑥_𝑛≤𝛼𝐿𝑝𝑥₀−𝑥^∗22𝑛,∀𝑥^∗∈Γ. (21)

Proof. Invoking Lemma 3 with 𝑥=𝑥^∗, 𝑦=𝑥_𝑘, and 𝜏=𝜏_𝑘, we obtain 2𝜏_𝑘𝑝𝑥^∗−𝑝𝑥_𝑘+1≥𝑥_𝑘+1−𝑥_𝑘²+2𝑥_𝑘−𝑥^∗,𝑥_𝑘+1−𝑥_𝑘=𝑥_𝑘+1−𝑥^∗2−𝑥_𝑘−𝑥^∗2, (22) which combined with (20) and the fact that 𝑝(𝑥^∗)=0, 𝑝(𝑥_𝑘+1)≥0 yields −2𝑝𝑝𝑥𝛼𝐿^𝑘+1≥𝑥_𝑘+1−𝑥^∗2−𝑥_𝑘−𝑥^∗2, (23) which implies 𝑥_𝑘+1−𝑥^∗≤𝑥_𝑘−𝑥^∗. (24) So {𝑥_𝑛} is a Fejér monotone sequence. Summing the inequality (23) over 𝑘=0,1,…,𝑛−1 gives −2𝛼𝐿𝑝^𝑛−1_𝑘=0𝑝𝑥_𝑘+1≥𝑥_𝑛−𝑥^∗2−𝑥₀−𝑥^∗2. (25) Invoking Lemma 3 one more time with 𝑥=𝑦=𝑥_𝑘 and 𝜏=𝜏_𝑘 yields 2𝜏_𝑘𝑝𝑥_𝑘−𝑝𝑥_𝑘+1≥𝑥_𝑘−𝑥_𝑘+1². (26) Since 𝜏_𝑘≥𝛽𝐿(𝑝) (see (20)) and 𝑝(𝑥_𝑘)−𝑝(𝑥_𝑘+1)≥0 (see (19)), it follows that 2𝑝𝑝𝑥𝛽𝐿^𝑘−𝑝𝑥_𝑘+1≥𝑥_𝑘−𝑥_𝑘+1². (27) Multiplying the last inequality by 𝑘 and summing over 𝑘=0,…,𝑛−1, we obtain 2𝛽𝐿𝑝^𝑛−1_𝑘=0𝑘𝑝𝑥_𝑘−(𝑘+1)𝑝𝑥_𝑘+1+𝑝𝑥_𝑘+1≥^𝑛−1_𝑘=0𝑘𝑥_𝑘−𝑥_𝑘+1², (28) which simplifies to 2𝑝𝑥𝛽𝐿−𝑛𝑝_𝑛+^𝑛−1_𝑘=0𝑝𝑥_𝑘+1≥^𝑛−1_𝑘=0𝑘𝑥_𝑘−𝑥_𝑘+1². (29) Adding (25) and (29) times 𝛽/𝛼, we get −2𝑛𝑝𝑝𝑥𝛼𝐿^𝑛≥𝑥_𝑛−𝑥^∗2+𝛽𝛼^𝑛−1_𝑘=0𝑘𝑥_𝑘−𝑥_𝑘+1²−𝑥₀−𝑥^∗2, (30) and hence, it follows that 𝑝𝑥_𝑛≤𝛼𝐿𝑝𝑥₀−𝑥^∗22𝑛,∀𝑥^∗∈Γ, (31) which yields lim_𝑛→∞𝑝𝑥_𝑛=0. (32) Since {𝑥_𝑛} is Fejér monotone, it is bounded. To prove the convergence of {𝑥_𝑛}, it only remains to show that all converging subsequences have the same limit. Suppose in contradiction that two subsequences {𝑥_𝑛𝑘} and {𝑥_𝑛𝑙} converge to different limits 𝑥 and 𝑥, respectively (𝑥≠𝑥). We are to show that 𝑥 is a solution of the MSFP. The continuity of 𝑝(𝑥) then implies that 0≤𝑝(𝑥)=lim_𝑘→∞𝑝𝑥_𝑛𝑘=lim_𝑛→∞𝑝𝑥_𝑛=0. (33) Therefore, 𝑝(𝑥)=0; that is, 𝑥∈𝐶=⋂^𝑡_𝑖=1𝐶_𝑖 and 𝐴𝑥∈𝑄=⋂^𝑟_𝑗=1𝑄_𝑗; that is, 𝑥 is a solution of the MSFP. Similarly, we can show that it is a solution of the MSFP. Now, by Fejér monotonicity of the sequence {𝑥_𝑛}, it follows that the sequence {‖𝑥_𝑛−𝑥‖} is bounded and nonincreasing and thus has a limit lim_𝑛→∞‖𝑥_𝑛−𝑥‖=𝑙₁. However, we also have lim_𝑛→∞‖𝑥_𝑛−𝑥‖=lim_𝑘→∞‖𝑥_𝑛𝑘−𝑥‖=0, and lim_𝑛→∞‖𝑥_𝑛−𝑥‖=lim_𝑙→∞‖𝑥_𝑛𝑙−𝑥‖=‖𝑥−𝑥‖, so that 𝑙₁=0=‖𝑥−𝑥‖, which is obviously a contradiction. Thus {𝑥_𝑛} converges to a solution of the MSFP (1). The proof is completed.

4. Numerical Tests

In order to verify the theoretical assertions, we present some numerical tests in this section. We apply Algorithms 5 and 7 to solve two test problems of [4] (Examples 13 and 14) and compare the numerical results of two algorithms.

For convenience, we denote the vector with all elements 0 by 𝑒₀ and the vector with all elements 1 by 𝑒₁ in what follows. In the numerical results listed in the following tables, “Iter.” and “Sec.” denoted the number of iterations and the CPU time in seconds, respectively. For Algorithm 7, “InIt.” denoted the number of total iterations of finding suitable 𝜏_𝑛 in (18).

Example 13 (see [4]). Consider a split feasibility problem as finding 𝑥∈𝐶={𝑥∈ℝ⁵∣‖𝑥‖≤0.25} such that 𝐴𝑥∈𝑄={𝑦=(𝑦₁,𝑦₂,𝑦₃,𝑦₄)^𝑇∈ℝ⁴∣0.6≤𝑦_𝑗≤1,𝑗=1,2,3,4}, where the matrix 𝐴=2−1323125212021−22−10−35. (34) The weights of 𝑝(𝑥) were set to 𝛼=0.9 and  𝛽=0.1. In the implementation, we took 𝑝(𝑥)<𝜀=10⁻⁹ as the stopping criterion as in [4].

For Algorithm 5, we tested 𝜏_𝑛=1.01𝐿(𝑝),1.1𝐿(𝑝),…,1.9𝐿(𝑝) and the numerical results were reported in Table 1 with different initial points  𝑥₀. (Since the number of iterations for 𝜏_𝑛=1.5𝐿(𝑝),1.6𝐿(𝑝),…,1.9𝐿(𝑝) was larger than those for 𝜏_𝑛≤1.4𝐿(𝑝), we only reported the results for 𝜏_𝑛≤1.4𝐿(𝑝).) We took 𝛾=1 and 𝜂=1.1 for Algorithm 7. Table 1 shows that Algorithm 5 was efficient when choosing a suitable  𝜏_𝑛 (𝜏_𝑛∈(𝐿(𝑝),1.1𝐿(𝑝)) was the best choice for the current example), while the number of iterations of Algorithm 5 was still larger than those for Algorithm 7.

Example 14 (see [4]). Consider the MSFP, where 𝐴=(𝑎_𝑖𝑗)_𝑁×𝑁∈ℝ^𝑁×𝑁 and 𝑎_𝑖𝑗∈(0,1) generated randomly: 𝐶_𝑖=𝑥∈ℝ^𝑁∣𝑥−𝑑_𝑖≤𝑟_𝑖𝑄,𝑖=1,2,…,𝑡,_𝑗=𝑦∈ℝ^𝑁∣𝐿_𝑗≤𝑦≤𝑈_𝑗,𝑗=1,2,…,𝑟, (35) where 𝑑_𝑖∈ℝ^𝑁 is the center of the ball 𝐶_𝑖, 𝑒₀≤𝑑_𝑖≤10𝑒₁, and 𝑟_𝑖∈(40,50) is the radius; 𝑑_𝑖 and 𝑟_𝑖 are both generated randomly. 𝐿_𝑗 and 𝑈_𝑗 are the boundary of the box 𝑄_𝑗 and are also generated randomly, satisfying 20𝑒₁≤𝐿_𝑗≤30𝑒₁, 40𝑒₁≤𝑈_𝑗≤80𝑒₁, respectively. The weights of 𝑝(𝑥) were 1/(𝑡+𝑟). The stopping criterion was 𝑝(𝑥)<𝜀=10⁻⁴ with the initial point 𝑥₀=𝑒₀∈ℝ^𝑁.

We tested Algorithms 5 and 7 with different 𝑡 and 𝑟 in different dimensional Euclidean space. In Algorithm 5, since a smaller 𝜏_𝑛 is more efficient than a larger one, we take 𝜏_𝑛=1.01𝐿(𝑝) in the experiment. We take 𝛾=1 and  𝜂=1.2 for Algorithm 7. For comparison, the same random values were taken in each test. The numerical results were listed in Table 2, from which we can observe the efficiency of the self-adaptive Algorithm 7, both from the points of view of number of iterations and CPU time.

Table 2 

Computational results for Example 14 with different dimensions and different numbers of 𝐶_𝑖 and 𝑄_𝑗.

Conflicts of Interests

There is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors express their thanks to Dr. Wenxing Zhang for his help in numerical tests. This research was supported by the National Natural Science Foundation of China (no. 11201476) and Fundamental Research Funds for the Central Universities (no. 3122013D017).